Abstract

The Major Histocompatibility Complex (MHC), as a family of highly polymorphic genes associated with immunity in the genome of the vertebrate, has become an important indicator for assessing the evolutionary potential of wildlife. In order to better protect Zootoca vivipara in the Greater Khingan Range and Lesser Khingan Range, to understand the genetic structure of Z. vivipara, and to explore the mechanism and phylogenetic relationship of the gene polymorphisms, the MHC molecular marker method was used to analyze Z. vivipara population. Forty-seven alleles were obtained from four populations. The four populations were highly polymorphic, rich in genetic information, and had significant genetic diversity. There were certain inbreeding phenomena. There was a high degree of genetic differentiation among populations, which was caused by genetic drift and natural selection. The sequence undergoes genetic duplication and recombination. The existence of trans-species polymorphism was found in the constructed phylogenetic tree. The present study provides a theoretical basis for species protection of Z. vivipara.

The Major Histocompatibility Complex (MHC) is a highly polymorphic receptor multigene family in the genome of the vertebrate, which is closely related to the immune system. Due to the need to identify multiple pathogens and adapt to the environment, MHC genes are polygenic and polymorphic, and are related to the survival and reproduction of the population, as well as immunity against disease. This is the unique advantage of MHC in population genetics and conservation genetics. Therefore, MHC has been considered to protect the latest genetic systems in genetics research. MHC molecular markers can overcome such limitations as maternal inheritance of mitochondrial genes and time consumption of microsatellites, so that MHC genes have become a hot genetic marker in population genetic structure and gene variation analysis, and have a broad application prospect in species conservation genetics [1,2]. Since the characteristics of MHC class I genes of Nerodia sipedon and Ameiva ameiva were first reported in 1992, research on reptile MHC genes has been increasing year by year [3]. Researchers often use MHC to study endangered animals and new species. For example, in the last decade, there have been many studies on Sphenodon punctatus. The researchers successfully obtained sequence information and analyzed their genetic structure and population dynamics [4–7]. Later, when analyzing the MHC of some endangered or newly discovered species [8,9], people further studied from the aspects of species relationship, survival and reproduction, and disease resistance [10–17].

Zootoca vivipara is a typical representative of the viviparous animals of the lizard family [18]. It is also a representative of the Palaearctic realm lizards [19]. Z. vivipara originated in the Mediterranean region. After the Glacial period, the natural selection made life more suitable in high latitudes, cold areas. Z. vivipara in China are all oviparous reproductive mode, distributed in Heilongjiang, Inner Mongolia and Xinjiang [20]. Its distribution range is relatively narrow, the population is small, and the distribution area is special. It is the southern limit of the distribution of the species, and the geographical interval of the distribution is large, which has certain research and protection significance. In recent years, with the development of molecular biology, the research on Z. vivipara has been deepened. In the origin and evolution of species, the discussion of the double-breeding model of oviparous-viviparous is a hot topic. Different genetic analysis has caused several differences [21–26]. At the same time, the use of molecular marker technology to explore sequence genetic polymorphism is quite concerned by researchers [27–33]. However, MHC molecular studies of Z. vivipara have not been reported.

Due to the destruction of the ecological environment caused by human disturbance, the habitat of the global Z. vivipara has been shrinking, and the number of Z. vivipara has gradually decreased [34]. At present, several populations in China are already in the Near Threatened status [35]. Therefore, using MHC genes to study the four geographic populations of Z. vivipara in China can better describe the Z. vivipara gene polymorphisms and population genetic dynamics, and further explain the effects of natural environment and human factors on phylogenetic relationships and existing distribution characteristics. By exploring its genetic polymorphism formation mechanism, the survival status and fitness of the Z. vivipara population can also be speculated immunologically. Furthermore, it reveals the evolutionary history of Z. vivipara population and the formation mechanism of systematic geographical pattern in the study area. This also provides theoretical support for the scientific protection of Z. vivipara species in the next step.

Materials and methods

Sample collection and molecular techniques

Thirty-four adult Z. vivipara males and females were randomly captured by hand-hunting method from four locations of the east and west of the Greater Khingan Range and Lesser Khingan Range: Ya Keshi county, Inner Mongolia autonomous region; Hu Ma county, Nen Jiang county and Sun Wu county, Heilongjiang province. The lizard’s tail tissue was collected and stored in liquid nitrogen, then the lizard was released back to its original site. Genomic DNA was extracted using the SanPrEP column animal genomic DNA extraction kit of Sangon Biotech (Shanghai) Co., Ltd., and total RNA was extracted by Takara’s TRIzol reagent. Three pairs of primers IMHCF/IMHCR, 2MHCF/2MHCR, and MHCIF/MHCIR were utilized to clone the mRNA fragment of MHC I, the DNA fragment of MHC I α2, and the DNA fragment of MHC I α3 [36]. SSCP analysis and transformation and cloning sequencing were performed separately.

Allele determination and sequence alignment

