Patterns of evolution of MHC class II genes of crows (Corvus) suggest trans-species polymorphism

Sample collection

We analyzed the data from 18 individuals of each species. American crows were sampled from nestlings (one bird per nest) in Yolo County, California (38°32′N, 121°45′W) on the campus of the University of California, Davis, in May and June 2012. Collection methods were approved by the Institutional Animal Care and Use Committee of the University of California, Davis (Permit Number: 16897). Approximately 50 µL of whole blood was preserved in Queens’s lysis buffer, and gDNA was extracted using a standard phenol/chloroform protocol followed by ethanol precipitation. Previously extracted gDNA from the blood of jungle crows and carrion crows was obtained from the Japanese National Museum of Nature and Science (Tsukuba, Japan) and the Yamashina Institute of Ornithology (Chiba, Japan). These samples were collected between 1994 and 2010 from the entire range of the Japanese archipelago, although most jungle crows were from the Tokyo area (16 from Tokyo, 1 each from Hokkaido and Okinawa), whereas carrion crow collection was more evenly distributed across the Japanese archipelago (6 from Kagoshima, 1 from Tsushima, 4 from Nara, 1 from Osaka and 6 from Tokyo; Fig. S1). One mL of whole blood for RNA extraction was collected from three juvenile jungle crows trapped at Ueno Zoo in Tokyo, Japan in 2012. The blood was flash frozen on site and stored at −80 °C. Jungle crow RNA was extracted from whole blood using the RNeasy Protect Animal Blood System (Qiagen, Hilden, Germany).

Amplification of MHC class IIB exon 2

In order to isolate MHC IIB exon 2, which codes for the PBR, jungle crow cDNA from two individuals (U2 and U3) was synthesized with the LongRange 2Step RT-PCR System (Qiagen, Hilden, Germany) using the manufacturer’s recommended amplification protocol. MHC IIB was then targeted in the cDNA using the primers MHC05 (Miller & Lambert, 2004), which anneals within exon 1 in several passerines (positions variable), and ComaIIbex3R, which was designed by aligning the conserved regions of IIB exon 3 from several passerines (positions 19–38). These primers amplified a 359 bp cDNA fragment. We verified that the amplicon was MHC IIB by cloning the fragment from both individuals (TOPO XL One Shot; Invitrogen, Carlsbad, California, USA). Twenty-four clones from each individual were selected and whole colony PCR was performed using the TOPO vector primers M13 F and M13 R. The 20 µL PCR contained a pluck of colony cells as template, 0.625 µM of each primer, 0.5 mM of each dNTP, 1.5 mM MgCl2, 1X PCR buffer and 0.5 U of ExTaq (Takara, Mountain View, California, USA) polymerase. The cycling conditions included an initial denaturation at 94 °C for 2 min followed by 30 cycles of 94, 55 and 72 °C each for 1 min and a final extension of 72 °C for 10 min. We sequenced each amplicon in both directions using M13 primers (Applied Biosystems 3130xl; Applied Biosystems, Carlsbad, Calfornia, USA). We performed a BLAST search and identified eight sequences that aligned with known passerine MHC IIB.

We targeted intron 1 by aligning the cDNA-derived amplicons using MUSCLE in Geneious 6.1.4 (Biomatters, San Francisco, California, USA) and a reverse primer, ComaIIbex2R, was designed in a conserved region of exon 2 and paired with MHC05. Using gDNA from two jungle crows (Coma5061 and Coma7509), a 516 bp fragment extending from exon 1 to exon 2 was amplified using 50–100 ng of gDNA, 1.0 µM of each primer, 0.5 mM of each dNTP, 1.5 mM MgCl2, 1X PCR buffer and 0.5 U of ExTaq (Takara, Mountain View, California, USA) polymerase (25 µL total volume). The cycling conditions included an initial denaturation at 94 °C for 30 s followed by 30 cycles of 94 (20 s), 60 (20 s) and 72 (45 s) °C and a final extension of 72 °C for 10 min. Each of these two amplicons was cloned and 10 colonies were sequenced in each direction using the cloning and whole colony PCR protocol described above. Sequences were aligned using MUSCLE and a forward primer flanking exon 2, ComaiF2, was designed for the (100%) conserved region at position 282–304 of intron 1. Next, we isolated intron 2 by pairing ComaiF2 with ComaIIBex3R using the same individuals and PCR, cloning, and whole colony PCR conditions described above. When amplifying from exon 2 to exon 3 in jungle crows, we discovered indels that resulted in large gaps at different locations in intron 2. 454 pyrosequencing is most efficient when sequences are of equal length, and for this reason, we designed a single degenerate reverse primer, ComaEx2RA, between positions 249–268 of exon 2. Finally, we confirmed that the primers ComaiF2 and ComaEx2RA amplified MHC IIB loci in all three crow species. Primer sequences and annealing temperatures are listed in Table S1.

