Estimation of the congruence between morphogenetic and molecular-genetic modules of gray voles Microtus S.L. variability along a climatic gradient

Full Text

INTRODUCTION

At present, the modular organization principle of multicellular organisms has a solid conceptual and empirical basis [1–8]. The main criterion for identifying modules is strong consistency of parts inside a module in combination with modules relative independence from each other [9–12]. It enables modular systems not only to function, but evolve in relative autonomy [13, 2, 14–16]. In broad terms, e.g., anatomic units (parts of body and bones), morphometric measurements and also other traits [8] may be considered as constituents (axes) of phenotypic modules.

The attempts to reveal the complexes of integrity-forming phenotypic traits have been taken since the beginning of the last century. The main prerequisite of such investigations consisted in that traits related genetically, ontogenetically or functionally correlate closer than others. Accordingly, the analysis was focused on the matrix of correlation between traits — the method of correlation pleiades [17], which determined the use of the principal components and other latent variables as a way of revealing stable blocks of phenotypic variability [10, 11, 18].

The internal consistency of parts within a module, combined with modules independence from each other, allows modular systems to function and evolve, to some extent, autonomously [2, 9–11, 13–16]. The integrity of morphogenesis imposes some limits on modules excess autonomization. Therefore, the “modularity concept” is closely connected with the concept of “morphological integration” [19, 20]. The traits integrated after E. Olson and R. Miller are such subsets of morphological traits which have a tendency to change jointly during ontogenesis and evolution, opposed to others far less related to each other. The division into subsets is made either formally on the correlation between the traits (r-group), or in the content — on the principle of traits joint functioning (F-group) [19]. Ideally, both divisions should correspond to each other.

One of the consequences of morphological integration is that the variability of organisms is concentrated in some limited directions (more often, in one or two) forming a “space of abilities”, whereas, in other directions, the variability is less [21–24]. It was shown in the work with the quantitative traits of rodent [23] and insectivorous [24] skulls that the first two directions of variability form the size–shape plane. The analysis of congruence with molecular-genetic traits (Cytb) allowed revealing the integrated blocks of two trait systems — morphological and molecular-genetic. Thereat, the quantitative estimates of the congruence between them became considerably higher at the elimination of the size component of variability.

The exponential growth of the research concerning the role of morphological modularity and integration in evolution [3, 6, 8, 10, 12, 25–30] has been observed since the beginning of the century and up to the present.

The presence of the corresponding genetic base is a necessary condition of modular evolution. The genetic base of integration suggests that a selection on some separate structural or functional peculiarities of an organism will always lead to a correlated response in its other parts, whereas the genetic base of modularity will provide a relative independence for evolutionary changes of different parts of an organism [2, 10, 13, 31]. It was particularly noted that inter-level integrations [10] should manifest themselves first in the congruent variability of different structural level modules.

Most of researchers accept that the integration concept is applicable at different biological levels [6, 10] including genetic and ecological integrations [13, 32], left and right sides integration of organisms body [32–37] and evolutionary integration of taxa [38–43].

For the last decades, more perspective methods, in particular, multidimensional scaling and PLS-methods that enable finding out deep joint latent variables belonging to different structural-functional systems and are presented by different types of traits, have become available. This gives an additional opportunity for the inter-level research of integration and modularity, which has just begun and the potential of which is still far from being exhausted [10].

In our work, the search of morphogenetic and molecular-genetic variability modules has been realized, and the quantitative estimation of the congruence between them in ten vole species for the group of genera Microtus in gradient geo-climatic conditions has been made. The choice of trait systems is conditioned by the data on the morphological variability m1 of ten gray vole species Microtus s.l. from the right and left lower jaw sides and the completeness of GenBank database on the mtDNA Cytb gene related to the investigated species. The working hypothesis consisted in that the degree of environmental extremeness will simultaneously affect both the stability of specific morphogeneses of investigated species, and the mtDNA Cytb gene properties.

MATERIALS AND METHODS

The morphotypic variability of the first lower molar (m1) of ten gray vole species Microtus s.l from the Siberian Zoological Museum of RAS, Zoological Museum of MSU and the Zoological Institute of RAS was studied. The description methods of morphotypic variability were published earlier [44]. It is known that the morphotypes from the right and left jaw sides may be different. As the morphotypes are connected with each other by chains of elementary transitions, the longer minimal chain is, the more different the morphogenetic trajectories of the morphotypes being at the ends of this chain are and, accordingly, the more seldom they can be realized on different sides of one individual [44–47].

In total, 5347 pairs of vole molar (m1) chewing surface contours were observed from the right and left lower jaw sides. Fifty-six morphotypes and their 222 pair combinations (direct and reverse combinations were considered different) were found. Twenty-five morphotypes were encountered (left and right, in total) less than five times and were excluded from the further analysis together with their pairs. Accordingly, 5306 chewing surface contours were left, the number of morphotypes being equal to 31, and the number of different combinations — 187 : 30 — symmetric and 157 — asymmetric. Symmetric combinations were encountered more often (Table 1). The fraction of symmetric and asymmetric pairs, common for all species, equaled 72.3 % and 27.7 %, respectively.

 

Table 1

Frequencies number of symmetric pairs of morphotypes m1 for each species and total for all species

Species	4H5	5H4	5H5	5H6	6H5	1K4	1K5	1K6	2K4	2K5	2K6	2K7	3K5	3K6	4K5	4K6	1M3	2M3	3M3	1M4	2M4	3M4	3M5	4M3	4M4	4M5	5M4	5M5	3P4	2Ma2	Sum
Alexandromys fortis	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	144	22	0	123	64	3	0	2	0	0	0	0	0	359
A. maximowiczii	0	0	1	0	0	0	3	0	0	8	1	0	45	8	4	5	0	20	166	0	4	100	10	3	18	2	0	0	0	0	398
A. middendorffii	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	65	0	0	113	16	2	38	1	0	0	0	0	237
A. mongolicus	0	0	0	0	0	0	2	0	0	6	3	0	6	9	0	0	0	4	15	0	14	67	1	0	0	0	0	0	0	0	127
A. oeconomus	1	0	0	0	0	24	828	56	20	819	51	1	16	4	0	0	0	30	2	0	3	0	0	0	0	0	0	0	0	0	1855
Lasiopodomys gregalis	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	4	7	68	2	14	213	6	0	1	0	0	0	1	1	317
Microtus agrestis	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	21	0	0	89	40	0	20	5	0	3	2	0	180
M. paradoxus	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	19	0	0	64	1	1	0	0	1	87
M. rossiaemeridionalis	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	70	2	0	108	4	0	1	0	1	188
Neodon juldaschi	2	3	86	2	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	94
Total	3	3	87	2	1	24	833	56	20	834	55	1	67	22	4	5	4	207	361	2	158	735	78	5	251	13	1	4	3	3	3842
Note: Modal morphotypes for each species according to [52] are shown in bold. Latin species names are given by [53], except absent from the fauna of Russia.

 

The following indices were used for the results interpretation: complexities and absolute asymmetry of the complexities (difference modulus) of closed fields number, labial, lingual sides of morphotypes, also minimal number of elementary transitions and Shannon entro- py [48].

Full nucleotide sequences of gene mtDNA Cytb (1140 bp) of the same gray vole species were chosen from the International GenBank [49] database. In total, 576 sequences were used (see Appendix).

The geographic coordinates and height above the sea level (also their ranks) of the places of collecting most representative samples of each species — total of 6 traits (Fig. 1) — were used as environmental characteristics. Besides, the average monthly climatic norms of these terrains — temperature, its amplitude and norm of precipitation (total of 36 traits) [50] — were analyzed.

 

Fig. 1. The locations of the most representative samples for each species. Species: 1 – A. fortis, 2 – A. maximowiczii, 3 – A. middendorffii, 4 – A. mongolicus, 5 – A. oeconomus, 6 – L. gregalis, 7 – M. agrestis, 8 – M. paradoxus, 9 – M. rossiaemeridionalis, 10 – N. juldaschi

 