Sequencing successful sequence peak maps were read with Chromas 2.6.5 software, and Blast homologous alignments (https://blast.ncbi.nlm.nih.gov/Blast.cgi) were performed on the NCBI webpage to determine positive clones, whether the sequence is a gene. The sequence alignment used Clustalx 2.1 with default parameters. The sequence alignment map was run with DNAMAN 9.0 and the predicted amino acid sequence was translated using MEGA 7.0 [37] and placed in the Expasy reading frame (http://web.expasy.org/translate/). The information of the alleles was acquired by sequence alignment and calculation of the polyacrylamide gel image, SSCP method. The two methods mutually verify that the sequence occurring simultaneously in two or more individuals is judged as an allele [38].

Calculation of genetic polymorphism

Nucleotide polymorphism (π) and variable sites (S) were calculated by using Dnasp [39]. Popgene 1.32 software was used to calculate the number of alleles (Na), the number of effective alleles (Ne), the expected heterozygosity (He), and the observed heterozygosity (Ho) and the Shannon information index (I). Polymorphic information content (PIC) [40], gene flow (Nm), Inbreeding coefficient (Fis), and population fixation index (Fst) were calculated by using Cervus 3.0 software. Then MEGA was used for Tajima’s D neutrality test. Amino acid sequence alignments were performed by using online COBALT (www.ncbi.nlm.nih.gov/tools/cobalt/). According to the recent MHC I gene research, the gene sequence of Z. vivipara relatives was downloaded from GenBank and rooted in human HLA-A and HLA-G. MHC gene homology relationship of four populations of Z. vivipara was detected by MEGA software. Then, based on the Poisson-based model (Neighbor-joining, NJ), the phylogenetic tree and branches are constructed for its nucleotide and amino acid sequences. The node confidence (BP) was acquired by 1000 repeated sampling using Bootstrap. In addition, the amino acid sequence of Z. vivipara relatives was downloaded from GenBank as an outer group, with human HLA-A, HLA-B, HLA-C, HLA-G, and HLA-I.

Selection and recombination detection

MEGA software calculated the nucleotide distance (Kimura 2-parameter model, K2P), amino acid distance (Poisson-corrected model), the ratio of non-synonymous mutation rate to synonymous mutation rate (dN/dS), self-expanding repetition value 1000 [41]. The selection site was tested using the Z-test of MEGA to calculate the sequence selection site and was detected using the FEL, MEME, FUBAR model (http://www.datamonkey.org/) [42]. Restructuring tests were performed using the GRAD model on the Datamonkey website (http://www.datamonkey.org/).

Results and discussion

Genetic diversity

The present study amplified MHC I genes of four populations of Z. vivipara in China for the first time, and cloned the complete exons 3 and 4 regions from the DNA and RNA molecule level for subsequent analysis. A total of 570 sequences were successfully cloned from 34 individuals of four populations of Z. vivipara. After comparison with Blast sequences, they were confirmed to be classical MHC class I functional genes [38]. The presence of non-canonical MHC class I genes was not found in this study, probably due to the small sample size of four populations and the small number of ten positive clones selected per individual, resulting in the failure to acquire all of the alleles. Among them, the acquired MHC I α2 has a nucleotide length of 201–223 bp, α3 has a nucleotide length of 181–184 bp, and α2 and α3 have an amino acid length of 111–178 aa. The number of alleles of four populations of Z. vivipara studied in this experiment differed, including 3–11 alleles. So, it was found that the same allele appeared differently in different individuals of different populations (Figure 1), and this difference also reflected the polymorphism of MHC class I genes [1,2]. This study used the number and calculation of alleles to speculate on the number of loci per Z. vivipara population. The results showed that there were at least three MHC I loci in Nen Jiang population (N). There were at least four MHC I loci in Sun Wu population (SW). There were at least six MHC I loci in Hu Ma population (H) and there were also at least six MHC I loci in Ya Keshi population (Y). Among them, the highest π of Y was 0.659, N had a lower π of 0.305, N had the most S in MHC class I genes (580), and H had the least S (306) (Table 1). It can be found that the nucleotide polymorphisms of MHC class I alleles of four populations of Z. vivipara were higher. Due to the current strict computational model, some loci may be ignored. However, the results of the analysis still show that there was a high degree of polymorphism at the locus level of four populations of Z. vivipara. It was speculated that MHC class I gene in Z. vivipara had the phenomenon of gene duplication and gene deletion in the evolution process, which is consistent with the ‘Birth-and-death model of evolution’ [43].

Sample map of Zootoca vivipara

Figure 1
Sample map of Zootoca vivipara

Pie charts represent the proportion of individuals with different allele numbers in various groups. The population codes are as follows: Y indicates Ya Keshi population; H indicates Hu Ma population; N indicates Nen Jiang population; SW indicates Sun Wu population.

Figure 1
Sample map of Zootoca vivipara

Pie charts represent the proportion of individuals with different allele numbers in various groups. The population codes are as follows: Y indicates Ya Keshi population; H indicates Hu Ma population; N indicates Nen Jiang population; SW indicates Sun Wu population.

Table 1
Nucleotide polymorphism analysis of MHC class I alleles of four populations of Z. vivipara
PopulationsNumber of lociNucleotide polymorphism, πNumber of variant sites, S
Nen Jiang (N) 0.305 580 
Sun Wu (SW) 0.478 372 
Hu Ma (H) 0.385 306 
Ya Keshi (Y) 0.659 529 
PopulationsNumber of lociNucleotide polymorphism, πNumber of variant sites, S
Nen Jiang (N) 0.305 580 
Sun Wu (SW) 0.478 372 
Hu Ma (H) 0.385 306 
Ya Keshi (Y) 0.659 529 

Alleles were named as Zovinj, Zovisw, Zovihm, and Zoviyks. The 11 mRNA sequences acquired were translated into amino acid sequences by MEGA software. MHC I sequences of 11 other animals were used as the outer group for comparison. When amino acid sequence alignments were performed, it was found that the putative antigen binding sites were all on the α2 domain and the β-2 microglobulin binding site was on the α3 domain. It can be observed in the difference in sequence that the polymorphism of the α2 domain was higher than that of the α3 domain, which may be related to the recognition of multiple antigens and the maintenance of protein structure. Using DANMAN to describe the homology relationship of sequences, it was found that the homology of different alleles in various groups was high. The homology of populations divided four populations into two categories: N and H had high homology, first clustered into one, and then clustered with SW. Finally, they were clustered together with Y, and the geographical distance from the populations were inconsistent. From the phylogenetic relationship, the four populations of Z. vivipara were mainly clustered in the NJ tree, and the self-supporting rate was higher, there was no significant difference. In the NJ tree at the DNA level, H was independent of one branch in the α2 domain, and SW was independent of one branch in the α3 domain. Respectively in H and SW appeared to branch, which may be due to the presence of positive and negative sequencing during cloning sequencing. At the same time, these sequences were under the same root, the difference was not significant. Then the phylogenetic analysis of the full amino acid sequence of α2 and α3 domains. It was found that four populations of Z. vivipara were clustered into two (Figure 2). This difference may be caused by the different evolutionary relationships between different domains caused by gene recombination [14]. It also indicated that MHC class I gene of Z. vivipara had high polymorphism as a whole. It is consistent with the previous analysis of the genetic diversity of Z. vivipara population in China [31,32]. In the above analysis of various levels of polymorphism, it was found that all four populations were highly polymorphic, which may be related to the ability to resist disease.

Amino acid phylogenetic tree of MHC class I gene of Z. vivipara based on NJ method

Figure 2
Amino acid phylogenetic tree of MHC class I gene of Z. vivipara based on NJ method

Bootstrap values from 1000 iterations were indicated above the branches.

Figure 2
Amino acid phylogenetic tree of MHC class I gene of Z. vivipara based on NJ method

Bootstrap values from 1000 iterations were indicated above the branches.

Genetic structure

In all populations, Na was generally greater than Ne. This indicated that alleles were unevenly distributed among the populations. The Ho of populations were generally less than He. This showed that there was a certain genetic variation [44]. PIC was generally greater than 0.5, belonging to the high polymorphic locus [40]. Using the Grant et al. method to estimate the demographic, it was found that the four populations of Z. vivipara were the result of long-term differentiation [45]. Compared with other reptilian animals, the polymorphism coefficient was found to be lower than that of the widely distributed species, but higher than several endangered species. It indicated that four populations of Z. vivipara may be more interfered, which is consistent with the near threatened status. The gene richness is directly proportional to I. Since Z. vivipara has I between 0.5 and 1, the gene was abundant [46]. Most Ho was less than He and Fis was greater than 0. Inbreeding was common within the population, and distant breeding occurred occasionally. This result was consistent with the mating method of Z. vivipara and the reduction in the number of individuals in the population. Tajima’s Test of Neutrality was greater than 0 indicating that the population was affected by balanced selection, and population contraction may have occurred, which was consistent with the increased survival risk of Z. vivipara populations (Table 2). Fst was generally was greater than 0.25, and there was a high degree of genetic differentiation between the populations [47]. Gene flow is an important factor affecting the genetic structure of the population. Meanwhile, the Nm between Z. vivipara populations was less than 1, genetic drift became the dominant factor in the differentiation of the population genetic structure [48], indicating less genetic exchange between the populations (Table 3). It was speculated that the effect of genetic drift causes the four populations of Z. vivipara to differentiate, and there may be some decline between the populations [49]. This may be related to habitat destruction and reduction.

Table 2
The genetic diversity index of four populations of Z. vivipara in MHC class I
PopulationsSitesNumber of alleles, NaNumber of effective alleles, NeObserved heterozygosity, HoExpected heterozygosity, HePolymorphic information content, PICShannon information index, IInbreeding coefficient, FisTajima’s D
Nen Jiang (N) Full-length 9.0000 7.1053 1.0000 0.8947 0.5550 0.7500 −0.1429 1.722 
 (α2) Exon3 3.0000 2.4000 0.0000 0.6000 0.4070 0.7000 0.2174 2.650 
 (α3) Exon4 3.0000 2.4000 0.0000 0.6000 0.5660 0.6250 0.1128 3.759 
Sun Wu (SW) Full-length 2.0000 1.8000 0.0000 0.3333 0.6670 0.6389 0.3333 4.295 
 (α2) Exon3 3.0000 2.5714 1.0000 0.4286 0.4287 0.7000 −0.2903 4.269 
 (α3) Exon4 2.0000 1.6667 0.0000 0.3333 0.3333 0.7130 0.3378 3.530 
Hu Ma (H) Full-length 2.0000 1.6667 0.0000 0.6000 0.8330 0.7045 0.2800 0.902 
 (α2) Exon3 6.0000 4.8462 1.0000 0.1538 0.7500 0.6974 −0.3097 3.383 
 (α3) Exon4 9.0000 7.1053 1.0000 0.8947 0.5858 0.7115 −0.0800 2.360 
Ya Keshi (Y) Full-length 2.0000 1.6667 0.0000 0.6000 0.3860 0.7045 0.0150 2.323 
 (α2) Exon3 6.0000 4.6364 0.0000 0.3636 0.6050 0.7130 0.0368 2.922 
 (α3) Exon4 3.0000 2.4000 0.0000 0.6000 0.7440 0.7000 0.0474 3.938 
PopulationsSitesNumber of alleles, NaNumber of effective alleles, NeObserved heterozygosity, HoExpected heterozygosity, HePolymorphic information content, PICShannon information index, IInbreeding coefficient, FisTajima’s D
Nen Jiang (N) Full-length 9.0000 7.1053 1.0000 0.8947 0.5550 0.7500 −0.1429 1.722 
 (α2) Exon3 3.0000 2.4000 0.0000 0.6000 0.4070 0.7000 0.2174 2.650 
 (α3) Exon4 3.0000 2.4000 0.0000 0.6000 0.5660 0.6250 0.1128 3.759 
Sun Wu (SW) Full-length 2.0000 1.8000 0.0000 0.3333 0.6670 0.6389 0.3333 4.295 
 (α2) Exon3 3.0000 2.5714 1.0000 0.4286 0.4287 0.7000 −0.2903 4.269 
 (α3) Exon4 2.0000 1.6667 0.0000 0.3333 0.3333 0.7130 0.3378 3.530 
Hu Ma (H) Full-length 2.0000 1.6667 0.0000 0.6000 0.8330 0.7045 0.2800 0.902 
 (α2) Exon3 6.0000 4.8462 1.0000 0.1538 0.7500 0.6974 −0.3097 3.383 
 (α3) Exon4 9.0000 7.1053 1.0000 0.8947 0.5858 0.7115 −0.0800 2.360 
Ya Keshi (Y) Full-length 2.0000 1.6667 0.0000 0.6000 0.3860 0.7045 0.0150 2.323 
 (α2) Exon3 6.0000 4.6364 0.0000 0.3636 0.6050 0.7130 0.0368 2.922 
 (α3) Exon4 3.0000 2.4000 0.0000 0.6000 0.7440 0.7000 0.0474 3.938 
Table 3
Genetic divergence (below diagonal) and gene flow (above diagonal) among four populations of Z. vivipara in the MHC class I
PopulationsNen Jiang (N)Sun Wu (SW)Hu Ma (H)Ya Keshi (Y)
Nen Jiang (N) — 0.8202 0.9912 0.4119 
Sun Wu (SW) 0.2835 — 0.7331 0.3807 
Hu Ma (H) 0.2302 0.2575 — 0.5878 
Ya Keshi (Y) 0.3981 0.4017 0.3202 — 
PopulationsNen Jiang (N)Sun Wu (SW)Hu Ma (H)Ya Keshi (Y)
Nen Jiang (N) — 0.8202 0.9912 0.4119 
Sun Wu (SW) 0.2835 — 0.7331 0.3807 
Hu Ma (H) 0.2302 0.2575 — 0.5878 
Ya Keshi (Y) 0.3981 0.4017 0.3202 — 

The genetic divergence (below diagonal) and gene flow (above diagonal).

Positive selection and recombination detection

In the present study, MHC class I gene was highly polymorphic, and its causes were analyzed and its maintenance mechanism was speculated. First, the genetic distance analysis showed that the ratio of the average nucleotide distance to the average amino acid distance of the four populations was generally less than 1, and the difference between each gene region was extremely significant (Table 4). Second, the ratio of the non-synonymous mutation rate to the synonymous mutation rate in different gene domains of different populations ω (dN/dS) was generally greater than 1, suggesting that MHC I genes of the four populations may be selected (Table 5). Then four methods were used to perform positive selection calculations on four different gene domains. Some amino acid sites that were positively selected existed near the putative antigen binding site [45]. Finally, these analyses indicate that Z. vivipara’s MHC class I genes were likely to have undergone positive selection. The pathogen-mediated selection of pathogens and parasites is one of the main reasons for maintaining high polymorphism, indicating that the positive selection was highly polymorphic with MHC gene of Z. vivipara and the regulation of pathogens by Z. vivipara (Table 6). Since the degree of variation of the antigen binding site is positively related to the type of antigen that can be recognized, which plays an important role in antigen recognition and presentation in the body’s immune response. This is similar to other vertebrate immune mechanisms [6,41,50]. Z. vivipara likes to live in a hidden environment, where it is mostly dark and humid, breeding many parasites and pathogens. This situation may be the reason why Z. vivipara was receiving positive selection. However, some positive selection mutation sites are not on the putative antigen-binding region, which may also be related to some important sites different from humans. The habitat of Z. vivipara is subject to certain human disturbances, and this effect may also alter the outcome of parasite-mediated selection [17].

Table 4
MHC I genetic divergence of nucleotide and amino-acid within four populations of Z. vivipara
PopulationsAve. nucleotide divergence/Ave. amino-acid divergence
Full-length(α2) Exon3(α3) Exon4
Nen Jiang (N) 0.314 0.606 0.470 
Sun Wu (SW) 0.210 0.417 0.679 
Hu Ma (H) 0.350 0.873 0.344 
Ya Keshi (Y) 0.406 0.458 0.380 
PopulationsAve. nucleotide divergence/Ave. amino-acid divergence
Full-length(α2) Exon3(α3) Exon4
Nen Jiang (N) 0.314 0.606 0.470 
Sun Wu (SW) 0.210 0.417 0.679 
Hu Ma (H) 0.350 0.873 0.344 
Ya Keshi (Y) 0.406 0.458 0.380 
Table 5
Ratio of non-synonymous mutation rate to synonymous mutation rate
PopulationsSitesNon-synonymous, dNSynonymous, dS(ω) dN/dSP-value
Nen Jiang (N) Full-length 0.291 0.303 0.960 0.020 
 (α2) Exon3 0.661 0.339 1.949 0.016 
 (α3) Exon4 0.667 0.333 2.000 0.060 
Sun Wu (SW) Full-length 0.285 0.283 1.007 0.050 
 (α2) Exon3 0.726 0.274 2.650 0.019 
 (α3) Exon4 0.769 0.231 3.333 0.048 
Hu Ma (H) Full-length 0.152 0.149 1.020 0.069 
 (α2) Exon3 0.706 0.294 2.406 0.011 
 (α3) Exon4 0.707 0.293 2.412 0.021 
Ya Keshi (Y) Full-length 0.349 0.259 1.347 0.028 
 (α2) Exon3 0.887 0.113 7.875 0.013 
 (α3) Exon4 0.686 0.314 2.182 0.020 
PopulationsSitesNon-synonymous, dNSynonymous, dS(ω) dN/dSP-value
Nen Jiang (N) Full-length 0.291 0.303 0.960 0.020 
 (α2) Exon3 0.661 0.339 1.949 0.016 
 (α3) Exon4 0.667 0.333 2.000 0.060 
Sun Wu (SW) Full-length 0.285 0.283 1.007 0.050 
 (α2) Exon3 0.726 0.274 2.650 0.019 
 (α3) Exon4 0.769 0.231 3.333 0.048 
Hu Ma (H) Full-length 0.152 0.149 1.020 0.069 
 (α2) Exon3 0.706 0.294 2.406 0.011 
 (α3) Exon4 0.707 0.293 2.412 0.021 
Ya Keshi (Y) Full-length 0.349 0.259 1.347 0.028 
 (α2) Exon3 0.887 0.113 7.875 0.013 
 (α3) Exon4 0.686 0.314 2.182 0.020 

The difference is significant when P<0.05.

Table 6
Positive selection site analysis of Z. vivipara
Gene sitesMethodSites predicted to be under positive selection
Full-length  10 11 14 15 17 23 24 28 29 30 34 35 40 41 43 44 56 60 64 74 85 104       Sum 
 MEGA                        11 
 FEL                                
 MEME              21 
 FUBAR                                  
(α2) Exon3  11 12 14 16 17 18 20 25 26 33 37 39 40 41 42 43 44 46 49 50 51 54 60 61 63 Sum 
 MEGA      29 
 FEL                                  
 MEME                              
 FUBAR                                  
(α3) Exon4  10 15 17 21 22 23 25 27 28 30 32 34 35 39 40 41 44 53 61 72 100              Sum 
 MEGA                     14 
 FEL                                
 MEME                               
 FUBAR                                  
Gene sitesMethodSites predicted to be under positive selection
Full-length  10 11 14 15 17 23 24 28 29 30 34 35 40 41 43 44 56 60 64 74 85 104       Sum 
 MEGA                        11 
 FEL                                
 MEME              21 
 FUBAR                                  
(α2) Exon3  11 12 14 16 17 18 20 25 26 33 37 39 40 41 42 43 44 46 49 50 51 54 60 61 63 Sum 
 MEGA      29 
 FEL                                  
 MEME                              
 FUBAR                                  
(α3) Exon4  10 15 17 21 22 23 25 27 28 30 32 34 35 39 40 41 44 53 61 72 100              Sum 
 MEGA                     14 
 FEL                                
 MEME                               
 FUBAR                                  

Recombinant analysis of MHC I gene sequences of four populations of Z. vivipara were performed using online software, and traces of three recombination effects (81, 149, and 161 points) were detected in the MHC sequence. And the difference was extremely significant (P<0.01). Recombination of MHC genes behaves differently in different species of vertebrates. In Osteichthyes and Amphibians, exon shulffing is caused by different evolutionary relationships of different domains [51]. In the study of the birds MHC gene, it is found that gene recombination and transformation occurred intralocus [52,53]. While in mammals only mini-casettes are recombined [54,55]. This difference is due to the different stages of the evolution of MHC I gene. Intralocus recombination, especially allelic recombination, is often found in reptiles MHC, which can produce high polymorphism of MHC and maintain it. According to the reported sequence characteristics of MHC gene in geckos, it is found that the polymorphism of MHC gene mainly comes from gene recombination [56]. There was also a gene recombination event in the high polymorphism MHC sequence of Ctenophorus decresii [14]. In this study, three genetic recombinations were found in the MHC sequence of Z. vivipara using the GRAD model, which was located between exon 3 and exons 3 and 4. It indicates that the sequence had been broken and the gene exchange during the evolution process. Therefore, genetic recombination also plays a role in maintaining the polymorphism of the MHC class I gene in Z. vivipara. However, due to the current lack of certain rules for the genetic recombination of non-mammalian MHC class I genes, different environmental factors and gene structures of different species may lead to different genetic recombination results. Moreover, no recombination phenomenon was found in the analysis of three mitochondrial gene regions and three nuclear gene loci in the Eurasian region [33], and it is not possible to determine whether this recombination event was prevalent in all populations of Z. vivipara. Finally, the small population size, number of individual samples, and number of alleles in the present study may also limit the discovery of recombination phenomena.

Cross-species polymorphism analysis

The phylogenetic NJ tree of MHC I α2 and α3 full amino acid sequences was constructed by using the amino acid sequence used in the sequence alignment as exogroups (Figure 2). The phylogenetic NJ tree of MHC I α2 and α3 was constructed by using the 22 related sequences downloaded from NCBI acquired as an outer group. The study found that the four populations of Z. vivipara were clustered from the phylogenetic tree of MHC I mRNA. So, there was no obvious differentiation of four populations of Z. vivipara, indicating that there were shared alleles between them. This is consistent with the results of previous genetic polymorphism. In addition, in the three phylogenetic trees, MHC I gene sequence of the Eremias przewalskii always clustered with Z. vivipara MHC I gene sequence. The examination of the whole experiment and analysis process was confirmed to be correct, indicating that there was a cross-species polymorphism in Z. vivipara MHC I gene. The species belonged to Zootoca spp. and Eremias spp. The phylogenetic relationship in the phylogenetic tree was greater than that of MHC I gene system, suggesting that MHC I gene already existed at least before the differentiation of Z. vivipara and E. przewalskii. The species differentiation of Z. vivipara has been identified as 4.5 million years ago (95% confidence interval: 2.6–6.1 million years ago) [22], and the species differentiation of Eremias spp. has been identified as 6.3 million years ago (95% confidence interval: 5.3–8.5 million years ago) [57]. It was speculated that Z. vivipara MHC I allele lineage was formed at least 6 million years ago. Therefore, it may also be related to the adaptation of the lizards to the environment of high latitude and high cold regions in the evolution process [58]. This situation is consistent with the theory that MHC is the original gene that appeared before the formation of existing species [59].

In the present study, experiments were performed on MHC genes of 34 Z. vivipara individuals in four populations. The results showed that four populations of Z. vivipara in China had high genetic polymorphisms. The gene replication and gene recombination experienced by MHC genes during evolution were maintained by positive selection. At the same time, Z. vivipara had inbreeding in the population, and gene exchange and higher degree of genetic differentiation had occurred between the populations by the combination of genetic drift and natural selection. The MHC gene of Z. vivipara even had trans-species polymorphism. Although the genetic polymorphism of the Z. vivipara MHC gene and its formation mechanism had been analyzed in the present study, the immunological mechanism of the MHC gene family has not been fully revealed. Since the conservation of Z. vivipara is continuing, the molecular immunology of Z. vivipara also needs to be further analyzed.

Competing Interests

The authors declare that there are no competing interests associated with the manuscript.

Funding

This work was supported by the Science Foundation of Heilongjiang Province [grant number C2016035].

Author Contribution

Wanli Liu: performed the experiments, analyzed the data and wrote the paper. Yufen Liu, Peng Liu and Wenge Zhao: conceived and designed the present study, and revised the paper.

Abbreviations

     
  • MHC

    major histocompatibility complex

  •  
  • PIC

    polymorphic information content

  •  
  • SSCP

    single strand conformation polymorphism

References

References
1.
Zinkernagel
R.M.
and
Doherty
P.C.
(
1974
)
Immunological surveillance against altered self components by sensitized T lymphocytesin lymphocytic choriomeningitis
.
Nature
251
,
547
548
[PubMed]
2.
Klein
J.
(
1986
)
The Natural History of the Major Histocompatibility Complex
.
John Wiley and Sons
,
New York
3.
Grossberger
D.
and
Parham
P.
(
1992
)
Reptilian class I major histocompatibility complex genes reveal conserved elements in class I structure
.
Immunogenetics
36
,
166
174
[PubMed]
4.
Miller
H.C.
,
Belov
K.
and
Daugherty
C.H.
(
2006
)
MHC class I genes in the Tuatara (Sphenodon spp.): evolution of the MHC in an ancient reptilian order
.
Mol. Biol. Evol.
23
,
949
956
5.
Miller
H.C.
,
Cookson
M.A.
and
Daugherty
C.H.
(
2007
)
Two patterns of variation among MHC class I loci in Tuatara (Sphenodon punctatus)
.
J. Hered.
98
,
666
677
[PubMed]
6.
Miller
H.C.
,
Miller
K.A.
and
Daugherty
C.H.
(
2008
)
Reduced MHC variation in a threatened tuatara species
.
Anim. Conserv.
11
,
206
214
7.
Miller
H.C.
,
Moore
J.A.
,
Nelson
N.J.
and
Daugherty
C.H.
(
2009
)
Influence of major histocompatibility complex genotype on mating success in a free-ranging reptile population
.
Proc. R. Soc. B Biol. Sci.
276
,
1695
1704
8.
Stiebens
V.A.
,
Merino
S.E.
,
Chain
F.J.J.
and
Eizaguirre
C.
(
2013
)
Evolution of MHC class I genes in the endangered loggerhead sea turtle (Caretta caretta) revealed by 454 amplicon sequencing
.
BMC Evol. Biol.
13
,
95
106
[PubMed]
9.
Zhai
T.
,
Yang
H.Q.
,
Zhang
R.C.
,
Fang
L.M.
,
Zhong
G.H.
and
Fang
S.G.
(
2017
)
Effects of population bottleneck and balancing selection on the Chinese alligator are revealed by locus-specific characterization of MHC genes
.
Sci. Rep.
7
,
5549
5556
[PubMed]
10.
Ansari
T.H.
,
Bertozzi
T.
,
Miller
R.D.
and
Gardner
M.G.
(
2015
)
MHC in a monogamous lizard-characterization of class I MHC genes in the Australian skink Tiliqua rugosa
.
Dev. Comp. Immunol.
53
,
320
327
11.
Santonastaso
T.
,
Lighten
J.
,
Oosterhout
C.V.
,
Jones
K.L.
,
Foufopoulos
J.
and
Anthony
N.M.
(
2017
)
The effects of historical fragmentation on major histocompatibility complex class II β and microsatellite variation in the Aegean island reptile, Podarcis erhardii
.
Ecol. Evol.
7
,
4568
4581
12.
Pearson
S.K.
,
Bull
C.M.
and
Gardner
M.G.
(
2017
)
Egernia stokesii, (gidgee skink) positively selected sites lack concordance with HLA peptide binding regions
.
Immunogenetics
69
,
49
61
[PubMed]
13.
Pearson
S.K.
,
Bull
C.M.
and
Gardner
M.G.
(
2018
)
Selection outweighs drift at a fine scale: lack of MHC differentiation within a family living lizard across geographically close but disconnected rocky outcrops
.
Mol. Ecol.
27
,
2204
2214
[PubMed]
14.
Hacking
J.
,
Bertozzi
T.
,
Moussalli
A.
,
Bradford
T.
and
Gardner
M.
(
2018
)
Characterisation of major histocompatibility complex class I transcripts in an Australian dragon lizard
.
Dev. Comp. Immunol.
84
,
164
171
15.
Brandley
M.C.
,
Young
R.L.
,
Warren
D.L.
,
Thompson
M.B.
and
Wagner
G.P.
(
2012
)
Uterine gene expression in the live-bearing lizard, Chalcides ocellatus, reveals convergence of squamate reptile and mammalian pregnancy mechanisms
.
Genome Biol. Evol.
4
,
394
411
[PubMed]
16.
Yuan
X.Y.
,
Zeng
X.M.
and
Guo
X.G.
(
2014
)
MHC class I exon 4 in the multiocellated racerunners (Eremias multiocellata): polymorphism, duplication and selection
.
Asian Herpetol. Res.
5
,
91
103
17.
Radwan
J.
,
Kuduk
K.
,
Levy
E.
,
Lebas
N.
and
Babik
W.
(
2014
)
Parasite load and MHC diversity in undisturbed and agriculturally modified habitats of the ornate dragon lizard
.
Mol. Ecol.
23
,
5966
5978
[PubMed]
18.
Bauwens
D.
and
Verheyen
R.F.
(
1987
)
Variation of reproductive traits in a population of the lizard Lacerta vivipara
.
Ecography
10
,
120
127
19.
Spellerberg
I.F.
(
1976
)
Adaptations of reptiles to cold
. In
Morphology and Biology of Reptiles
vol.
3
, (
d’a Bellaires
A.
and
Cox
C.B.
, eds), pp.
261
285
,
Linnean Society Symposium Series
20.
Zhao
E.M.
,
Zhao
K.T.
and
Zhou
K.Y.
(
1999
)
Fauna Sinica Reptilia
, vol.
2
, p.
251
,
Science Press
,
Beijing
21.
Odierna
G.
,
Aprea
G.
,
Capriglione
T.
and
Puky
M.
(
2004
)
Chromosomal evidence for the double origin of viviparity in the 21. European common lizard, Lacerta (Zootoca) Vivipara
.
Herpetol. J.
14
,
157
160
22.
Cornetti
L.
,
Menegon
M.
,
Giovine
G.
,
Heulin
B.
and
Vernesi
C.
(
2014
)
Mitochondrial and nuclear DNA survey of Zootoca vivipara across the eastern Italian alps: evolutionary eelationships, historical demography and conservation implications
.
PLoS ONE
9
,
e85912
[PubMed]
23.
Cornetti
L.
,
Belluardo
F.
,
Ghielmi
S.
,
Giovine
G.
,
Ficetola
G.F.
,
Bertorelle
G.
et al.
(
2015
)
Reproductive isolation between oviparous and viviparous lineages of the Eurasian common lizard Zootoca vivipara in a contact zone
.
Biol. J. Linn. Soc.
114
,
566
573
24.
Cornetti
L.
,
Ficetola
G.F.
,
Hoban
S.
and
Vernesi
C.
(
2015
)
Genetic and ecological data reveal species boundaries between viviparous and oviparous lizard lineages
.
Heredity
115
,
517
526
[PubMed]
25.
Dupoué
A.
,
Rutschmann
A.
,
Le Galliard
J.F.
,
Clobert
J.
,
Angelier
F.
,
Marciau
C.
et al.
(
2017
)
Shorter telomeres precede population extinction in wild lizards
.
Sci. Rep.
7
,
169
176
[PubMed]
26.
Cornetti
L.
,
Griffith
O.W.
,
Benazzo
A.
,
Panziera
A.
,
Whittington
C.M.
,
Thompson
M.B.
et al.
(
2017
)
Candidate genes involved in the evolution of viviparity:a RAD sequencing experiment in the lizard Zootoca vivipara (Squamata:Lacertidae)
.
Zool. J. Linn. Soc.
183
,
196
207
27.
Heulin
B.
,
Surget-groba
Y.
,
Guiller
A.
,
Guillaume
C.P.
and
Deunff
J.
(
1999
)
Comparisons of mitochondrial DNA (mtDNA) sequences (16s rRNA gene) between oviparous and viviparous strains of Lacerta vivipara: a preliminary study
.
Mol. Ecol.
8
,
1627
1631
[PubMed]
28.
Boudjemadi
K.
,
Martin
O.
,
Simon
J.C.
and
Estoup
A.
(
1999
)
Development and cross-species comparison of microsatellite markers in two lizard species, Lacerta vivipara and Podarcis muralis
.
Mol. Ecol.
8
,
518
520
[PubMed]
29.
Horreo
J.L.
,
Pelaez
M.L.
,
Suarez
T.
,
Breedveld
M.C.
,
Heulin
B.
,
Surget-Groba
Y.
et al.
(
2018
)
Phylogeography, evolutionary history and effects of glaciations in a species (Zootoca vivipara) inhabiting multiple biogeographic regions
.
J. Biogeograph.
45
,
1616
1627
30.
Takeuchi
H.
,
Takeuchi
M.
and
Hikida
T.
(
2013
)
Low genetic diversity in the Japanese population of Zootoca vivipara (Squamata: Lacertidae) revealed by mitochondrial DNA
.
Curr. Herpetol.
32
,
66
70
31.
Zhao
D.
,
Liu
Y.F.
,
Zhao
W.G.
and
liu
P.
(
2017
)
Genetic diversity of females and offsprings based on microsatellite marking in different populations of Zootoca Vivipara
.
Nat. Sci. J. Harbin Normal Univ.
33
,
50
53
32.
Liu
H.
,
Liu
Y.F.
,
Liu
P.
and
Zhao
W.G.
(
2018
)
Sequence polymorphism analysis of mitochondrial Cyt b gene and D-loop region in Zootoca vivipara
.
Heilongjiang Anim. Sci. Vet. Med.
14
,
200
202
33.
Horreo
J.L.
,
Pelaez
M.L.
,
Suarez
T.
and
Fitze1
P.S.
(
2018
)
Development and characterization of 79 nuclear markers amplifying in viviparous and oviparous clades of the European common lizard
.
Genetica
146
,
115
121
[PubMed]
34.
Braverman
I.
(
2016
)
Anticipating endangerment: the biopolitics of threatened species lists
.
BioSocietier
12
,
132
157
35.
Zhao
W.G.
(
2018
)
Rare and Endangered Wild Animals and Plants in Heihe. 1
, pp.
44
45
,
Science Press
,
Beijing
36.
Murphy
B.F.
,
Thompson
M.B.
and
Belov
K.
(
2009
)
Evolution of viviparity and the maternal immune system: major histocompatibility complex (MHC) class I genes in skinks
.
Orbit Univ. Sydney Undergrad. Res. J.
1
,
1
17
37.
Tamura
K.
,
Peterson
D.
,
Peterson
N.
,
Stecher
G.
,
Nei
M.
and
Kumar
S.
(
2011
)
MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods
.
Mol. Biol. Evol.
28
,
2731
2739
38.
Zhao
M.
,
Wang
Y.Z.
,
Shen
H.
,
Li
C.L.
,
Chen
C.
,
Luo
Z.H.
et al.
(
2013
)
Evolution by selection, recombination, and gene duplication in MHC class I genes of two Rhacophoridae species
.
BMC Evol. Biol.
13
,
113
[PubMed]
39.
Librado
P.
and
Rozas
J.
(
2009
)
DnaSP v5: a software for comprehensive analysis of DNA polymorphism data
.
Bioinformatics
25
,
1451
1452
[PubMed]
40.
Botstein
D.
,
White
R.
,
Skolnick
M.
and
Davis
R.W.
(
1980
)
Construction of a genetic linkage map in man using restriction fragment length polymorphisms
.
Am. J. Hum. Genet.
32
,
314
331
[PubMed]
41.
Hughes
A.L.
and
Yeager
M.
(
1998
)
Natural selection at major histocompatibility complex loci of vertebrates
.
Annu. Rev. Genet.
32
,
415
435
[PubMed]
42.
Murrell
B.
,
Wertheim
J.O.
,
Moola
S.
,
Weighill
T.
,
Scheffler
K.
and
Pond Kosakovsky
S.L.
(
2012
)
Detecting individual sites subject to episodic diversifying selection
.
PLoS Genet.
8
,
e1002764
[PubMed]
43.
Nei
M.
,
Gu
X.
and
Sitnikova
T.
(
1997
)
Evolution by the birth-and-death process in multigene families of the vertebrate immune system
.
Proc. Natl. Acad. Sci. U.S.A.
94
,
7799
7806
[PubMed]
44.
Nei
M.
(
1978
)
Estimation of average heterozygosity and genetic distance from a small number of individuals
.
Genetics
89
,
583
590
[PubMed]
45.
Halbert
N.D.
,
Grant
W.E.
and
Derr
J.N.
(
2005
)
Genetic and demographic consequences of importing animals into a small population: a simulation model of the Texas State Bison Herd (USA)
.
Ecol. Model.
181
,
263
276
46.
Thorpe
J.P.
(
1982
)
The molecular clock hypothesis: biochemical evolution, genetic differentiation and systematics
.
Annu. Rev. Ecol. Syst.
13
,
139
168
47.
Wright
S.
(
1965
)
The interpretation of population structure by F-statistics with special regard to systems of mating
.
Evolution
19
,
395
420
48.
Slatkin
M.
(
1985
)
Rare alleles as indicators of gene flow
.
Evolution
39
,
53
65
[PubMed]
49.
Zheng
X.Z.
,
Xu
H.F.
and
Lu
H.J.
(
1997
)
Advances in research on genetic heterogeneity of animal populations
.
Biodiversity
5
,
210
216
50.
Takahata
N.
(
1990
)
Allelic genealogy under overdominant and frequency-dependent selection and polymorphism of major histocompatibility complex loci
.
Genetics
124
,
967
978
[PubMed]
51.
Bos
D.H.
and
Waldman
B.
(
2006
)
Evolution by recombination and trans-species polymorphism in the MHC class I gene of Xenopus laevis
.
Mol. Biol. Evol.
23
,
137
143
[PubMed]
52.
Bahr
A.
and
Wilson
A.B.
(
2012
)
The evolution of MHC diversity:evidence of intralocus gene conversion and recombination in a single-locus system
.
Gene
497
,
52
57
[PubMed]
53.
Miller
H.C.
and
Lambert
D.M.
(
2004
)
Gene duplication and gene conversion in class II MHC genes of New Zealand Robins (Petroicidae)
.
Immunogenetics
56
,
178
191
[PubMed]
54.
Hughes
A.L.
,
Hughes
M.K.
and
Watkins
D.I.
(
1993
)
Contrasting roles of interallelic recombination at the HLA-A and HLA-B loci
.
Genetics
133
,
669
680
[PubMed]
55.
Yeager
M.
and
Hughes
A.L.
(
1999
)
Evolution of the mammalian MHC: natural selection, recombination, and convergent evolution
.
Immunol. Rev.
167
,
45
58
[PubMed]
56.
Radtkey
R.R.
,
Becker
B.
,
Miller
R.D.
,
Riblet
R.
and
Case
T.J.
(
1996
)
Variation and evolution of class I MHC in sexual and parthenogenetic Geckos
.
Proc. Biol. Sci.
263
,
1023
1032
[PubMed]
57.
Guo
X.G.
,
Dai
X.
,
Chen
D.L.
,
Papenfuss
T.J.
,
Ananjeva
N.B.
,
Melnikov
D.A.
et al.
(
2011
)
Phylogeny and divergence times of some racerunner lizards (Lacertidae: Eremias) inferred from mitochondrial 16SrRNA gene segments
.
Mol. Phylogenet. Evol.
61
,
400
412
58.
Zhao
W.G.
(
2008
)
The Amphibia and Reptilia Fauna of Heilongjiang Province. 1
, pp.
133
147
,
Science Press
,
Beijing
59.
Chen
F.F.
,
Pan
L.
,
Geng
Z.Y.
,
Liu
X.L.
and
Yu
W.Y.
(
2010
)
Origin, evolvement and resistance mechanism of polymorphism of MHC molecules
.
Chinese J. Anim. Vet. Sci.
41
,
1061
1067

Author notes

*

These authors contributed equally to this work and should be considered as co-first authors.

This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY).