Bioinformatic analysis of the coding region of the melatonin receptor 1b gene as a reliable DNA marker to resolve interspecific mammal phylogenetic relationships

1.   Introduction

The problem of living beings taxonomy that would adequately reflect their phylogenetic history has always been at the forefront for evolutionary biologists and taxonomists due to the controversial points that often arise in the classification process. Resolving such points is one of the tasks of the modern phylogenetics. Until recently, the construction of classification systems was based on the comparison of morpho-anatomical features or karyological data but the high level of homoplasy, the variety of morphological characteristics, and numerous chromosomal variations make the work with such features more difficult. Nowadays, taxonomic research is more often based on the concept of "molecular clocks", the use of molecular markers, and the phylogenetic trees construction, providing an opportunity to clarify or even significantly supplement traditional classifications of living organisms ^[1].

1.1.   Phylogenetics and the molecular clock

The presence of a certain number of differences in the amino acid and nucleotide sequences of hemoglobin in different animal species was first noted while studying these sequences in the middle of the 20^thcentury. It was found that hemoglobin molecules of humans and horses differ by 18 amino acids, mice and horses by 22 amino acids, humans and sharks by 79 amino acid residues. Subsequently, similar differences were found as well for the nucleotide sequence encoding this protein ^[2].

This observation led to the idea that each stage of evolution is marked by a certain number of changes in the primary gene or protein sequences. Zuckerkandl & Pauling proposed the "molecular clock" theory trying to explain these unusual results. According to this theory, the frequency of amino acid and nucleotide substitutions correlates with the timeline of evolution ^[2]. Certain genes and proteins can play the role of "small" and "big" hands of the "molecular clock" covering the desired part of the "evolutionary dial" depending on the sequence substitution rate. Since then, biomolecules have become not only objects but also full-fledged tools of molecular-genetic and taxonomic research.

However, large-scale taxonomic studies based on a single gene or protein alone can lead to problems in interpreting the results as the frequency of substitutions in different sequences always varies to some extent.

Moreover, from a phylogenetic point of view the nucleotide sequence is much more informative than the amino acid one primarily due to the possibility of synonymous mutations affecting DNA but having no effect on the encoded amino acid sequence ^[3] (see Figure 1).

Figure 1.  From phylogenetic point of view DNA contains more information than protein. The figure shows a schematic representation of a mutational process where five substitutions occurred in a part of DNA that did not affect the coded sequence of a protein due to the degeneracy of the genetic code.

DownLoad: Full-Size Img PowerPoint

Effective phylogenetic DNA markers for systematization members of the animal kingdom can include mitochondrial and nuclear genes as well as non-coding sequences having the different limits of resolution in molecular phylogenetic analysis.

1.2.   Variety and the limits of resolution of DNA markers application for animal kingdom systematics

DNA markers can be conventionally divided into mitochondrial and nuclear genetic markers ^[2]. Mitochondrial DNA (mtDNA) sequences that have proven to be effective markers for phylogeny ^[4] include the D-loop ^[5], cytochrome b (Cytb) gene ^[2], 12S and 16S rRNA genes ^[5], some regions of the cytochrome oxidase C gene ^[2], and sequences encoding NADH dehydrogenase (ND) subunits ^[5].

The mitochondrial genome is appealing to specialists in the gene systematics field primarily because of its small size and the possibility to select universal primers for amplification that considerably facilitates the analysis. The GenBank database contains partial or even complete mtDNA sequences from almost all modern animal species ^[2]. Nowadays, the exclusivity of mtDNA as a source of universal phylogenetic markers is questioned more and more often. This is due to the high rate of nucleotide substitutions which is higher than in nuclear DNA (nDNA) ^[4], as well as the high level of "silent substitutions" resulting from synonymous mutations ^[3,5]. It is known that the average rate of nucleotide substitutions in mtDNA is sometimes 10 times higher than in nDNA, and is estimated at 2–4% per one million years ^[5,6]. At the same time, "silent substitutions" appear 100 times more frequently in the mitochondrial genome ^[7].