A new algorithm for analyzing the congruence and combining the data based on the matrices of distances between the species, i.e. the DJ–method, was proposed in a number of our publications [23, 24, 51]. If a set of objects has two or more description systems, then it is possible to calculate the pair similarity/dissimilarity index matrix between objects on each of description system. The similarities/dissimilarity matrix can be transformed into the object–trait matrix. It is necessary in order to map the existing set of objects for each description system as the cloud of points in the Euclidean space. Thereafter, the cloud of points (objects) is analyzed for the congruence of their mutual location, e.g. assisted with the PLS–methods [54], generalized Procrustean analysis [55] or common correlation methods. The algorithm was realized in the Jacobi 4 package [56] of applied programs. In this work, the frequencies table of all pairs of morphotypes — symmetric and asymmetric — is normalized at the sum on the rows. The Cavalli-Sforza and Edwards matrix [46, 57] of Euclidean distances was calculated between all species. As this matrix was obtained on the base of differences in the occurrence of common pairs of morphotypes in different species, it reflects the variability of specific morphogenesis [45, 46]. On this basis, further on in the text, we will call it the morphogenetic matrix of distances.

The Cytb nucleotide sequences divided into codons were analyzed in the work. The transformation of nucleotide sequences from the GenBank into the Excel format and splitting into codons were carried out using MEGA7 and Excel packages.

It was believed earlier that, as synonymous nucleotide mutations do not change any amino acid they code, they have no consequences for the protein function or viability of an organism. However, as a result of a comparative analysis of sequences, the non-random synonymous codon frequencies distribution was revealed later in different organisms (individuals and species). This phenomenon was called the “Synonymous Codon Usage Bias” (SCUB) [58]. In our research, at the calculation of evolutionary distances, we used the model in which the substitution of not a nucleotide, but a codon [59, 60] is considered to be an elementary event. These substitutions can be both synonymous and nonsynonymous (replacements).

The matrix of p-distances that reflects the total number of codon substitutions for each pair of sequences was calculated between sequences. If amino acids are used instead of codons, then the p-distance matrix will reflect the number of only nonsynonymous substitutions of codons. The difference between them gives the p-distance matrix of only synonymous substitutions of codons. Any p-distance matrix is a matrix of squared Euclidean distances [51]. Therefore, square roots were extracted from the p-distances for all three molecular-genetic matrices. Then distance matrices were transformed into object–trait matrices with the method of principal coordinates [61]. This way, points in molecular-genetic Euclidean spaces were matched corresponding to all objects. PAST and Jacobi 4 packages were used to calculate distance matrices and principal components with the method of principal coordinates. Note that the method of principal coordinates is the direct method of calculating principal components in case when a matrix of distances between objects is Euclidean [61–64]. The centroid, which represents the species [65], was calculated for each set of points belonging to one species. Further on, matrices of pairwise Euclidean distances were calculated between centroids of all species.

All interspecific distance matrices, including the morphogenetic one, were processed with the method of principal coordinates. As a result of such processing, we obtained a limited set of modules (principal components) that describe their variability in morphogenetic and molecular–genetic spaces [66] for ten gray vole species under investigation.

In order to minimize the number of compared correlation coefficients, both totalities of geo–climatic factors (each separately), after the centering, normalization and calculating of Euclidean distance matrices, were processed with the method of principal coordinates. Part of calculations was carried out independently in different packages to control the correctness of the obtained results. Pearson linear correlation coefficients were calculated in packages Statistica 8, PAST and Jacobi 4, and their significance — in Statistica 8 and PAST packages (Student t-test for Pearson correlation coefficient, df = n-2). Next, the correlations of variability modules with the obtained geo-climatic gradients were analyzed.

RESULTS

Morphogenetic modules of variability. On the plane of the first two principal components of morphogenetic distances matrix (32.7 and 23.8 % of total variance), the species mutual location is determined by the differences of their specific morphogeneses (Fig. 2).

 

Fig. 2. The configuration of ten species of the gray voles Microtus s.l. on the plane of the first two PCs of the morphogenetic distances matrix

 

On the first principal morphogenetic component (mgPC1), N. juldaschi and A. oeconomus are drastically different from all the rest species due to the low average complexity of their morphotypes. Such interpretation of this variability module follows from its high correlations with the traits that characterize the complexity of m1 morphotypes (Table 2).

 

Table 2

Pearson correlation coefficients (×1000) between morphotypes m1 traits of ten species of the gray voles Microtus s.l. and two first PCs of the morphogenetic distances matrix

