Genetic variation and selection in the major histocompatibility complex Class II gene in the Guizhou pony

Animal collection and DNA isolation

A total of 50 blood samples were collected from GZP in Ziyun County, Anshun, Guizhou Province, China. All ponies used in our study were 4 to 8 years old. All animal procedures were approved by the Institutional Animal Care and Use Committee of Guizhou University (Approval number EAE-GZU-2018-P007). The GZP were randomly selected and were all well-developed and in good health, with heights ranging from 102 to 118 cm and weights between 210 to 265 kg. Blood samples were collected from the jugular vein and were kept in EDTA Na2. All samples were stored at −20 °C until DNA extraction. Genomic DNA was extracted from blood samples using the SQ Blood DNA Kit (OMEGA, USA). The nucleic acid concentration of the extracted genomic DNA was calculated by determining OD260/OD280, and detected by 0.7% agarose gel electrophoresis.

PCR amplification, cloning, and sequencing

The exon 2 regions of the ELA-DQA, DQB, DRA, and DRB genes were amplified from genomic DNA using PCR with specific primers. We amplified 246 bp of the DRA using the equid-specific primers DRA-F and DRA-R (Albright-Fraser et al., 1996), 246 bp of the DQA using the primers DQA-F and DQA-R (Fraser & Bailey, 1998), 276 bp of the DRB using the primers DRB-F and DRB-R (Fraser & Bailey, 1996), and 230 bp of the DQB using the primers DQB-F and DQB-R (Mashima, 2003). All primers were synthesized by the Bio-Engineering Company (Shanghai, China) (Table 1). The total PCR volume was 20 µL, and contained 10 µL of 2 × PCR Mixture (0.1 U Taq Plus Polymerase/µL, 500 µM dNTP each, 20 mM Tris–HCl (pH8.3), 100 mM KCl, 3 mM MgCl₂), 0.4 µL of upstream/downstream primers (10 µmol/L), and 1 µL templates. PCR amplification was carried out with initial denaturation at 95 °C for 5 min, followed by 30 cycles (95 °C for 30 s, 58 °C for 30 s, and 72 °C for 30 s), and a final extension at 72 °C for 10 min. PCR products were extracted and purified using the Gel Extraction Kit (OMEGA, USA), and were ligated into pGEM^®-T vectors and transformed into E. coli competent cells. Twenty positive clones of each sample were removed with a sterile toothpick and were detected using the Sanger sequencing method (Invitrogen, China). Alleles were confirmed if the same allele was observed in at least two different individuals.

DNA sequence polymorphism analysis

The base composition of the DRA, DRB, DQA and DQB genes was analyzed using MEGA7 software (Kumar, Stecher & Tamura, 2016). Standard descriptive diversity indices for each locus within the GZP were calculated using MEGA7 software, including the variable sites (V), parsim-info sites (P), singleton sites (S), and the transition/transversion bias ratio (R). It was important to ascertain whether the variability was uniformly distributed or was confined to small segments of the variable regions when determining the nature of the variable region. The Wu–Kabat variability index was calculated using the formula by Wu and Kabat Wu & Kabat (1970) with respect to amino acids at peptide-binding pockets. The variation of amino acids was calculated by the mutation rate (variability = number of different amino acids at a certain position/frequency of the most common amino acids at this position) (Wu & Kabat, 1970). Selection was estimated using MEGA7 software in terms of the relative rates of nonsynonymous (d_N) and synonymous (d_S) mutations, according to Nei and Gojobori’s method with the Jukes and Cantor (JC) correction (Nei & Gojobori, 1986). The selection Z-Test (P < 0.05) was performed for all sites under the null hypothesis of neutrality (d_N = d_S) and the alternative hypotheses of non-neutrality (d_N ≠ d_S), positive selection (d_N >  d_S), and purifying selection (d_N <  d_S).

Site-specific selection analyses and protein 3D structure analysis