Until recently, the high frequency of substitutions in mtDNA was seen as an advantage for phylogenetic studies ^[2,4] but studies in recent years suggest otherwise. The high frequency of substitutions in the mitochondrial genome enhances the "saturation effect" in the nucleotide sequence shielding the "phylogenetic" signal ^[5]. Some mitochondrial genes reach a "saturation" plateau after 15–20 million years at just 25% variation ^[5,8]. Further on, substitutions are concentrated at positions that were once already substituted, which is a potential source of homoplasy ^[5]. As a result, the mtDNA of species that diverged 60–80 million years ago may not differ from species that separated from a common ancestor no more than 20–30 million years ago ^[5].

Because of the slower substitution rate, site saturation in nDNA is slower than in mtDNA, which is why it is more and more often that different regions of nDNA are used to reconstruct the evolution of taxa of different orders. This applies both to individual areas and to their combinations ^[5]. Nowadays, both coding and non-coding regions of the nDNA are used in the reconstruction of phylogenies such as nucleolus organizer regions (NORs) ^[9], alpha-2B adrenergic receptor gene (A2AB) ^[5,10], von Willebrand factor gene (vWF) ^[5,10], interphotoreceptor retinoid-binding protein gene (IRBP) ^[11], microsatellite repeats ^[12,13], introns ^[5], etc. Phylogenetic reconstructions with nuclear and mitochondrial DNA regions served as a basis for the compilation of various tables and nomograms of DNA marker boundaries in molecular phylogenetic analysis ^[5] (see Figure 2).

Figure 2-A.  Limits of resolution of mitochondrial DNA markers application in molecular phylogenetic analysis (Sup - superorder, O - order, F - family, G - genus, Sp - species, P - population, Sub - subfamily, I - individual).

DownLoad: Full-Size Img PowerPoint

Figure 2-B.  Limits of resolution of nuclear DNA markers application in molecular phylogenetic analysis (Sup - superorder, O - order, F - family, G - genus, Sp - species, P - population, Sub - subfamily, I - individual). The use of multiple satellite repeats as part of a multilocus analysis increases the resolution of satellite nDNA up to the level of the individual.

DownLoad: Full-Size Img PowerPoint

The marker regions of nDNA, just like other types of markers, have a number of shortcomings that are largely reflected in the phylogenies reconstructed on their basis. Non-coding regions may have excessive mutational variability, whereas coding genes, while changing more slowly, may be subject to convergent evolution. As a result, the number of informative markers that can reliably reflect the evolutionary process is not so big. We can conclude that there is a need to expand the range of sequences that could add to the list of effective phylogenetic markers and thereby increase the reliability of evolutionary history reconstruction, which is partly the aim of this work. One of the criteria for the new marker search was the highest conservatism of the target nucleotide sequence in most representatives of the animal kingdom as well as in those processes in which the sequence or its protein product is involved. In process of reviewing the relevant scientific literature and bioinformatic sequences annotated in databases it became obvious that the melatonin receptor family satisfies the above-mentioned requirements.

1.3.   The melatonin receptors

Melatonin (N-[2-(5-methoxyindol-3-yl)ethyl]acetamide) is a unique indolamine found in all taxa of living organisms from bacteria to higher vertebrates. It is thought to be one of the oldest signaling molecules ^[14,15]. Melatonin mediates the regulation of a variety of biological functions in animals through three different subtypes of membrane G-protein-coupled receptors: Mel1a (MT1, MTNR1A) ^[16], Mel1b (MT2, MTNR1B) ^[17] and Mel1c (MTNR1C) ^[18]. MTNR1A and MTNR1B receptors are expressed in various mammalian tissues, whereas additional receptors MTNR1C have only been identified in fish, amphibia, and some birds ^[18,19]. Interestingly enough, peak MTNR1C expression occurs during the day ^[20] rather than night in contrast to MTNR1A ^[21,22] and MTNR1B ^[23,24]. In a cell nucleus melatonin has an affinity to retinoic acid-related orphan receptors - ROR 1α, ROR 2α, RZR, and through these proteins it is able to directly modify the transcriptional activity of many genes ^[25,26].

Melatonin and its receptors are involved in numerous physiological and pathophysiological processes such as regulation of blood pressure ^[27], circadian rhythms ^[28], retinal function ^[29,30], oncogenesis ^[31], reproductive seasonality ^[32], ovulation ^[33,34], and regulation of stem cell proliferation ^[35,36] and differentiation ^[37,38].