Trait	mgPC1	mgPC2
Total complexity of the right morphotype	971**	–188
Total complexity of both morphotypes	970**	–191
Total complexity of the left morphotype	969**	–194
Closed fields complexity of the left morphotype	903**	300
Closed fields complexity of the right morphotype	903**	297
Absolute asymmetry of morphotypes complexity	–067	693*
Minimal number of elementary transitions	–060	691*
Shannon entropy	251	682*
Absolute asymmetry of lingual side complexity	019	637*
Fraction of symmetric individuals	143	–632*
Note: *p < 0.05, **p < 0.001.

 

The investigated vole species are arranged along the second principal morphogenetic component (mgPC2) in the order of their increase m1 morphotypes asymmetry and the decrease of symmetric individuals fraction (see Table 2). As asymmetry is traditionally considered as an index of ontogenesis instability [67], it follows from Fig. 2 that species M. paradoxus, M rossiaemeridionalis and N. juldaschi are characterized by the most stable morphogenesis, and A. fortis and A. mongolicus — by the most unstable morphogenesis.

Molecular-genetic modules of variability. It was found out in the investigation process that the principal components of synonymous distances matrix demonstrate the highest congruence with one of the morphogenetic modules of variability out of three considered matrices (total, nonsynonymous and synonymous interspecific distances). Therefore, only they are discussed below. The species mutual location on the plane formed by these principal components (Fig. 3) well agrees to the results obtained with the methods of molecular phylogenetics [68], despite that our result was obtained on the model of replacing synonymous codons, but not single nucleotides.

 

Fig. 3. The configuration of ten species of the gray voles Microtus s.l. on the plane of the first two PCs of the synonymous distances matrix

 

At calculating the codon frequency correlations with the first principal components of synonymous distances matrix (genPC1), the highest correlations were found for two codons: GCA (r = –0.911, n = 10, p = 0.00025) and ACC (r = 0.902, n = 10, p = 0.00036). Thus, there is a significant increase of the ACC codon frequency and a decrease of the GCA codon frequency along the first principal component of genus Alexandromys voles, compared to those of genus Microtus.

Environmental factors. The first principal components of the matrices of geographic (62.2 % of total variance) and climatic (53.4 % of total variance) traits for the locations of capturing the most representative samples of each species highly correlate (r = 0.872, n = 10, p = 0.001) (Table 3). Totally, they build up the geo-climatic gradient “south-west–north-east”. There is a decrease in the average temperature of all months and winter-spring precipitation (from February to April) and an increase of summer precipitation (from June to September) along the gradient.

The congruence of the modules of morphogenetic, molecular–genetic and environmental variabilities was studied with the help of correlations between the first two components of the corresponding matrices (Table 3).

 

Table 3

Pearson correlation coefficients (×1000) between first two PCs of the morphogenetic (mgPC), molecular-genetic (genPC), geographic (geoPC), and climatic (climPC) variability

PCs	mgPC1	mgPC2	genPC1	genPC2	climPC1	climPC2	geoPC1	geoPC2
mgPC1	–	0	–058	–229	338	049	343	–242
mgPC2	0	–	847**	–260	695*	325	598	536
genPC1	–058	847**	–	0	788**	148	759*	472
genPC2	–229	–260	0	–	–047	–003	133	014
climPC1	338	695*	788**	–047	–	0	872**	368
climPC2	049	325	148	–003	0	–	–159	690*
geoPC1	343	598	759*	133	872**	–159	–	0
geoPC2	–242	536	472	014	368	690*	0	–
Note: *p < 0.05, **p < 0.001.

 

As is seen in Table 3, one of the strongest correlations (r = 0.847, n = 10, p = 0.002) was obtained between the first molecular-genetic and the second components of morphogenetic variability (Fig. 4).

 

Fig. 4. Correlation (r = 0,847, n = 10, p = 0,002) between PC1 of the synonymous distances matrix and PC2 of the morphogenetic distances matrix

 

The first component of molecular-genetic variability also correlates with the geo-climatic gradient, largely with the first climatic component (r = 0.788, n = 10, p = 0.0068) (Fig. 5). A weaker correlation (r = 0.695, n = 10, p = 0.026) is observed between the second component of the morphogenetic matrix and the first climatic component.

 