We estimated the nonsynonymous and synonymous substitutions in the overall domain, ABS, and non-ABS for the DQA, DQB, DRA and DRB alleles. We assessed the positive selection using CodeML subroutine in the PAML program (Yang, 2007), which was more sensitive than other methods for assessing selection at the molecular level (Anisimova, Bielawski & Yang, 2001). The PAML program used the maximum likelihood estimation to examine heterogeneity in rates of ω = d_N/d_S among codons (Bielawski & Yang, 2003). The PAML program was able to better detect the molecular evidence of selection compared to other programs (Anisimova, Nielsen & Yang, 2003). We assessed heterogeneity in ω (ω  <  1: purifying selection, ω = 1: neutral evolution, ω >  1: positive selection) across the four alleles (DQA, DQB, DRA and DRB) to identify codons under positive selection. The observed ω value followed six models in PAML: M0 (one ratio, average ω across all sites), M1a (nearly neutral), M2a (positive selection), M3 (discrete), M7 (beta), and M8 (beta and omega) (Yang et al., 2000). We used the online SWISS-MODEL program (Biasini et al., 2014; Waterhouse et al., 2018) (https://swissmodel.expasy.org/interactive) to make predictions about the DQA, DQB, DRA and DRB protein structures.

Phylogenetic allele networks

We constructed a median-joining haplotype network to infer the phylogenetic relationships among the sequence haplotypes (Bandelt, Forster & Röhl, 1999) using the maximum parsimony in Network 4.6.1 (http://www.fluxus-engineering.com/sharenet.htm). The haplotype median networks of DQA, DQB, DRA and DRB between GZP and known horse species (Eqca, E.callabus; Eqpr, E.przewalski; Eqki, E.kiang; Eqgr, E.grevyi; Eqas, E.asinus; Eqbu, E.burchelli; Eqze, E.zebra; Eqhe, E.hemionus) from GenBank were plotted using NetWork 4.6. The frequency information and population proportion of the alleles were incorporated into the visualization of the network. Sequences from the horse, including E. callabus, E. przewalskii, E. burchellii, E. asinus, were incorporated to evaluate the distance from the Guizhou pony’s haplotypes (Table 2).

Analysis of nucleotide diversity

184 alleles were identified from 1,000 sequencing clones, with 118 effective alleles selected from the total. Of the 118 alleles, there were 18 novel DQA alleles (GenBank accession number: MT304744 –MT304761), 38 new DQB alleles (MT304705 –MT304743), 22 new DRA alleles (MT304762 –MT304783) and 28 new DRB alleles (MT304784 –MT304811) (Table S1). The alignment results are listed in Table S1 for the effective alleles from DQA, DQB, DRA, DRB and the sequences of JQ254059, AF034122, AJ575295, and AF144564. A considerable sequence diversity within the genus was revealed based on the DQA, DQB, DRB alignment results. The nucleotide diversity in DRA was much lower than in DQA/B and DRB in GZP, which is comparable with the level of nucleotide diversity in DRA from other species in the Equus genus (Kamath & Getz, 2011). Within the GZP, the genetic diversity was much higher in DQA, DQB, and DRB than in DRA and the ratio (variable sites/length) was the lowest at the DRA locus (15.04%) and highest at the DQB locus (46.08%).

Analysis of nucleotide compositions

The GC contents of DQB and DRB were higher than those of DQA and DRA (Table S2). The content of G+C (48.10%) was slightly lower than that of A+T (51.90%) at DQA, and the content of G+C (48.10%) was lower than that of A+T (51.90%) at DRA, which revealed that DQA and DRA had lower GC percentages. However, the base composition of G+C (64.20%) was higher than that of A+T (35.80%) in DQB alleles, and the base composition of G+C (61.90%) was more than that of A+T (38.10%) in DRB alleles, revealing that DQB and DRB had a higher GC content. The R (transitions/transversions) was 1.357 and 2.241 in the DQA and DRA alleles, respectively. However, there was an R of 0.778 and 0.573 in the DQB and DRB alleles, respectively. Our results revealed that the DRA locus was more well-conserved than the other loci.

The amino acid composition analyses

We determined that the exon 2 of the DQA, DQB, DRA, and DRB nucleotide sequences encoded 82, 76, 81, and 79 amino acid sequences, respectively. The underlined residues in Fig. 1 indicated an assumed ABS, based on the HLA equivalents (Brown et al., 1988; Brown et al., 1993), and may contact the antigen peptides (Figs. 1A–1D). There were 38 (45.12%), 51 (67.10%), 16 (19.75%), and 49 sites (62.02%) that were variable in the predicted amino acid sites of the DQA, DQB, DRA and DRB of GZP populations, respectively. For the ABS, 17 of 21 sites (80.95%), 12 of 18 sites (66.67%), 5 of 20 sites (25.00%), and 13 of 14 sites (92.85%) were diverse at the DQA, DQA, DQB, DRA and DRB loci. The amino acid compositions of DQA, DQB, DRA, and DRB at ABS were calculated using MEGA 7 software (Fig. 2). There were more polar R-amino acids at the DQA locus (46.30%), and included Gly, Cys, Ser, Tyr, Thr, Asn, and Gln (Fig. 2). The non-polar R-amino acids at the DRA locus had the highest percentage (58.86%), and included Ala, Leu, Val, Trp, Ile, Phe, Pro, and Met (Fig. 2). The largest proportion of charged R-amino acids (30.70%) was located at DQB, and consisted of Arg, Lys, His, Glu, and Asp (Fig. 2).

Global selection analyses

The Wu-Kabat variability index was not used to select all of the variable amino acids (Wu & Kabat, 1970). A total of fifteen amino acids at the DQA locus were strongly selected at residues 10, 17, 18, 21, 23, 30, 46, 51, 52, 58, 60, 61, 62, 65 and 72, with the highest variability occurring at residue 60 (Fig. 3). Many polymorphic sites were observed at the DQB locus, with eight high mutation loci at residues at 16, 27, 38, 46, 47, 57, 61, and 65 (Fig. 3). 12 residues were found to be polymorphic at the DRA locus, including residues at 12, 14, 15, 19, 29, 39, 47, 49, 63, 64, 67 and 69 (Fig. 3). Amino acid residues at six different positions in the DRB locus had high values (more than 30) on the Wu-Kabat variability index, and the strongly selected amino acid regions were found at residues 1, 2, 4, 5, 6, 7, 8, 12, 19, 28, 36, 47, 48, 50, 56, 58, 61, 62, 65 and 69 (Fig. 3). A comparison of the d_N/d_S ratio averaged across the whole coding region suggested that positive selection occurred at loci DQB (d_N/ d_S = 1.127, p = 0.322) and DRB (d_N/d_S = 1.228, p = 0.202), and purifying selection appeared at DQA (d_N/d_S = 0.779, p = 0.143) and DRA (d_N/d_S = 0.560, p = 0.069) (Table 3). All codon sites were not statistically significant according to the Z-tests (p > 0.05, Table 3). The estimates of d_N/d_S suggested that DQA and DRA were not affected by the positive selection at the genetic level.

Site-specific selection analyses

It is unlikely for selection to act uniformly across MHC genes over evolutionary time. Selection was more likely to occur at specific codons based on their functional role. The rate of nonsynonymous substitutions for the ABS (d_N = 0.594 ± 0.132) exceeded the number of synonymous substitutions four times (d_S = 0.128 ±  0.080) at the DRB (Table 3). Our results are in agreement with those observed in the Argentine Creole horse, which exhibited rates of nonsynonymous substitutions more than four times the number of synonymous substitutions at exon 2 of ELA-DRB (Díaz et al., 2001). The ABS rates of synonymous substitutions and nonsynonymous substitutions for DQA and DQB were similar (d_N = 0.330 ± 0.088, d_S = 0.287 ± 0.110; d_N = 0.206 ± 0.079, d_S = 0.133 ±  0.076, respectively) (Table 3). The ABS sites at the DRA exhibited less nonsynonymous substitutions (d_N = 0.017 ±  0.008) than synonymous substitutions (d_S = 0.028 ± 0.022) with a d_N/d_S ratio of 0.607 (Table 3). Z-tests performed separately on ABS were significant for DRB (p = 0.001) providing evidence for positive selection at these sites. We could not reject the null hypothesis of neutral evolution at the non-ABS site (Table 3). The Z-tests by site type at the DQA, DQB, and DRA sites could not reject the null hypothesis of neutrality (p > 0.05). In contrast, the non-ABS sites showed more synonymous substitutions than nonsynonymous substitutions with d_N/d_S ratios of 0.581, 0.895, 0.520, and 0.678 at DQA, DQB, DRA and DRB, respectively (Table 3).

The results from the selection analyses in PAML revealed different levels of selection for the four loci (Table 4). The variable evolutionary rates across the codon sites (M3) fit our data better than the M0 model and models M2a and M8 had higher log-likelihoods than positive selection (M1a and M7). The M2a and M8 models implied that approximately 2% of sites may be under positive selection at the DQA site (ω = 8.583, ω = 8.425) (Table 4). The posterior means of ω were estimated across the DQA codons under positive selection models and predicted fourteen sites (positions 10, 17, 18, 30, 46, 51, 52, 60, 61, 65, 68, 69, 70, 72) that may be under selection (ω > 1), nine (10, 30, 46, 51, 52, 61, 65, 68, and 72) of which were also putative ABS, based on the HLA equivalents (Fig. 4). However, the discrete model (M3: 3 discrete evolutionary rate classes) had the highest log-likelihood and estimated that only 6.6% of codon sites had ω values greater than one (ω = 10.031) with the remaining 93.4% of sites being assigned ω values close to 0 (Table 4).

However, the posterior means of ω across DQB codon sites estimated by M2a (ω = 6.240) and M8 (ω = 6.373) predicted that only five codons were under significant positive selection (positions 16, 27, 47, 57, 61). These five codons were also known as putative ABS based on the HLA equivalents (Fig. 4). The M3 model at the DQB estimated that approximately 8% of codon sites had ω values greater than one (7.3% with ω = 1.283; 0.6% with ω = 6.935) with the remaining 92% of sites being assigned ω values close to 0 (ω = 0.054) (Table 4).

The M3 model estimated that only 6.6% of codon sites had ω values greater than one (ω = 10.031) with the remaining 93.4% of sites being assigned ω values close to 0 for the DRA (Table 4). Moreover, the posterior means of ω across the DRA codon sites estimated by M2a (ω = 10.286) and M8 (ω = 10.323) predicted that nine codons (positions 12, 14, 15, 16, 18, 19, 49, 64, 68) were under significant positive selection. However, only two sites (19 and 49) were known as putative ABS based on the HLA equivalents (Fig. 4).

The M3 model estimated that 14.6% of codon sites had ω values greater than one (ω1 = 1.351, ω2 = 6.823) at the DRB (Table 4), which was higher than the other MHC class II codons (DQA, DQB and DRA). The posterior means of ω across the DRB codon sites were estimated by the M2a (ω = 5.972) and M8 (ω = 5.961) models and predicted that twelve codons (positions 1, 2, 23, 28, 47, 48, 58, 61, 62, 65, 69, and 77) were under significant positive selection, seven (2, 28, 48, 58, 61, 69, and 77) of which were also putative ABS based on the HLA equivalents (Fig. 4).

Evolutionary analysis

We could not determine the genealogy of DQA, DQB, DRA and DRB due to the presence of loops in the network (Figs. 5A–5D). Some alleles were more likely to be ancestral based on their internal position in the network and a greater frequency of mutational connections. These alleles seemed more likely to be ancestral at the DQA locus, including DQA1, DQA3, Eqca17, and Eqca18. Allele DQA1 appears to be ancestral for most alleles, namely DQA12, DQA13, DQA14, DQA15, DQA9, Eqca10, Eqbu6, Eqca20, Eqas1, Eqas2, Eqbu5, Eqca19, Eqbu4, and Eqbu12. Three haplotypes of Przewalski’s horse (Eqpr3, Eqpr4 and Eqpr2) were separated from DQA3 by two mutational steps and were most closely related to GZP haplotypes. Meanwhile, the DQA1 allele was shared among four species (Eqca15, Eqgr1, Eqbu2 and Eqze1), DQA2 was shared with Eqca14, and DQA3 was shared between two species (Eqca16 and Eqbu20). At the DQB locus, allele DQB1 appears to be ancestral for most alleles, including DQB31, DQB23, DQB32, DQB10, DQB30, DQB33, DQB29, DQB3, DQB24, DQB34, DQB40, Eqas7 and Eqas4. We found that haplotypes DQB5 and DQB13 were shared between Eqca1 and Eqca7, respectively. Allele DRA5 was shared between Eqbu7 and Eqca5, and DRA1 was shared with DRA3 for the DRA locus. Interestingly, allele DRA1 seemed more likely to be ancestral, containing twenty alleles, including DRA21, DRA15, DRA13, DRA8, DRA18, DRA7, DRA14, DRA6, DRA20, DRA19, DRA17, DRA16, DRA2, DRA11, DRA10, DRA12, Eqca2, Eqca6, Eqca7, and Eqca8. Haplotypes Eqhe and DRA9 were separated from DRA1 by as two mutations step as are most closely related GZP haplotypes. Allele DRA5 seemed more likely to be ancestral, Eqbu, Eqze, Eqgr, Egas were separated from DRA5 by one or two mutational step and are most closely related to the GZP haplotypes. Most DRB alleles were dispersed throughout the whole network, and there was a closer genetic relationship between GZP and other horse species. Wild ass haplotypes, Eqas3, Eqas4 and Eqas6, were separated from DRB28 by one mutational step, as are most closely-related GZP haplotypes. Furthermore, the haplotypes DRB2 (Eqpr1) and DRB3 (Eqpr2) were shared by GZP and Przewalski’s horse. The haplotypes DRB1 (Eqca5), DRB2 (Eqca12), DRB4 (Eqca7), DRB5 (Eqca1), DRB15 (Eqca2) and DRB23 (Eqca8) were shared between GZP and the European horse.