The MTNR1B gene in representatives of different families of living creatures is contained in different autosomes. For example, in humans and higher primates this gene is located on the long (q) arm of chromosome 11, while in representatives of the Felidae family it is located on chromosome D1. In all mammals, this nucleotide sequence encodes one of the most affine forms of melatonin receptor MT2 ^[39]. The product of this gene is an integral G-protein coupled transmembrane receptor. It is believed that namely MT2 receptors are involved in the regulation of seasonal and circadian rhythms ^[40]. However, recent studies have revealed the involvement of melatonin receptors not only in the circadian system but also in the homeostatic regulation of energy balance, which seems to be related to the influence of melatonin-MTNR1B complex on the transmission of insulin and leptin hormonal signals in hypothalamus. In other words, MTNR1B is an important modulator of some intracellular signal transduction, affecting lipid and glucose metabolism ^[41,42].

MTNR1B and MTNR1A also seem to be convenient markers as they occur in representatives of almost all taxonomic groups of eukaryotic organisms and have a fairly low level of genetic polymorphism, possibly due to the archaism and invariance of circadian processes in various representatives of the animal world. Due to this invariance, even single mutations in these genes lead to the development of pathological conditions such as type 2 diabetes ^[43,44] and autism ^[45].

2.   Materials and methods

2.1.   Sequences

Review of literature data shows that most melatonin receptors have in common the presence of seven obligate transmembrane domains in their protein molecule structure. This structural stability in different living organisms shows that relevant amino acid and/or nucleotide sequences are potential markers for phylogenetic studies. However, the review of relevant literature and our own studies have shown that not all types of melatonin receptor are suitable for such purposes. For example, the article by Li, et.al. ^[20] shows that MTNR1C receptors are not specific to mammals at all and for this reason cannot fully reflect the dynamics of the real evolutionary process. Data from the same article as well as the results of our preliminary studies (not shown) showed that MTNR1A cannot be effective phylogenetic markers as well due to its low statistical support and the apparent divergence of phylogenetic reconstructions based on them from established ones ^[20].

As a result, it became clear that MTNR1B receptors were the most suitable for phylogeny construction, although not without pitfalls. Preliminary dendrograms based on the amino acid sequences of these receptors (not shown) were characterized by very low statistical support (bootstrap index < < 50%) and frequent occurrence of nonsense branches, which indirectly indicated that synonymous substitutions were present in exons and that amino acid sequences of this receptor cannot be satisfactory phylogenetic markers, unlike nucleotide ones.

Based on the literature ^[46,47], the coding sequence of the cytochrome b gene (CDS cyt-b) was chosen for the present study as a control phylogenetic marker as it has proven itself in molecular phylogenetic analysis thanks to its broad taxonomic range enough to comprise different potential indicator taxa (Figure 2-A).

Thus, reference coding nucleotide sequences mtnr1b and cyt-b were selected for 34 different mammalian species (with monotremes as an outgroup) to accomplish the task. All sequence information is annotated in the GenBank database hosted by NCBI (http://www.ncbi.nlm.nih.gov/). A complete list of the sequence identifiers involved in this work is given in Table 1.

Table 1.  Nucleotide sequence identifiers of mtnr1b and cyt-b of different species annotated in GenBank.

	Species	Nucleotide sequence mtnr1b ID	Nucleotide sequence сyt-b ID
Marsupials and placentals	Bos taurus	NM_001206907.1	NC_006853.1
	Bubalus bubalis	XM_006053939.1	NC_006295.1
	Camelus ferus	XM_006182375.1	NC_009629.2
	Canis lupus	XM_844629.2	NC_002008.4
	Ceratotherium simum	XM_004427448.1	NC_001808.1
	Eptesicus fuscus / serotinus	XM_008149162.1	NC_022474.1
	Equus caballus	XM_001917051.1	NC_001640.1
	Equus przewalskii	XM_008515596.1	NC_024030.1
	Felis catus	XM_003992620.2	NC_001700.1
	Galeopterus variegatus	XM_008591633.1	NC_004031.1
	Gorilla gorilla	XM_004051965.1	NC_001645.1
	Homo sapiens	NM_005959.3	NC_011137.1
	Leptonychotes weddellii	XM_006732963.1	NC_008424.1
	Lipotes vexillifer	XM_007449662.1	NC_007629.1
	Loxodonta africana	XM_003415638.1	NC_000934.1
	Macaca mulatta	XM_001084265.2	NC_005943.1
	Microtus ochrogaster / levis	XM_005347416.1	NC_008064.1
	Monodelphis domestica	XM_001369486.1	NC_006299.1
	Myotis davidii	XM_006771944.1	NC_025568.1
	Mus musculus	NM_145712.2	NC_006915.1
	Mustela putorius	XM_004763604.1	NC_020638.1
	Odobenus rosmarus divergens	XM_004410754.1	NC_004029.2
	Pan paniscus	XM_003813777.1	NC_001644.1
	Pan troglodytes	XM_522146.5	NC_001643.1
	Panthera tigris	XM_007079066.1	NC_010642.1
	Papio anubis	XM_003910543.2	NC_020006.2
	Physeter catodon	XM_007107131.1	NC_002503.2
	Pongo abelii	XM_002822349.1	NC_002083.1
	Pteropus alecto	XM_006907469.1	NC_023122.1
	Sarcophilus harrisii	XM_003764327.1	NC_018788.1
	Sus scrofa	XM_003129761.2	NC_000845.1
	Trichechus manatus	XM_004391269.1	NC_010302.1
	Vicugna pacos	XM_006206752.1	NC_002504.1
Outgroup monotremes	Ornithorhynchus anatinus	XM_003430777.1	NC_000891.1

---

 | Show Table

DownLoad: CSV

2.2.   Multiple and pairwise alignments

Multiple and pairwise sequence alignment was performed using the heuristic algorithm ClustalW integrated into the MEGA 6.0 computer software package for molecular genetic and evolutionary analysis ^[48,49]. Such a method consists in a vector alignment procedure, components of which are two values - the total weight of the matches and the total number of deleted fragments. This algorithm requires setting a number of parameters, such as GEP (Gap Extension Penalty), GOP (Gap Opening Penalty), delay divergent sequences which ensure priority alignment of more similar DNA sites, and transition weight or fraction (A↔G or C↔T). In the present work, the GEP and GOP parameters for DNA were set to 6.66 and 15, respectively. At the same time, a delay divergent sequences rate was 30% and a transition weight setting at a mean of 0.5 ^[49].

2.3.   Phylogenetic analysis

Phylogenetic analysis was performed using the MEGA 6.0 software ^[50]. Phylogenetic trees were constructed using minimal evolution (ME), neighbor-joining (NJ), and maximum likelihood (ML) methods ^[51,52,53]. The bootstrap method ^[54] was used to assess the statistical validity of the groupings obtained by ME, NJ and ML methods. Significance of reconstructions was estimated using the bootstrap method using 1000 replicates. The dendrograms show the percentages of support for significant branches. Support scores below 50% are considered statistically insignificant and are usually not shown on the reconstructed tree. The relative divergence time was calculated automatically in MEGA 6.0.

3.   Results

3.1.   Phylogenetic analysis of the Cyt-b control DNA sequence

The coding sequence of the Cyt-b gene averaged 1143 bp for all species taken in the study, ranging from 1137 bp for Loxodonta africana to 1149 bp for Monodelphis domestica. 60% of the positions within 1143 bp were informative for phylogenetic analysis.

Multiple and pairwise sequence alignment for all 34 species provided a matrix of 1149 sites with an insignificant end GAP that cannot affect further analysis and therefore does not require elimination. The frequencies of each of the 4 nucleotides (a, t, g, c) vary slightly among different species (see Table 2.). Thus, the calculated ratio of frequency of transitions to transversions (see the Table 2.) for purine (k1) and pyrimidine bases (k2) has allowed to calculate the general average transition/transversion ratio (R), which was equal to 0.397, where

Table 2.  Frequency of occurrence of each of the 4 nucleotides (a, t, g, c), ratio of frequencies of transitions to transversions for purines (k1) and pyrimidines (k2). Lnl parameter (maximum log likelihood). Mean evolutionary distance (K) and mean frequency of evolutionary substitutions (λ).

The sequence	a%	c %	g%	t%	k1	k2	Lnl	K	λ (10-9)
CDS mtnr1b	17.68	31.45	25.37	25.50	5.19	3.57	-7772.99	0.200	0.43
CDS cyt-b	29.57	30.52	12.86	27.05	2.53	0.01	-18599.07	0.311	0.62

---

 | Show Table

DownLoad: CSV

Based on the distance matrices computed with MEGA 6.0, the evolutionary distance mean (K) and the frequency of evolutionary substitutions (λ) for the cytochrome b gene were calculated to be 0.311 units and 0.62 x 10^-9 nucleotide substitution per site per year, respectively, where

Equally important was the definition of the natural logarithm of the probability function

where p is the probability and the binomial coefficient gives the number of ways to order h successes out of n trials. The calculation of the logarithm is available in the MEGA 6.0 software. This logarithm is commonly indicated as the maximum log likelihood (Lnl) and has a predictive value, indicating the degree of validity of the phylogenetic model build on the basis of a marker. The higher the Lnl parameter (the less negative it is) - the more adequate is the proposed model (see Table 2.).

Phylogenetic trees constructed from cytochrome b gene sequences have similar topology when the NJ, ML, and ME methods are used. The final results of phylogenetic analysis are illustrated by phylogenetic trees constructed using the ML method (Figure 3-4). At least 8 distinct clusters are identified on the tree (Figure 3B). Most of the identified clusters have satisfactory support but their composition and mutual arrangement at the placental level are ambiguous.

Figure 3.  Phylograms constructed by ML analysis of mtnr1b (A) and cyt-b (B). Monotremes were taken as outgroup; red asterisks indicate nodes with branch support < 50% and topology instabilities. Branch lengths are proportional to the number of nucleotide substitutions. Scale bar refers to a phylogenetic distance of 0.05 nucleotide substitutions per site. Numbers on the branches indicate bootstrap percentage after 1000 replications in constructing the tree.

DownLoad: Full-Size Img PowerPoint

Figure 4.  Condensed cladograms obtained from the phylogenetic analysis of CDS mtnr1b (A) and cyt-b (B) using ME, NJ, and ML methods. Monotremes were taken as outgroup; numbers above the nodes correspond to bootstrap percentage after 1000 replications for ML method, numbers below the nodes correspond to bootstrap indices for ME, NJ, respectively.

DownLoad: Full-Size Img PowerPoint

The carnivores, primates, rodents, and marsupials form their own well-supported and satisfactorily supported clades. The Cetartiodactyla group is satisfactorily supported only in the case of Bovidae and Cetacea, while evolutionary branches of the Suiformes and Tylopods have an unstable position, which is why the consensus cladogram (Figure 4B) identifies them as separate clusters. The greatest distance separates the platypus, as a member of the monotremes, from all other mammalian species (up to 0.44). The smallest genetic distances were found between E. cabalus and E. przewalskii (0.001), P. paniscus and P. troglodytes (0.05). The distance separating P. paniscus and H. sapiens, H. sapiens and P. troglodites, P. paniscus and G. gorilla, and P. troglodytes and G. gorilla is also relatively small (≤ 0.12).

3.2.   Phylogenetic analysis of MTNR1B DNA sequences

The coding sequence of the MTNR1B gene averaged 1100 bp, ranging from 1014 bp for Ornithorhynchus anatinus to 1380 bp for Galeopterus variegatus. Multiple and pairwise sequence alignment for all 34 species provided a matrix of 1476 sites, 357 of which were flanking GAPs (286 positions at the beginning and 71 at the end of the matrix). After removal of the terminal GAPs, 1119 sites were available for analysis, about 60% of which were informative for phylogenetic analysis, just as in the case of cytochrome.

The calculated values of the frequency of occurrence of each of the 4 nucleotides, the ratio of frequencies of transitions to transversions (k_1-2), the average value of evolutionary distance (K), the frequency of evolutionary substitutions (λ), and the Lnl parameter are summarized in Table 2.

Based on the data on frequencies of transitions to transversions (k_1-2), the overall ratio of transitions to transversions (R) was calculated, which did not exceed a value of 2.12 for mtnr1b.

The phylogenetic analysis is generally in agreement with the obtained results using cytochrome b sequence: trees built from the CDS mtnr1b sequences have identical topology when using the NJ, ML, and ME methods. At least 8 separate clusters are also identified on that tree (Figure 3A). Most of the selected clusters have high statistical support. Their composition and mutual arrangement are more stable compared to controls, but some topological uncertainties are still present at the level of such superorders as Laurasiatheria and Euarchonta.

Cetartiodactyla, odd-toed ungulates, carnivores, primates, rodents, afrotherians, and marsupials form their own well-supported clades. In this case, the Cetartiodactyla group finds strong support by grouping the ruminants, cetaceans, suiformes and tylopods into a single cluster (Figure 3A). Bats (chiroptera) are often grouped with the order Carnivora, but due to low bootstrap support (<50) they are given as paraphyletic groups in the consensus cladogram (Figure 4A). As in the case of the calculations for the control marker cyt-b, the greatest distance separates the monotremes from all other animal species (up to 0.56). The negligible genetic distances are found between E. cabalus and E. przewalskii (0.001), P. paniscus and P. troglodytes (0.005), H. sapiens and G. gorilla (0.009), P. paniscus and H. sapiens (0.011), H. sapiens and P. troglodytes (0.014), P. paniscus and G. gorilla (0.015), P. troglodytes and G. gorilla (0.018), P. tigris and F. catus (0.018), B. taurus and B. bubalis (0.018), M. mulatta and H. sapiens (0.032), P. anubis and P. paniscus (0.033), P. anubis and P. troglodytes (0.036), M. mulatta and G. gorilla (0.036), P. abelii and P. anubis (0.039), C. simum and E. przewalskii (0.075), C. simum and E. cabalus (0.077).

4.   Conclusions

For this work we selected the coding sequences of the melatonin receptor 1B gene (MTNR1B) distinctive to placental mammals, marsupials and monotremes. Multiple alignments of the protein-coding sequences of the selected gene for all 34 species included an average of 1100 positions containing both homologous regions and regions with deletions and nucleotide substitutions. The alignment showed that CDS mtnr1b combines both highly conserved regions and phylogenetically informative variable regions.

The coding sequence of the cytochrome b gene from the same biological sources was used as a classic phylogenetic marker for comparison. In both cases, the monotremes class was used as the outgroup. It has been experimentally found that it is reasonable to optimize full-length mtnr1b CDSs to determine evolutionary relationships. Optimization should be reduced to the exclusion of flanking GAPs which occur after alignment. Deeper optimization, which truncates mtnr1b to single regions encoding variants of the receptor transmembrane domain, significantly reduces the quality of the constructs, leading to nonsense branches (not shown).

Distance matrices generated with the program MEGA 6.0 were used to calculate the average frequency of evolutionary substitutions (λ) for mtnr1b and cyt-b, reflecting the magnitude of both genes' polymorphisms. In the former case, the polymorphism frequency was 1.4 times lower (see Table 2.).

Phylogenetic trees constructed based on mtnr1b using different methods (ML, NJ, ME) showed a similar topology (Figures 3–4), confirming the phylogenetic informativity of the selected nucleotide sequence.

Comparison of cyt-b and mtnr1b as DNA markers for phylogeny shows the superiority of mtnr1b in terms of topology stability, bootstrap support indices, and infraclass level resolution (Figures 3–4).

Thus, mtnr1b provides more confident clustering and better phylogenetic signal to resolve nodes at the level of superorder categories (placental mammals, marsupials, monotremes) compared to the control gene.

Chronograms and condensed consensus cladograms with satisfactory degrees of bootstrap support were generated based on mtnr1b phylogenies, clustering all mammalian species involved in the present study into nine different monophyletic groups (clades) that had evolved from earlier forms during the geological history of Earth (Figures 4–5): 1) Cetartiodactyla, 2) Perissodactyla, 3) Carnivora, 4) Chiroptera, 5) Paenungulata (Afrotheria), 6) Euarchonta, 7) Glires (Rodentia), 8) Marsupialia, 9) Monotremata.

Figure 5.  Consensus chronogram of mammalian phylogenetic relationships with time scale based on comparative analysis of CDS mtnr1b. The vertical branch along which there is instability in the topology of the nodes of different clades - the "line of instability"- is depicted by red color.

DownLoad: Full-Size Img PowerPoint

It is evident from the obtained dendrograms (both experimental and control) that it is the few node-based clades of some placental mammals that have the lowest statistical support (<50%), while statistical support to nodes at other taxonomic levels is higher and more stable (see Figure 3). The works of other authors, who have taken other markers as a basis for constructing mammalian phylogenies, also show low stability of the positions of nodes uniting some clades of placentals and high stability of topologies of other levels ^[55,56].

Such topological instability should not be negatively interpreted, as it itself may be a unique marker reflecting the pattern of speciation that prevailed at various stages of evolution. In order to explain this result, it is worth recollecting the different nature of the process of divergent evolution ^[57,58,59]. We can generalize that divergence can be competitive and non-competitive, depending on the root causes, ultimately affecting the topology of monophyletic groups in the building a phylogenetic tree. While competitive divergence is characterised by a high diversity of species and populations, limited number of ecological niches, and slow evolutionary processes, non-competitive divergence is characterised by a critically low number of species, large number of distinct ecological niches, and high intensity of the evolutionary processes.

Most probably non-competitive divergence is closely related to sympatry ^[60] and may result from global catastrophes (climate change, asteroid impact, etc.) leading to the mass death or extinction of most living beings and the total release of ecological niches. Members of the remaining small number of species or even populations no longer compete with one another but seek to occupy the vacated ecological niches, even in the case of their partial overlap. In this case, processes such as food specialization, population fluctuations, and fixation of spontaneous mutations occur much faster than in competitive evolution, in which intraspecific and inter-population struggle is reduced to just the opposite - trying to stay within a particular niche by all means.

Hypothetically, these two types of the divergent evolution can have an extremely strong influence on the statistical support indicators of phylogenetic tree topologies. The instability of the clade topology of some placental mammals in the dendrograms constructed in this work (Figures 3–5) suggests that mammalian divergence could have been non-competitive around 70 million years ago. Such a thesis seems logical, because it is known that the so-called Cretaceous–Tertiary (K–T) mass extinction of species occurred precisely in this period ^[61,62].

The concepts of competitive and non-competitive divergence can also be closely linked to the phenomenon of gene drift and the bottleneck effect. Assuming that the populations of most living beings, including mammals, were reduced to a critical minimum by forces unrelated to natural selection (such as asteroid impact, tectonic activity or pandemics), the result should have been the accidental elimination of large numbers of individuals and a drastic decrease in the genetic diversity, which is the basis for the bottleneck effect. When populations began to flourish again, they retained some genes that were specific to the surviving individuals and not to the original species generally.

In order to assess whether the divergence on clades in the constructed cladograms is a random or non-random process, we can use Tajima's statistical test ^{[63,64,65,66]}. The calculation of this function is available as part of the MEGA 6.0 software package. We tested unstable nodes and nodal hoards lying before and after the line with node instability. In all cases, monotremes were used as the external group. As a result, Tajima's D-statistics values for nodes lying before and after the "line of instability" were strongly positive (D > > 0), while nodes lying on the "line of instability" had weakly negative values (D < 0).

If the standard interpretation for D values is taken into account, the results may indicate that divergence in the nodes before and after the "line of instability" was of a balanced competitive selection, while divergence in the nodes lying on the line itself may be associated with a sharp non-competitive populations rise after factors, previously limiting their growth, were randomly removed, i.e. there are signs of passing through a bottle-neck.

Judging by the number of nodule clades (see Figure 5), in that epoch there were only a few primitive mammal populations remained on Earth, which had begun developing new ecological niches that had been devastated by the disappearance of competing species of mammals and reptiles. In analogous niches, of which there were a great many at that time, evolution may well have been parallel with increasing frequencies of parallelisms and backward substitutions, so that some similar traits appeared independently without preventing these groups from "mixing" at a certain evolutionary stage. As a result, very often we encounter increasing levels of "phylogenetic noise" and topological uncertainty when reconstructing phylogenies at the level of nodes uniting clades of placental mammals, which is not observed at higher (infraclass) and lower (order, species) levels (Figure 4–5).

The analyses presented in the paper demonstrated the applicability of protein-coding regions of the MTNR1B gene as a molecular marker for the phylogenetic study of the evolutionary relatedness among various groups of organisms. The trees constructed in this work, with a few exceptions, are consistent with those established by various molecular markers ^[67,68,69] and paleontological or morphological comparisons ^[70,71,72], indicating great potential of the chosen marker to resolve relationships of long branches of highly divergent Mammalian clades.

Acknowledgments

Funding: The work was performed within the framework of the Program for Basic Research in the Russian Federation for a long-term period (2021-2030) [No. 122022800499-5].

Conflicts of interest

The authors declare no conflict of interest.