Fig. 5. Correlation (r = 0,788, n = 10, p = 0,0068) between PC1 of the climatic variability and PC1 of the synonymous distances matrix

 

The first principal component of the morphogenetic variability does not correlate with anything.

DISCUSSION

The present investigation is a continuation of a series of our articles published on the search of congruent variability of traits from different trait systems, also on the development of the adequate algorithms of such search [23, 24, 69]. An attempt of searching the congruence between the morphotypic and molecular-genetic variability of gray voles Microtus s.l. by analyzing the corresponding matrices of distances with the help of Mantel test [69] was taken earlier with the same material. The test showed the insignificant correlation between the matrices (r = 0.23, p = 0.13, Np = 106), from which a conclusion on the considerable independence of these trait systems of the studied species was made. Nevertheless, this issue has not been cleared out completely, as the 2B–PLS-analysis demonstrated the presence of congruent variability between the indices of morphotypes asymmetry and entropy, and the frequencies of gene Cytb nucleotides and dinucleotides along the geo-climatic gradient. To specify the investigation results, it was decided to use a bit other approach: namely, to split all the total variability of traits of different trait systems into autonomous units, i.e. modules that may contribute to revealing those of them that change jointly. At present, the modular approach to the analysis of variability is widely practiced worldwide. Different criteria (functional, ontogenetic, statistical, evolutionary, etc.) are used to reveal modules. In the present research, we used a statistical approach to pointing out variability modules [10, 11, 18], which are the principal components of two interspecific matrices of Euclidean distances.

The first matrix (morphogenetic) was built up based on the differences in the Microtus s.l. occurrence frequencies of the common pairs of ml; the second matrix (molecular-genetic) — based on the synonymous p-distances between the Cytb codon sequences of the same species. In this article, we consider two first modules for each trait system, as there was no congruent variability found beyond them. They are directions in the corresponding Euclidean spaces. The experience in this field shows that, in most cases, latent variables can be interpreted informatively [70–72].

Concerning the modular approach, ontogenesis is a totality of morphogeneses, which can be presented as discrete blocks or modules that can be changed relatively independently, other events of ontogenesis being unaffected [73, 74]. As, in our research, the interspecific morphogenetic distances matrix was obtained on the base of the difference in the occurrence frequencies of common pairs of m1 morphotypes, both symmetric and asymmetric, it reflects the variability of specific morphogeneses. Calculating the principal components of this matrix allows its decomposition into independent constituents, i.e. modules which are a result of different morphogenetic events. In this case, these modules were interpreted by us as the average morphogenetic complexity, which reflects the complexity of the molar prismatic design, and the asymmetry of the phenotypic complexity which reflects the instability of the morphogenesis of this organ.

The species mutual location in Fig. 2 is determined by the differences of their specific morphogeneses. Two species N. juldaschi and A. oeconomus are characterized by the low average complexity of their morphotypes (separated on the first principal component), which distinguishes them from the rest species that stretch along the second component forming pairs. Thus, the pair of M. paradoxus, M. rossiaemeridionalis has the most stable morphogenesis, and the pair of A. fortis, A. mongolicus has the least stable morphogenesis, as the latter have the highest asymmetry indices of morphotypic complexity and the fraction of asymmetric individuals. Interesting is that no significant correlations with a direction of molecular-genetic or geo-climatic variability were found for the first principal component, which is the main direction of morphogenetic variability and reflects the variability of the average m1 morphotypic complexity in the investigated species.