454 pyrosequencing

Fusion primers were synthesized (Fasmac, Atsugi City, Kanagawa, Japan) by ligating the standard Roche multiplex identifiers (MIDs) to both the forward and reverse MHC IIB primers and the 454 adaptor primers (the Roche “Basic Amplicon” design). In order to minimize the formation of PCR artifacts, we reduced the cycle number and primer concentration and eliminated the final extension step during amplicon generation (Lenz & Becker, 2008). PCRs contained 25–50 ngs of template, 0.5 µM of each fusion primer, 0.5 mM of each dNTP, 1.5 mM MgCl2, 1X PCR buffer, 5% DMSO and 0.5 U of ExTaq (Takara, Mountain View, California, USA) polymerase (25 µL total volume). The cycling conditions included an initial denaturation at 94 °C for 30 s followed by 28 cycles of 94 (20 s), 60 (20 s) and 72 (45 s) °C.

PCR amplicons were purified using the Agencourt AMPure XP System (Beckman Coulter, Brea, California, USA) and recovered DNA was checked for quality and primer dimers on a 1.5% agarose gel. Samples were then quantified using the Quant-iT PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, California, USA) and pooled in equimolar amounts according to the Roche GS Junior Titanium Amplicon Library Preparation Method Manual. The quality and concentration of the pooled samples was checked on a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California, USA) and then sequenced using the Roche GS Junior Titanium Sequencing System at the Japanese National Museum of Nature and Science in Tsukuba, Japan.

Amplicons were bi-directionally sequenced, and reads that passed the initial 454 Roche Junior quality filter were de-multiplexed using jMHC (Stuglik, Radwan & Babik, 2011). jMHC is especially efficient for de-multiplexing bi-directionally sequenced double-tagged amplicons because it bins only those sequences that contain both forward and reverse primers and their associated MIDs with no ambiguous characters (“N”s) in the primers, MIDs or target sequence. Sorted FASTA files (excluding primers and MIDs) were then imported into Geneious and aligned with published MHC class IIB sequences. Sequences containing indels were removed from the data set at this time.

Bioinformatics variant validation and genotyping

Variants were accepted as “true” if they occurred in at least three reads (r > 2) in both of two independent PCRs of the same individual. The criteria of (r > 2) is derived from the estimate of a read containing at least one sequence error for a 171 bp fragment of (HLA) DRB exon 2 of ∼0.11 (Galan et al., 2010). Thus, the probability that an identical error occurs three times in the same sample is ∼10⁻⁸ (Galan et al., 2010). At this stage, the most likely types of artefacts remaining in the data set are either single base substitutions or fragment chimeras generated during PCR. Although both of these artifacts are randomly generated, if these errors occur early in the PCR, they could be replicated and represented at relatively high frequencies in individual amplicons. To further reduce the probability of erroneously including single base substitutions, we classified sequences as true variants only when they differed by at least two nucleotide substitutions from more common variants. Galan et al. (2010) reported an average chimera frequency of 0.06 and some chimeras occurred at frequencies >0.1; thus, elimination of chimeras using frequency thresholds may erroneously inflate variant estimates. It is highly unlikely that two independent PCRs would each generate identical chimeras (where as many as 10 different variants could recombine) or single base substitution errors three times (r > 2). By validating variants with a second PCR for each individual (rather than a second PCR within or across individuals) chimeras are eliminated, yet variants that are rare in the population are not erroneously eliminated from the data set. The total number of reads that were included at this stage was then designated as “net reads” for the purposes of genotype confidence level calculations.

We used the program ‘Negative Multinomial’ to calculate the number of sequences that are necessary for amplifying all the variants at least three times with a confidence level of 99% (Galan et al., 2010). While the program is limited to eight possible variants, the relationship between n, net reads, r, the number of times each variant must be sampled, and m, the maximum number of variants per amplicon, is linear. Our cloned sequences indicated a maximum of 20 variants per individual which corresponds to a minimum of 175 reads for 99.9% confidence; however, this assumes equal sampling of variants. Because amplification bias is likely when co-amplifying 10 loci, we expected this bias to be reflected in over- and under-representation of some variants (Burri et al., 2014). Therefore, we increased the minimum required net reads for genotyping to 275, which represents a frequency of 1.1% (amplification bias of ∼ one order of magnitude) for a variant that has exactly three copies in an amplicon represented by 275 net reads.

Tests for selection and MHC IIB phylogenetics

For each species, we generated an alignment of validated variants, and identified the codons that have been shown to comprise the PBR of HLA in humans (Brown et al., 1993). While most MHC studies rely on the HLA PBR codons identified in humans to estimate selection, it has not been shown that these same codon positions apply universally or that selection is limited to peptide binding codons. Therefore, in order to identify specific codons that are likely under selection in crows, we performed a Wu–Kabat analysis on the amino acid alignment of verified alleles in all three species. A Wu–Kabat plot predicts amino acids that are likely subject to selection by identifying positions of high variability. Variability is calculated by dividing the number of amino acids at a position by the frequency of the most common amino acid (Wu & Kabat, 1970). For each species, we identified codons/amino acids as likely under selection if the Wu–Kabat variability metric was ⩾ twice the mean of Wu–Kabat values for all sites (Wu & Kabat, 1970; Bos & Waldman, 2006). These sites are denoted WuK-PS (PS = polymorphic site). We tested for selection on the HLA PBR sites by calculating the ratios of non-synonymous (dN) and synonymous mutations (dS) using a Z test (modified Nei-Gojobori-Jukes Cantor correction) in MEGA v. 6.0 (Tamura et al., 2011). For each species, individual codons were tested for selection using the HYPHY package in MEGA. This method, known as the counting method, estimates the number of nonsynonymous and synonymous mutations that have occurred at each codon by using a maximum likelihood reconstruction of the ancestral state of each sequence. The test statistic dN–dS is used for detecting codons that have undergone selection and for positive values the probability of rejecting the null hypothesis of neutral evolution (p-value) is calculated (Suzuki & Gojobori, 1999; Pond & Frost, 2005). Significant p values are denoted as positively selected sites. This method is more computationally compatible with large data sets than empirical or hierarchal Bayesian approaches (Yang et al., 2000), yet has been shown to provide nearly identical results with a sufficient number of sequences (Pond & Frost, 2005).

We constructed nucleotide phylogenies from five different partitions of exon 2: variable sites identified by the Wu-Kabot plot (WuK-PS, 33bp), non-WuK-PS sites (213 bp), HLA PBR sites (66 bp), non-HLA PBR sites (180 bp) and the contiguous 246 bp fragment of exon 2. For each analysis, 24 different nucleotide substitution models were tested for best fit (ML) in MEGA and the Jukes-Cantor (JC) Model + G (Gamma) + I (invariant sites) was chosen based on the Bayesian Information Criterion (BIC) score.

DNA phylogenetics and species divergence

To ensure that the mtDNA phylogenies used in Haring et al. (2012) were consistent with our populations, we tested a subset of carrion crows (N = 4), jungle crows (N = 6) and American crows (N = 4) at the same mtDNA region using the primers CR-Cor+ and Phe-Cor- (Table S1) and the same PCR conditions as Haring et al. (2012), followed by Sanger sequencing as described above. We aligned these sequences with four mtDNA sequences from each species that were used in the original phylogenetic construction by Haring et al. (2012). The Hasegawa-Kishino-Yano (HKY) model (Lenz et al., 2013) had the lowest BIC scores in the best-fit substitution model (ML). HKY ML trees with 500 bootstrap replicates were constructed in MEGA v. 6.0.