Next, we obtained the correlation of the second component of the morphogenetic matrix with the “south-west–north-east” geo-climatic gradient (especially with the first climatic principal component (Table 3), which confirms the mophotypes asymmetry at the increase of stress in the “organism–environment” interaction (Table 2). As the first climatic principal component is formed by exclusively the contribution of temperatures, we believe temperature is the main interspecific variability factor of morphogeneses [75–77].

However, the second morphogenetic component demonstrated an even higher correlation with the first component of interspecific molecular-genetic distances matrix calculated on the model of the synonymous substitutions of Cytb codons (Table 3). The fact of SCUB established for a big number of species of living organisms serves as base for the calculation. The ratio of synonymous codons in mRNA is an important fact that influences the translation speed and protein synthesis kinetics [78–81]. Besides, the differences of synonymous codons, actually, serve as additional code of controlling protein folding in vivo and, hence, they are one more source of protein conformation variability [82].

The reasons for SCUB remain the subject of discussions. In particular, the hypotheses of random mutations and directional selection were formulated [83]. The result obtained within our research, namely, the correspondence of the first principal component of interspecific distances matrix to the geo-climatic gradient “south-west–north-east” (Table 3), is the evidence in favor of the selective formation mechanism of the observed SCUB. A significant increase of the ACC codon frequency and a decrease of the GCA codon frequency occur along the geo-climatic gradient in the voles of genus Alexandromys, compared to those of genus Microtus. The obtained SCUB may be conjugated with the fast adaptive changes of the secondary protein structure conditioned by the impact of environmental factors and that, within the present research, is just an assumption.

Acknowledgments

The research was supported by the Fundamental Scientific Research Program (FSR) of State Academies of Sciences for 2013–2020, project № VI.51.1.9. (AAAA-A16-116121410119-4) and RFBR grants №№ 17-04-00269-a, 19-07-00658-a. The material from the collection of Siberian Zoological Museum, Novosibirsk, of the Institute of Systematics and Ecology of Animals SB RAS (state task № AAAA-A17-117101070002-6), supported by the Program of Bioresource Collections Development of the Ministry of Science of Russia, was used.

APPENDIX

GenBank accession numbers of sequences:

218952. agrestis: AF119271, AF159402, AY167149–AY167213, FJ619746–FJ619786, GU563195–GU563299, JX284248–JX284283, KF218851–KF218952.
218953. fortis: AF163894, EU126809, EU870632, HQ123607, HQ123608, JF906125, KJ081871–KJ081954.
218954. gregalis: AF163895, AY513803, KF839591, KF839592.
218955. juldaschi: AY513808, EF599112, EF599113.
218956. maximowiczii: FJ986303, FJ986311, FJ986312, FJ986319, KJ857275–KJ857280, KJ857287–KJ857291.
218957. middendorffii: AF163898, FJ986314, FJ986315, HM119493.
218958. mongolicus: FJ986304, FJ986305, FJ986309, FJ986310.
218959. oeconomus: AB372193–AB372207, AF163902, AY219981–AY220045, DQ452134–DQ452142.
218960. paradoxus: KC953622, KC953624.
218961. rossiaemeridionalis: AY513819–AY513823, U54472–U54478, U54493, U54495.

About the authors

Vera Yu. Kovaleva

Federal Institute of Systematics and Ecology of Animals of the Siberian Branch of the Russian Academy of Sciences

Author for correspondence.
Email: vkova@ngs.ru
ORCID iD: 0000-0003-1685-4820

DSc, Sen Sci Res, Laboratory of Animal Ecology

Russian Federation, 11, Frunze street, Novosibirsk, 630091

Alexandr A. Pozdnyakov

Federal Institute of Systematics and Ecology of Animals of the Siberian Branch of the Russian Academy of Sciences

Email: al_ap@mail.ru

PhD, Sen Sci Res, Laboratory of Animal Ecology

Russian Federation, 11, Frunze street, Novosibirsk, 630091

Yuri N. Litvinov

Federal Institute of Systematics and Ecology of Animals of the Siberian Branch of the Russian Academy of Sciences

Email: Lyun13@yandex.ru

DSc, Deputy Director, Laboratory of Animal Ecology

Russian Federation, 11, Frunze street, Novosibirsk, 630091

Vadim M. Efimov

Federal Institute of Systematics and Ecology of Animals of the Siberian Branch of the Russian Academy of Sciences; Novosibirsk State University; Tomsk State University

Email: vmefimov@ngs.ru

DSc, Prof. of Cytology and Genetics Department; Prof. of Vertebrate Zoology and Ecology Department

Russian Federation, 11, Frunze street, Novosibirsk, 630091; 1, Pirogova street, Novosibirsk, 630090; 34a, Lenina prospect, Tomsk, 634050

References

Supplementary files