Using the mtDNA sequences, we estimated divergence times for the three crow species: the divergence of carrion crow/American crow clade from jungle crow, and the divergence between carrion crows and American crows. Divergence times to the most recent common ancestor were estimated in BEAST v. 2.0 (Bouckaert et al., 2014). We used the HKY site model and a Yule process speciation prior for the branching rates. We tested the assumption of a molecular clock using the maximum likelihood method which compares the ML value for the given topology with and without the molecular clock constraints under the HKY model (MEGA v. 6.0). We ran two separate analyses using two different mtDNA substitution rates (per site, per My) previously estimated for crows (C. macrorynchos = 0.0567555 and C. coronoides = 0.0603431) (Nabholz, Glemin & Galtier, 2009). The Markov chain Monte Carlo (MCMC) analyses were run for 10⁷ generations (10,000 iteration burn in); the mean and 95% highest posterior density interval (HPD) for divergence times were calculated in Tracer v. 1.6 (Bouckaert et al., 2014).

Pyrosequencing

The goal of the first 454 run was to estimate the number of IIB loci and the minimum read number required for accurate genotyping (99.9% confidence). Six individuals from each species were sequenced and a total of 125,839 reads passed the initial 454 filter. After excluding reads of less than 270 bp, there were 76,891 reads with an average length of 340 nucleotides. Amplicons had between 8,384 and 483 initial reads (mean = 1,213) that were sorted by jMHC. These samples had an average of 376 net reads (range 171–3,211). The maximum number of variants (differing by at least 2 nucleotides) found in a single individual was 20, which suggests that crows have up to 10 IIB loci per haplotype. The required minimum read number was calculated as 211; however, we increased the cutoff to 275 net reads to account for amplification bias (see methods).

For the second 454 run, a total of 107,279 reads passed the initial 454 filter. After excluding reads of less than 270 bp, there were 90,203 reads with an average length of 340 nucleotides. A total of 18 individuals (with replicates) from each species met the criteria for variant validation of at least 275 net reads and at least three identical reads per variant in both of two independent PCRs of the same individual. All DNA sequences were deposited in NCBI Genbank: Accession numbers: MHC IIB: TBA; crow mtDNA: KM246294–KM246308; crow nuclear intronic: KM246309–KM246323

MHCIIB variation and selection

The 248 bp fragment was conserved at the final two base pair positions in all three species and was trimmed to 246 bp. The number of validated MHC IIB nucleotide variants for each species was: jungle crows = 89, carrion crows = 81 and American crows = 67. Individuals had a range of 7–20 variants, indicating that the MHC IIB in crows exhibits copy number variation (Table 1), although it is possible that the range of variant number could be explained by other phenomena, including null alleles, allele sharing across loci, homozygous vs. heterozygous loci or the variant validation methodology (where sequences with indels were excluded from the data set). Few variants were shared among species and there was no discernable pattern of allele sharing relating to phylogeny or sympatry: Four unique variants were shared among American crows and jungle crows, three among American crows and carrion crows, three among carrion crows and jungle crows and four variants were shared by all three species.

The codon-based Z test for selection (non-synonymous/synonymous mutations) averaged across the PBR yielded a significant p value for American crows and carrion crows but not for jungle crows (Table 1). The Wu-Kabot plot identified 11 amino acid sites that were highly variable in all three species: 6, 23, 25, 32, 33, 52, 61, 62, 65, 73 and 81. Eight of these sites (73%) corresponded to known peptide binding sites in HLA (Fig. 1). Among the three species, jungle crows exhibited the most amino acid variation, with five hypervariable (>20 Wu–Kabat index) amino acid sites (Fig. 1). The HYPHY test for selection identified nine positively selected sites that were shared among the three species (2, 4, 33, 52, 55, 66, 68, 72 and 81), but only 5 of these (56%) corresponded to the HLA PBR (Fig. 1).

In total, we recovered 172 unique amino acid sequences among the three species (70 in jungle crows, 58 in carrion crows and 44 in American crows). We detected no species-specific motifs or motifs that were shared exclusively between two species. Within-species amino acid diversity was relatively high (0.245 ± 1.5), similar to that which we calculated for scrub jays (Aphelocoma coerulescens), another corvid species (23.3%; N = 11; NCBI (U23958-65; 72,73,75) (Edwards, Grahn & Potts, 1995). An alignment (82 sites) of common amino acid fragments from the three crow species with four other passerine species, and a falconiform species for comparison, highlights the high level of amino acid diversity found among and between species of passerines (Fig. S2).

MHC IIB phylogenies

We observed a similar pattern among the five different partitions of exon 2: MHC IIB variants formed several supported (bootstrap values ≥50), interspecific clusters in trees comparing putatively functional sites. Interspecific clusters were also observed when comparing putatively neutral sites, but there were fewer of these supported clusters, and in general they were less inclusive. This pattern was most striking in trees constructed using the contiguous 246 bp exon 2 fragment (Fig. 3). For instance, a tree that was constructed by comparing non-synonymous substitutions (Fig. 3A) shows six well supported clusters, each containing all three crow species, whereas a tree constructed by comparing synonymous substitutions of the same fragment (Fig. 3B), shows only three such clusters. Supported clustering at putatively functional sites was also more inclusive in these exon 2 trees, with 44% of all variants (105/237) forming supported interspecific clades compared to just 18% (32/237) in the tree using synonymous sites. The trees constructed from WuK-PS and non-Wuk-PS sites (Figs. S3A and S3B) as well as those constructed from the HLA PBR and non-PBR sites (Figs. S4A and S4B) displayed similar, though less distinct patterns. Among all partitions, the most common (supported) clustering was among all three species. We did not observe biased clustering of variants among the sympatric carrion crows and jungle crows; or among the more recently diverged, allopatric, American crows and carrion crows.

Supertyping

Among the three species, there were nine shared positively selected sites in exon 2. Five of these nine sites corresponded to HLA PBR sites (Fig. 1). From the 140 unique positively selected sites amino acid variants a total of eight supertypes were identified (Fig. 4). Only one supertype (Supertype 4) was limited to carrion crows and jungle crows, whereas the rest of the supertypes were shared among the three species. The supertypes generated by the WuK-PS and the HLA PBR codons were not well supported and thus excluded from our analysis.

mtDNA phylogeny and divergence estimate

The CR-Cor+ and Phe-Cor- primers amplified an 836–920 bp fragment of the mtDNA CR region in the three species of crows. The ML tree placed carrion crows and American crows in a separate clade from jungle crows with 100% bootstrap support (Fig. 2), thus confirming that the jungle crows and carrion crow populations from this study conform to the phylogeny of Haring et al. (2012). The ML-based test of the molecular clock failed to reject the null hypothesis of equal evolutionary rate throughout the tree (p = 0.85); thus we used a strict clock setting for the MCMC iterations in the BEAST program. Using the mtDNA substitution rate of 0.0567555/site/My the estimated mean divergence time to MRCA of all species was 0.8992 Ma with a 95% HPD interval of 0.7055, 1.1147 Ma (effective sample size (ESS) = 8581.7, median = 0.8942, standard deviation = 0.107). The estimated divergence to the MRCA of carrion crows and American crows was 0.4372 Ma with a 95% HPD interval of 0.3036, 0.5801 Ma (ESS = 7682.7, median = 0.4319, standard deviation = 0.0723) (Fig. 2). Divergence estimates using the higher substitution rate (C. coronoides = 0.0603431/site/My) yielded similar values: the MRCA of all species = 0.8509 Ma, HPD interval = 0.6578, 1.0486 Ma; MRCA of carrion crows and American crows = 0.4143 Ma, HPD interval = 0.2881, 0.5458 Ma.