Genetic variation at the MHC DRB1 locus is similar across Gunnison's prairie dog (Cynomys gunnisoni) colonies regardless of plague history.

Yersinia pestis was introduced to North America around 1900 and leads to nearly 100% mortality in prairie dog (Cynomys spp.) colonies during epizootic events, which suggests this pathogen may exert a strong selective force. We characterized genetic diversity at an MHC class II locus (DRB1) in Gunnison's prairie dog (C. gunnisoni) and quantified population genetic structure at the DRB1 versus 12 microsatellite loci in three large Arizona colonies. Two colonies, Seligman (SE) and Espee Ranch (ES), have experienced multiple plague-related die-offs in recent years, whereas plague has never been documented at Aubrey Valley (AV). We found fairly low allelic diversity at the DRB1 locus, with one allele (DRB1*01) at high frequency (0.67-0.87) in all colonies. Two other DRB1 alleles appear to be trans-species polymorphisms shared with the black-tailed prairie dog (C. ludovicianus), indicating that these alleles have been maintained across evolutionary time frames. Estimates of genetic differentiation were generally lower at the MHC locus (F ST = 0.033) than at microsatellite markers (F ST = 0.098). The reduced differentiation at DRB1 may indicate that selection has been important for shaping variation at MHC loci, regardless of the presence or absence of plague in recent decades. However, genetic drift has probably also influenced the DRB1 locus because its level of differentiation was not different from that of microsatellites in an F ST outlier analysis. We then compared specific MHC alleles to plague survivorship in 60 C. gunnisoni that had been experimentally infected with Y. pestis. We found that survival was greater in individuals that carried at least one copy of the most common allele (DRB1*01) compared to those that did not (60% vs. 20%). Although the sample sizes of these two groups were unbalanced, this result suggests the possibility that this MHC class II locus, or a nearby linked gene, could play a role in plague survival.

Introduction

Infectious diseases can play an important role in shaping genetic variation and fitness in free‐living populations (Anderson and May 1979; Keller and Waller 2002; Keesing et al. 2010). For example, diseases may decrease allelic variation at immune system loci via selective sweeps (de Groot et al. 2002; Carnero‐Montoro et al. 2012) and reduce variation at neutral loci as a result of genetic bottlenecks (Maruyama and Fuerst 1985; Luikart and Cornuet 1998). Conversely, a cyclic pattern of rare allele advantage under pressure by coevolving pathogens (an evolutionary arms race) should result in the maintenance of a greater number of alleles over time due to balancing selection (Bodmer 1972; Potts and Wakeland 1990; Slade and McCallum 1992). Empirical evidence from a growing number of species suggests that increased genetic diversity at neutral and adaptive loci is associated with reduced susceptibility of individuals to infectious diseases and parasites (O'Brien and Evermann 1988; Siddle et al. 2007; Altermatt and Ebert 2008). Therefore, characterizing extant genetic diversity, particularly at functionally relevant immune system genes, can contribute to a more comprehensive understanding of the complex dynamics between genetic diversity and disease resistance in free‐living populations. Among populations of conservation concern, introduced diseases are considered one of the top extinction dangers (Altizer et al. 2003). Thus, understanding the impact of infectious diseases in wild systems is critical for species management and protecting human public health (Daszak 2000).

An introduced disease organism with significance for many rodents in western North America is Yersinia pestis, the causative agent of plague. This bacterial pathogen originated in Asia and was introduced to several U.S. ports from China around 1900, at which time it became established in commensal rodents and fleas in San Francisco (Link 1955; Gage and Kosoy 2005). The pathogen soon jumped to native rodents, especially colonial ground squirrels (Eskey and Haas 1940), and spread eastward from California. Plague began affecting prairie dogs (Cynomys spp.) in the 1930s and led to major population reductions in the five species found in western North America (Hoogland et al. 2004). Prairie dog decline has become a critical conservation concern for black‐footed ferrets (Mustela nigripes) because they are highly specialized predators of prairie dogs that are completely dependent on large persistent colonies for survival (Biggins et al. 2006).

Given the repeated occurrence of plague epizootics in the western United States and their high mortality in prairie dog populations, it is reasonable to infer that plague has been a strong selective force on prairie dogs and other rodents over the past 80 years (Thomas et al. 1988). All prairie dog species appear to be highly susceptible to Y. pestis, and mortality levels reach nearly 100% in colonies of Gunnison's prairie dog (C. gunnisoni) (Cully 1997; Cully et al. 1997), Utah prairie dog (C. parvidens) (Cully and Williams 2001), and black‐tailed prairie dog (C. ludovicianus) (Pauli et al. 2006; Cully et al. 2010). White‐tailed prairie dog (C. leucurus) colonies seem to experience a lower mortality rate (Anderson and Williams 1997; Griffin et al. 2010), possibly because they are less social than other prairie dogs and do not spread the pathogen as efficiently (Gasper and Watson 2001). Gunnison's prairie dog appears to be particularly susceptible to local colony extirpations resulting from plague (Eskey and Haas 1940; Girard et al. 2004), and up to 70% of monitored C. gunnisoni colonies in the state of Arizona were extirpated between the early 1980s and 2001 (Wagner et al. 2006). The range of C. gunnisoni lies in the Four Corners area, a region of the southwestern United States that has been a major focus of plague outbreaks in the past century (Ari et al. 2010).

Two recent experimental challenge studies in C. gunnisoni and C. ludovicianus suggest that some resistance to Y. pestis may exist (Rocke et al. 2012; Busch et al. 2013). Survival of wild‐caught C. ludovicianus collected in plague‐endemic areas of Colorado and Texas (50% and 60% survival, respectively) was much higher during exposure to the fully virulent Y. pestis strain CO92 than individuals originating from a colony in western South Dakota (5% survival) that had not previously experienced plague (Rocke et al. 2012). A similarly high survival rate (50–70%) was found during an experimental challenge with Y. pestis CO92 in C. gunnisoni from two populations in Arizona (Espee Ranch and Aubrey Valley) (Busch et al. 2013). A strong antibody response does not appear to be the primary mechanism of survival in either species and suggests the possibility of a genetic component for plague resistance.

Both the Aubrey Valley (AV) and Espee Ranch (ES) populations are closely monitored because they are release sites for black‐footed ferrets (Fig. 1). Plague is known to occur at ES, and the most recent outbreaks occurred in 2009 (Busch et al. 2013) and 2014 (T. Rocke, unpublished data). In contrast, the AV population has not experienced any observable mortality from plague, even during years when plague epizootics affected nearby colonies (Van Pelt 1995). In 1996, plague occurred within just 6–10 km of AV at an adjacent colony complex near Seligman, AZ (the SE population). Despite apparent differences in plague occurrence, both AV and ES prairie dogs showed higher survivorship than expected during experimental challenge (Busch et al. 2013), which suggests that many individuals in both colonies may have some level of resistance to the disease.

Figure 1

Map of the three Cynomys gunnisoni colonies from this study. (A) Location of the three study colonies; AV = Aubrey Valley; SE = Seligman; ES = Espee Ranch. (B) Inset map showing the study area in Arizona. The AV and ES colonies (66 km apart) are reintroduction sites for black‐footed ferrets (Mustela nigripes). Additional details of collection locations are available in a previous publication (Busch et al. 2011).

A recent study of plague in Madagascar has suggested that multiple MHC regions may be under selection in black rat (Rattus rattus) populations that are regularly exposed to endemic Y. pestis (Tollenaere et al. 2012). Class I and II MHC genes are among the most polymorphic genes in the vertebrate genome (Klein and Figueroa 1986; Edwards and Hedrick 1998), and variability in this region often reflects evolutionarily relevant processes (Sommer 2005). The second exon of the class II DRB gene often bears a signature of selection because it comprises a substantial portion of the antigen‐binding region of this molecule (Klein et al. 1993). Variation in MHC genes has been linked to individual fitness in the presence of parasites (Paterson et al. 1998; Thoss et al. 2011; Kloch et al. 2013), but little is known about the specific contribution made by the class II DRB locus during Y. pestis infection.

In this study, we investigated patterns of genetic variation in C. gunnisoni colonies where plague is known to occur repeatedly, ES and SE, and the AV colony, where plague appears to be absent. We compared genetic structure at microsatellite markers with that of the MHC class II DRB1 locus. Other studies have used comparisons of neutral markers versus MHC genes to detect recent selection acting on the MHC gene region (Boyce et al. 1997; Landry and Bernatchez 2001; Worley et al. 2006). Assuming that plague is a strong selective pressure on C. gunnisoni populations and that this disease has a similar effect on all colonies (e.g., nearly 100% mortality) when it is present, we hypothesized that this MHC locus would exhibit appreciably lower levels of genetic differentiation (F ST) than observed at microsatellite markers. Neutral genetic markers should show greater levels of genetic structure due to population fluctuations and the increased strength of genetic drift on small populations. Likewise, we predicted that the AV population will have greater variation at the MHC locus, because this population has not been subjected to the multiple plague events that could result in reduced diversity at other colonies (ES and SE). We also generated MHC genotypes of all C. gunnisoni individuals from the previous experimental plague challenge study (Busch et al. 2013) to explore the association of specific DRB1 alleles with individual survivorship.

Material and Methods

Sample collection and DNA extraction

We sampled three wild populations of C. gunnisoni in northern Arizona to evaluate genetic diversity at the MHC class II DRB1 locus and 12 microsatellite repeat markers. Multiple tissues such as spleen, liver, lung, muscle, skin, and blood were collected from 94 prairie dogs. Animals were collected during two separate sampling sessions (Appendix 1). In the first sampling session (2005–2006), 34 prairie dogs were sampled for a previous study (Busch et al. 2011) using animals shot by local hunters at Aubrey Valley (AV2006, n = 16) and a neighboring colony at Seligman, Arizona (SE, n = 18). In the second session (2009), prairie dogs were live‐trapped from two colonies where black‐footed ferrets have been released in Arizona: AV (n = 30) and Espee Ranch (ES) (n = 30). The 60 animals from 2009 are designated as AVC or ESC in this study because they were taken into captivity and used in a plague challenge experiment for C. gunnisoni (Busch et al. 2013). Important details of survivorship (Y. pestis dose, fate, sex, age, etc.) for these 60 individuals are provided in Appendix 1. All tissue samples, including those collected in the field, were quickly frozen in liquid nitrogen and stored at −80°C prior to DNA or RNA extraction. We also made use of skin samples from black‐tailed prairie dogs (C. ludovicianus) live‐trapped in Phillips County, Montana after a plague epizootic in 2008 (generously donated by Randy Matchett from the Charles M. Russell National Wildlife Refuge). The C. ludovicianus samples were used as a species comparison for our initial description of MHC DRB1 alleles from C. gunnisoni. All work was approved by the Northern Arizona University Institutional Animal Care and Use Committee under protocol 05‐005.

Spleen or muscle samples were used to extract genomic DNA (gDNA) from C. gunnisoni. Tissues were disrupted with bead milling for 20 sec following a previous protocol (Allender et al. 2004), and gDNA was extracted using the Qiagen DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA) according to the manufacturer's instructions. DNA samples were diluted to 20 ng/μL working stocks for downstream PCR methods.

MHC sequencing

We used two sequencing‐based approaches to initially isolate and characterize the MHC class II DRB1 gene, which has not yet been described for any prairie dog species. First, we used a traditional approach to amplify, clone, and sequence the exon 2 fragment of the DRB1 locus using HLA primers GH46F/GH50R (Erlich and Bugawan 1990) (Appendix 2). PCRs were carried out in 10 μL volumes containing the following reagents (given in final concentrations): 20–40 ng of gDNA template, 1X PCR buffer, 2 mM MgCl2, 0.2 mM dNTPs, 2 U Platinum® Taq DNA polymerase (Invitrogen, Carlsbad, CA), and 0.2 μM of each primer. PCRs were thermocycled according to the following conditions: 95°C for 10 min to release the Plat Taq ® antibody, followed by 35 cycles of 94°C for 15 sec, 65°C for 20 sec, 72°C for 30 sec, and a final extension step of 72°C for 5 min. PCR products were electrophoresed on a 2% agarose gel for 1 h at 100V. If amplification was faint for an individual, additional PCRs were obtained before cleaning with a QiaQuick PCR Purification Kit (Qiagen, Valencia, CA).

To determine the specific DRB1 alleles present in heterozygous individuals, purified PCR product was ligated into the pGEM‐T vector (Promega, Madison, WI) and transformed into Escherichia coli JM109 cells. Cells were plated onto ampicillin/IPTG/X‐Gal agar and cultured at 37°C for 18 h before picking positive clones (white or light blue in color). PCR screens were performed on 24 colonies using the same conditions as for the initial PCR, except that primers M13F and R were used at an annealing temperature of 55°C. Eighteen positive clones from each sample were sequenced to ensure detection of each individual allele (Babik 2010). Sequencing reactions were performed with M13 primers using a BigDye® v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA), and the products were separated on an ABI3130 automated sequencer. Sequences were aligned using the Lasergene software package (DNASTAR, Madison, WI).

Our second approach to MHC class II gene discovery in C. gunnisoni used a Roche 454 platform to pyrosequence PCR amplicons using GH46F/GH50R primers modified with sequencing adaptors and bar codes (Appendix 2). These tailed primers were used in a 10‐μL PCR with the following conditions: 40 ng gDNA, 1 U Platinum® Taq DNA polymerase, 1X PCR buffer, 2 mM MgCl2, 0.2 mM dNTPs, and 0.4 μM each primer. A stringent touchdown procedure was required to obtain clean amplification of the DRB1 gene. Thermocycler conditions started at 95°C for 10 min, followed by eight cycles of 94°C for 45 sec, 66°C for 30 sec, and 72°C for 30 sec. In the next ten cycles, we decreased the annealing temperature (T a) by 0.4°C per cycle (decrement from 66°C to 62.4°C). The T a was then kept at 62°C for the remaining 17 cycles (35 cycles total), with a final extension of 72°C for 10 min. We pooled equimolar amounts of these PCR products and shipped them to the Genome Resource Center at the University of Maryland for 454 sequencing. Raw 454 sequence reads for each prairie dog were sorted by bar code and clustered into unique alleles. These allele clusters were aligned in SeqMan Pro (DNASTAR, Madison, WI) to generate consensus sequences for each prairie dog.

Validation of DRB expression

It is important to validate the expression of newly described DRB1 alleles, especially in nonmodel organisms. Therefore, we used a Total RNA Purification Kit (Norgen Biotek, Thorold, ON) to extract RNA from approximately 40 mg of spleen tissue from a subset of C. gunnisoni (n = 10); these ten validation individuals are highlighted in the first column of Appendix 1 with a “V” notation. We treated the total RNA with amplification grade DNase I (Invitrogen, Grand Island, NY) to degrade possible residual gDNA that may have been present in the RNA extractions. To confirm the absence of gDNA, we diluted the RNA 1/10 with molecular grade H2O for PCR and used the exon 2 primers GH46F/GH50R (Appendix 2). Taq polymerase is incompatible with RNA templates, and therefore, PCR failure indicated the absence of gDNA. Positive control reactions with C. gunnisoni gDNA were included, and all reactions were electrophoresed on 2% agarose as described above to validate the degradation of DNA.

Undiluted total RNA samples were then used as templates to synthesize cDNA in 10 μL reactions using Smartscribe™ reverse transcriptase (Clonetech, Mountain View, CA) and two oligos: 3′ SMART CDS Primer II (5′‐AAGCAGTGGTATCAACGCAgtacTTTTTTTTTTTT TTTTTTTTTTTTTTTTTTVN‐3′) and SMART II G‐linker (5′‐AAGCAGTGGTATCAAC GCAgtacGCGGG‐3′). The incubation time for first strand synthesis was extended to 90 min, and the reaction product was diluted with 40‐μL molecular grade H2O. Resulting cDNA was quantified on a NanoDrop 8000 spectrophotometer (Thermo Scientific, Waltham, MA) and diluted to 20–40 ng/μL for PCR. Successful cDNA synthesis was confirmed using the exon 2 primer PD1DRB_F (5′‐CGTTTCCTGGAGCAAGTTTCACATG‐3′) and a reverse primer anchored in exon 3, DRBexon3‐R2 (5′‐AGACCAGGAGGTTGTGGTGCTG‐3′) (Appendix 2). PCRs were carried out in 10 μL volumes containing the following reagents (given in final concentrations): 20–40 ng of cDNA template, 1X PCR buffer, 2 mM MgCl2, 0.2 mM dNTPs, 0.2 μM of each primer, 1 U Platinum® Taq DNA polymerase, and 0.2 μM of each primer. PCRs were thermocycled according to the following conditions: 95°C for 10 min to release the Plat Taq ® antibody, followed by 35 cycles of 94°C for 15 sec, 65°C for 20 sec, and 72°C for 30 sec, with a final extension step of 72°C for 10 min. PCR products were electrophoresed on a 2% agarose gel for 1 h at 100V. Amplicons were cloned into the pGEM‐T vector (Promega) and sequenced as described above. Each unique allele that was validated at this stage was named according to conventions for nonmodel species (e.g., MhcCygu‐DRB1*01 and MhcCygu‐DRB1*02) (Klein et al. 1990).

In the process of sequencing exons 2 and 3 from cDNA, we discovered that certain alleles (Cygu‐DRB1*02, Cygu‐DRB1*05, and Cygu‐DRB1*06) had accumulated multiple single nucleotide polymorphisms (SNPs) in the exon 2 priming site for GH50R, which was originally designed for human DRB sequences (Appendix 3). The accumulation of up to six SNPs in these alleles created an amplification bias that led to allelic dropout in standard PCR. Therefore, we designed and validated two new reverse primers (CyguDRB1_exon2_R1: 5′‐CTCTCCTCTCCACAGTGAAGCTCTCA‐3′ and CyguDRB1_exon2_R2: 5′‐CTCTCCGCTCCACAGCGAAGCTCTCA‐3′) that, when combined and used with GH46F, will allow unbiased amplification of all six MhcCygu‐DRB1 alleles (PCR conditions in Appendix 2).

Quantitative PCR assays for MHC

To provide rapid and accurate identification of all known DRB1 alleles we discovered in C. gunnisoni, we designed six quantitative PCR (qPCR) assays. Briefly, the allele typing process requires three steps (all methods provided in Appendix 2):

In the first step, four assays are used to amplify all six known alleles from this study. Alleles Cygu‐DRB1*02, Cygu‐DRB1*03, and Cygu‐DRB1*04 can be unambiguously identified with these four qPCR assays. However, alleles Cygu‐DRB1*01, Cygu‐DRB1*05, and Cygu‐DRB1*06 require 1–2 additional assays to accurately score the allele.

Next, the subset of individuals that appear to have at least one Cygu‐DRB1*01 allele are tested with an assay to distinguish allele 01 from 05 or 06.

If needed, a third PCR with competitive forward primers is used to unambiguously resolve alleles 05 and 06.

The specific PCR outcomes expected for each diploid genotype (21 unique combinations) are described in Appendix 4. These assays were first validated on cDNA, then run against the gDNA of all 94 prairie dogs and compared to the genotypes scored from Sanger and 454 sequences (Appendix 1).

MHC analyses

To assess genetic diversity at the MHC‐DRB1 locus within these C. gunnisoni populations, we calculated allele frequencies, observed heterozygosity (H O), expected heterozygosity (H E), and tested for Hardy–Weinberg equilibrium using a Markov Chain implemented in FSTAT 2.9 (Goudet 1995, 2001). This same program was used to estimate the fixation index, F ST, for a global mean and all possible pairwise comparisons. Significance testing in FSTAT used 1000 permutations per pairwise comparison. Because we had two population samples from AV separated by three years, the 16 prairie dogs collected in 2006 (AV2006) were excluded in all calculations that produced a mean value (e.g., mean H O and global F ST).

The MHC‐DRB1 allele sequences found in C. gunnisoni were analyzed for evidence of selection acting on the MHC. To do so, our six MhcCygu‐DRB1 allele sequences were aligned with a human sequence, translated to amino acids (Fig. 2), and conserved sites and putative antigen‐binding sites (ABS) were inferred from the human sequence, in accordance with previous studies (Kaufman et al. 1994; Bondinas et al. 2007; Yuhki et al. 2008). We calculated rates of synonymous (d S) and nonsynonymous (d N) substitutions per site at 89 amino acid positions and the 21 putative ABS using the Nei‐Gojobori method and applying a Jukes–Cantor correction for multiple substitutions at a site (Nei and Gojobori 1986) in DnaSP version 5 (Librado and Rozas 2009). A Z‐test was computed in MEGA v. 6.0.6 to test for positive selection acting on the second exon of this DRB1 locus (Tamura et al. 2013). We also implemented maximum likelihood and Bayesian methods to test for positive selection using PAML 4.8 (Yang 1997, 2007), as inferring selection from distance‐based measures of d N/d S has been shown to be problematic (Kryazhimskiy and Plotkin 2008; Wolf et al. 2009). A full description of model parameters and tests is provided in Appendix 5.

Figure 2

Translated amino acid alignment of MHC‐ alleles from Gunnison's prairie dog, Cynomys gunnisoni (MhcCygu‐, alleles 01–06: GenBank accessions KU695893‐KU695896 and KR338362‐KR338363) and black‐tailed prairie dog, C. ludovicianus (MhcCylu‐, alleles 01–05: KR338364‐KR338368). Dots represent identity to the Cygu‐*01 amino acid sequence at the top. Dashes indicate missing data. Shaded positions show putative antigen‐binding sites (ABS) after (Yuhki et al. 2008). None of the amino acid sites are predicted to be under positive selection in either prairie dog species. Other ground squirrels are included in the alignment, including spotted suslik (Spermophilus suslicus) Spsu‐: EF569186 (Biedrzycka and Radwan 2008; Biedrzycka et al. 2011); Eastern woodchuck (Marmota monax) Mamo‐A: KJ675569 (Moreno‐Cugnon et al. 2015); European marmot (Marmota marmota) Mama‐: JQ837902 (Kuduk et al. 2012), Abert's squirrel (Sciurus aberti) Scab‐: AAA42356; and human (Homo sapiens) ‐: NP_002115.

To further explore the evolutionary history of these alleles, we included five MHC‐DRB1 alleles isolated from black‐tailed prairie dogs (C. ludovicianus). These MhcCylu‐DRB1 alleles were discovered by PCR amplicon sequencing, but their expression has not been validated at this time. All statistical analyses for the C. ludovicianus genetic data were performed as described above for C. gunnisoni. A maximum likelihood phylogeny was constructed for the DRB1 locus of both species of prairie dog using MEGA 6.0.6. The best nucleotide substitution model (AICc = 4081) was found to be Jukes–Cantor with a gamma distribution (JC + G). Confidence estimates for nodes on the tree were generated from 1000 bootstrap replicates.

Because recombination within sequences can lead to the false identification of positive selection (Anisimova et al. 2003), we tested our allele sequences for evidence of recombination in both prairie dog species using three commonly used measures of linkage disequilibrium (r 2, D', and G4) in the program PERMUTE (Wilson and McVean 2006).

Microsatellite marker comparisons

To establish a baseline level of genetic variation among the three prairie dog populations, we genotyped all individuals using a set of 12 microsatellite markers (Appendix 6) that were originally developed for black‐tailed prairie dogs (Jones et al. 2005) and ground squirrels (May et al. 1997; Stevens et al. 1997). All details of PCR amplification and multiplexing are provided in Appendix 6. Microsatellite alleles were electrophoresed on an AB3130 sequencer and scored with GENEMAPPER v4.0 software (Applied Biosystems, Foster City, CA). All individuals were run twice at all markers, and no genotyping errors were observed.

To assess neutral genetic structure of the study populations, we calculated H O and allele frequencies of our microsatellite data using the program GENALEX 6.4. An exact test for Hardy–Weinberg equilibrium using a Markov Chain was implemented in FSTAT 2.9. Global and pairwise estimates of F ST were obtained using 1000 permutations per pairwise comparison in FSTAT. A global estimate of F ST from the DRB1 MHC locus was compared to the global 95% confidence interval (CI) from microsatellite markers.

A second point of comparison was made between microsatellite loci and the DRB1 locus using the software fdist2 (Beaumont and Nichols 1996). First, a neutral distribution for microsatellite loci was generated, and then outlier loci were removed sequentially from the analysis, starting with those furthest from the 95% CI contour line. fdist2 was rerun iteratively with an adjusted mean F ST until all loci followed a neutral distribution. Loci that fell outside of the 95% CI of the neutral distribution were considered to be putatively non‐neutral outliers under positive selection. fdist2 was run with the following parameters: 30 samples per population, four populations, and an expected F ST of 0.083 between all populations. The model was run for 50,000 iterations, and the 95% CI of the neutral distribution was calculated. To investigate the relative strength of demographic versus selective effects on MHC allele frequency distributions, we used both neutral microsatellite loci and the adaptive MHC locus to estimate genetic relationships among these populations. fdist2 was rerun without the presumptively non‐neutral microsatellites as described above, and with an adjusted F ST of 0.045, as determined from the iterative initial runs of the model. The 95% CI of the neutral distribution was generated and plotted against the pairwise F ST values of microsatellites and the MHC DRB1 locus to determine whether the MHC DRB1 locus conformed to neutral expectations simulated in fdist2.

Because most C. gunnisoni colonies have experienced plague in recent years (including ES and SE), we tested for evidence of genetic bottlenecks in all four population samples. We used the 12 microsatellite loci in the program bottleneck (Piry et al. 1999) to test for heterozygosity excess (Cornuet and Luikart 1996) or a mode shift in the frequency of rare alleles (Luikart et al. 1998). We ran the heterozygosity excess method under two different mutation models, the stepwise mutation model (SMM) and infinite allele model (IAM), and tested significance with the Wilcoxon signed‐rank test to evaluate deviations from a 50:50 ratio of markers with heterozygosity deficiency:excess.

Survival analyses of Y. pestis challenge

We used survival data gathered from 60 C. gunnisoni individuals from a previously published Y. pestis challenge study (Busch et al. 2013) to determine whether survival in prairie dogs from ESC and AVC was associated with MHC heterozygosity and/or the presence/absence of specific DRB1 alleles. Distributions of the Kaplan–Meier survival curves for the two populations were compared using the logrank test (Harrington and Fleming 1982) in the “survival” package in R (R Core Team RCT 2013).

Results

MHC variation

The DRB1 exon 2 locus that we characterized in C. gunnisoni appears to be a single‐copy, classically expressed MHC gene with the potential to be shaped by natural selection. We detected six unique alleles (MhcCygu‐DRB1*01 to 06: GenBank accessions KU695893–KU695896 and KR338362–KR338363) among the 94 C. gunnisoni analyzed in this study (Appendix 1). BLAST comparisons against the GenBank database yielded close similarity to MHC class II DRB1 sequences from other ground squirrels, especially the spotted suslik, Spermophilus suslicus (Biedrzycka and Radwan 2008). At the C. gunnisoni exon 2 locus, we recovered a maximum of two alleles per individual and genotype frequencies were in Hardy–Weinberg equilibrium for all four population samples. This result provides compelling evidence that the C. gunnisoni exon 2 locus is a single copy. We successfully validated expression of all six alleles using spleen mRNA from a subset of 10 C. gunnisoni. Each exon 2 allele is an open‐reading frame (ORF) that encodes a class II beta subunit molecule with no insertions, deletions, or stop codons (Appendix 3). The translated amino acid alignment of this exon 2 locus (Fig. 2) has a combination of variable and conserved sites expected for a classical DRB1 locus (Kaufman et al. 1994; Bondinas et al. 2007; Yuhki et al. 2008). Together, these lines of evidence provide robust validation that this class II gene is a fully functional MHC gene in C. gunnisoni.

Alleles 05 and 06 were rare and occurred in just nine heterozygous individuals (9.6% of the 94), marked with a superscripted “m” in Appendix 1. Both the Sanger and 454 sequencing platforms missed every instance of these two alleles, and we did not discover the extent of the allelic dropout problem until all 94 individuals had been screened with our six qPCR assays. The two new reverse primers that we designed (CyguDRB1_exon2_R1 and CyguDRB1_exon2_R2) at the extreme 3′ end of exon 2 will robustly amplify all six DRB1 alleles in C. gunnisoni (Appendix 2).

Contrary to our prediction, the AV population did not exhibit greater variation at the MHC DRB1 locus when compared to ES or SE (the two plague sites). In fact, the greatest number of alleles (n = 6) was found at ES. There were only modest differences in allele frequency among these colonies (Appendix 7). Allele MhcCygu‐DRB*01 was by far the most common allele, with a frequency of 67–86%. It was also the only allele shared by all three populations. Three of the six alleles (04, 05, and 06) were fairly rare, with frequencies < 0.1. The level of H O in this MHC locus ranged from 0.27 to 0.57 (mean H O = 0.38 ± 0.1) and was virtually the same as H E (0.25–0.52) (Appendix 6). H O was not significantly different among the population samples after excluding AV2006 to avoid using redundant samples from AV (Kruskal–Wallace χ 2 = 3, df = 3, P > 0.1).

The colonies that we sampled generally exhibited a low level of nucleotide diversity at the DRB1 locus, particularly at four highly similar alleles (01, 03, 05, and 06). The uncorrected nucleotide distance of these four alleles ranged from 0.015 to 0.055, and together, they reach a very high combined frequency (80–92%) across all population samples. Most of the nucleotide diversity in the C. gunnisoni DRB1 locus is accounted for in the two remaining alleles, 02 and 04 (uncorrected distance 0.13–0.24). These two alleles appear to be trans‐species polymorphisms (Klein 1987; Takahata and Nei 1990) due to their close similarity to the DRB1 alleles of black‐tailed prairie dogs (KR338364–KR338368) and other ground squirrels (Fig. 3). The five black‐tailed prairie dogs from Montana yielded a similar pattern of nucleotide diversity to Gunnison's prairie dog, with three alleles (MhcCylu‐DRB*01, 03, and 05) showing a close relationship (uncorrected distance = 0.011–0.027) and two additional alleles (MhcCylu‐DRB*02 and 04) contributing the majority of the observed nucleotide diversity (uncorrected distance = 0.10–0.22).

Figure 3

Maximum likelihood tree of MHC‐ alleles based on exon 2 nucleotide sequence. Diamonds (♦) denote Gunnison's prairie dog, Cynomys gunnisoni, and squares (■) denote black‐tailed prairie dog, C. ludovicianus. Nodes with bootstrap confidence values ≥90% are labeled, and the scale bar units are single nucleotide polymorphisms (SNPs). GenBank accession numbers include KU695893‐KU695896 and KR338362‐KR338363 (MhcCygu‐ alleles 01–06); KR338364‐KR338368 (MhcCylu‐ alleles 01–05); JQ837902‐JQ837908 (European marmot (Marmota marmota) alleles Mama*01–07); KJ675569‐KJ675576 (Eastern woodchuck (Marmota monax) alleles Mamo‐1A*01–08) and KJ675577‐KJ575583 (alleles Mamo‐1B*01–07); EF569186‐EF569200 and HM461912‐HM461913 (spotted suslik (Spermophilus suslicus) alleles Spsu*01–18); XM_005338933 (thirteen‐lined ground squirrel (Ictidomys tridecemlineatus) allele Ictr*01); M97616 (Abert's squirrel (Sciurus aberti) allele Scab*01); NM_002124 (Human (Homo sapiens) allele HLA*01).

We found no evidence of positive selection acting on the DRB1 locus during the evolutionary history of both prairie dog species. A slightly higher rate of nonsynonymous substitutions to synonymous substitutions (d N/d S = 1.17) was observed for antigen‐binding sites (ABS) in C. ludovicianus only. Z‐tests of positive selection were nonsignificant for both species, whether considering all positions or only ABS sites (P > 0.05 for all tests). In codon‐based tests (PAML), models of neutrality (the null hypothesis) were not rejected in favor of positive selection models for C. gunnisoni. The model incorporating positive selection fit our C. ludovicianus data better than models that did not incorporate this factor based on AIC values; however, the likelihood ratio test was not significant at P = 0.05 (LR value = 5.812; Appendix 5). When comparing all five codon‐based models for C. ludovicianus, scores are similar for model 2a (the model assuming that positive selection acts on a subset of sites) and model 8 (all parameters of model 2a plus variation in ω is allowed that follows a beta‐distributed pattern of substitution rates). However, model 2a has fewer parameters and is therefore the most parsimonious model that explains our data.

Of the three commonly used measures of linkage disequilibrium employed by the program PERMUTE (r 2, D', and G4), only the r 2 estimate was consistent with recombination in C. gunnisoni (P = 0.017). All three measures were nonsignificant in C. ludovicianus.

Genetic variation at neutral versus MHC loci

Using 12 microsatellite markers, we observed moderate levels of neutral genetic diversity within the three prairie dog colonies, regardless of whether plague had affected the colonies or not (Appendix 6). H O (range 0.51–0.57) and H E (range 0.50–0.59) reached similar levels at SE, ES, and both samples of AV. All markers were in Hardy–Weinberg equilibrium in the four population samples, with the exception of marker IGS‐1 in the SE population (Appendix 6; P = 0.02). Direct comparisons of Ho (excluding the AV2006 subset) were not significantly different (Wilcoxon rank sum W = 88.1; P > 0.1). The total number of alleles per locus was low (range 2–6), and no difference was observed in the mean number of alleles per population sample (3.4–4.1) (mean Wilcoxon rank sum W = 96.2; P > 0.1; Appendix 6). These low‐to‐moderate diversity measures may indicate that all populations have experienced fairly similar demographic histories, regardless of whether plague was present or absent in recent decades.

The level of population structure (F ST) estimated by the microsatellite markers (Table 1) was low‐to‐moderate and did not seem to associate with the presence of plague. The mean F ST (excluding the AV2006 subset) was 0.098 (95% CI: 0.035–0.149). No evidence of differentiation over time (F ST = 0) was observed between the two AV population samples from 2006 and 2009 (AV2006 vs. AVC in Table 1). Genetic differentiation was generally highest in pairwise comparisons that involved ESC (F ST = 0.09–0.12), which is expected given its greater geographic distance (~66 km) from the SE and AV colonies. The SE colony is <6 km distant from AV at its closest border, but was significantly differentiated from both the AV2006 and AVC samples (F ST = 0.05 and 0.08).

Table 1

Population structure estimates (F ST) for all pairwise population comparisons using microsatellite loci (top panel below diagonal) and the MHC DRB1 locus (bottom panel below diagonal). Asterisks denote P‐value significance levels (*0.05, **0.005, ***0.0005, ns is not significant). The top panel also includes 95% confidence intervals for microsatellite F ST estimates above the diagonal. Population samples include Aubrey Valley (AV 2006 and AVC); Seligman 2006 (SE); and Espee Ranch (ESC). Global estimates do not include the AV 2006 population sample because this would cause a downward bias from using two AV samples. The “C” subscript denotes AVC and ESC samples collected in 2009 that were used in a previous laboratory plague challenge experiment (Busch et al. 2013)

	AV 2006	SE 2006	ESC	AVC
MSAT F ST comparisons (with 95% CI above diagonal)
AV 2006	–	0.003–0.122	0.020–0.183	−0.012–0.022
SE 2006	0.054***	–	0.036–0.148	0.031–0.146
ESC	0.093***	0.091***	–	0.039–0.226
AVC	0.002 (ns)	0.080***	0.123***	–
Global	0.098 (0.035–0.149)
DRB1 F ST comparisons (below diagonal)
AV 2006	–
SE 2006	0.071**	–
ESC	−0.005 (ns)	0.051**	–
AVC	−0.004 (ns)	0.03*	0.018*	–
Global	0.033

The MHC locus generally displayed lower levels of genetic differentiation (global F ST = 0.033) compared to neutral microsatellite markers (Table 1). The colonies known to suffer from plague outbreaks (SE and ESC) displayed about half the level of differentiation at the DRB1 locus (F ST = 0.051, P < 0.05) compared to microsatellites (F ST = 0.091, P < 0.0002). Interestingly, colonies AVC and ESC displayed very low differentiation at DRB1 (F ST = 0.018, P < 0.05), which was significantly less than the 95% CI of the microsatellite loci for these two colonies (F ST = 0.12; CI = 0.03–0.22). The relatively narrow range of F ST variation at DRB1 compared to microsatellite loci could be an indication of past selective pressure acting on the MHC. However, these results do not suggest a simple correlation with plague presence/absence in these prairie dog colonies.

We also used a simulation technique in fdist2 to compare observed versus expected F ST values from neutral microsatellites against the DRB1 locus (Fig. 4). As the DRB1 locus is subject to selection pressures (e.g., plague or other disease), we would expect to see F ST values that are closer to zero than the neutral microsatellite loci. Although the observed F ST estimates for the MhcCygu‐DRB1 locus were low, they all fell within the lower 95% CI of the simulated fdist2 distributions estimated from microsatellite markers (Fig. 4). The relatively narrow range of F ST values at DRB1 compared to the FDIST2 simulation may suggest that past selection has acted on the MHC, but the selective signature seems to be partially masked by demography. The outlier analysis identified three microsatellite markers that are potentially non‐neutral loci: IGS‐BP1, IGS‐6, and A111 (labeled individually in Fig. 4). We therefore reran the FSTAT analysis after removing these three loci and still recovered a global F ST estimate (F ST = 0.065; CI = 0.04–0.09) that was significantly greater than the DRB1 locus (0.033). It is important to note that the fdist2 result shown in Figure 4 is conservative because the simulation yielded broader 95% CIs than the FSTAT analysis.

Figure 4

Plot of expected heterozygosity () versus that tests whether the Cynomys gunnisoni MHC locus conforms to the level of genetic structure expected for neutral loci. In this analysis, the observed values from the MHC locus are directly compared to expected values from microsatellites generated with four simulated populations in fdist2 (middle line). The top and bottom lines define the 95% confidence intervals (CIs) of the simulated distribution of relative to expected heterozygosity () for 12 presumably neutral microsatellite loci. values outside of the 95% CIs for neutral loci are suggestive of a locus under selection pressure. Three potentially non‐neutral loci are individually labeled (IGS‐BP1, IGS‐6, and A111). Only three of the four Cynomys gunnisoni populations were included in the MHC pairwise comparisons because AV was sampled twice (AV2006 and AV2009). We therefore excluded AV2006 to avoid introducing bias into this analysis. Triangles represent pairwise MHC values; the circle represents the global MHC value; crosses represent pairwise microsatellite values.

In bottleneck tests, only the ES colony showed strong evidence of a recent population reduction (Appendix 8). Simulations under the IAM led to expected heterozygosity (H E) excess in 11 of 12 markers at ES (P = 0.002). The SMM model yielded nine markers with H E excess at ES (P = 0.09). The only other colony with a bottleneck signature was SE, with 11 of 12 markers showing H E excess under the IAM (P = 0.034) but not under the SMM (P = 0.733). It is worth noting that the IAM is a less conservative model and has been known to identify heterozygosity excess in nonbottlenecked populations (Luikart and Cornuet 1998). The ratio of markers with H E deficiency versus excess was much more balanced at AV in both sample years, which is consistent with the demographic stability observed at AV in recent decades. In all colonies, the frequency of rare alleles matched that of a normal L‐shaped distribution. Hence, no bottleneck signal was observed using this method.

Survival analyses

Comparison of Kaplan–Meier survival curves (Fig. 5) among C. gunnisoni subjected to Y. pestis challenge showed that survival rates were significantly higher in individuals that possessed at least one copy of MhcCygu‐DRB1*01 compared to individuals that did not possess this allele (60% vs. 20%) (χ 2 = 5.8, df = 1, P = 0.02). An important caveat to this result is that the sample sizes were very unbalanced: 55 prairie dogs from ESC and AVC had at least one copy of allele 01, but only five individuals did not carry it. MHC heterozygosity at the DRB1 locus was not a significant factor in determining survival during plague challenge (χ 2 = 0.1, df = 1, P = 0.79), and of the 34 prairie dogs that were homozygous for allele 01, 19 survived and 15 died (Appendix 1). In contrast, every prairie dog carrying a copy of allele 04 (n = 3) died when challenged with plague, including a 0104 heterozygote.

Figure 5

Kaplan–Meier survival curves showing that Cynomys gunnisoni individuals possessing at least one MhcCygu‐*01 allele survived significantly longer after challenge with a single injection of Yersinia pestis strain CO92 compared to individuals not possessing this allele (χ 2 = 5.8, df = 1, P = 0.02). Plague challenge survivorship data were collected in a previous study (Busch et al. 2013).

Discussion

Yersinia pestis is endemic in western North America and a significant pathogen of native rodents and associated wildlife, such as black‐footed ferrets. The four species of prairie dogs in the United States have experienced recurring selection pressure from plague since the 1930s, with mortality rates approaching 100% at affected colonies (Cully et al. 1997; Cully and Williams 2001; Hoogland et al. 2004). It is possible that survival to plague is due to a genetic mechanism (as opposed to an antibody‐based mechanism) and, if so, the colonies that repopulate from a handful of resistant individuals may show improved survival during subsequent plague outbreaks. Two recent experimental challenge studies in C. gunnisoni and C. ludovicianus revealed that 50–70% of prairie dogs from plague‐endemic areas can survive a single injection of the highly virulent Y. pestis strain CO92 (Rocke et al. 2012; Busch et al. 2013). Survival was not correlated with the initial delivery dosage of Y. pestis in either prairie dog species because some animals succumbed to the lowest dose (50 bacteria) while others survived the highest dose (50,000 bacteria). This type of variation is necessary for developing an evolutionary response to a highly virulent pathogen such as Y. pestis. Here, we discuss the association of MHC alleles in prairie dogs experimentally challenged with plague and present the first description of MHC diversity in C. gunnisoni.

Variation at the DRB1 locus in C. gunnisoni was lower than in most other rodent species, but similar to that observed in ground squirrels. With just six alleles, C. gunnisoni allelic richness is extremely low compared to mice and rats (de Bellocq et al. 2008) but about the same as seen in recently fragmented populations of alpine marmot (Kuduk et al. 2012) and spotted suslik (Biedrzycka and Radwan 2008). The low number of alleles could possibly be due to recent bottlenecks in C. gunnisoni, which like all prairie dog species, have experienced declines from plague, poisoning and habitat loss. Indeed, we found evidence of recent bottlenecks at the two plague sites (SE and ES), although not at AV. Fortunately, the nucleotide diversity among existing alleles in C. gunnisoni (0.015–0.24) and C. ludovicianus (0.011–0.22) remains high and is similar to that reported for other ground squirrel species (Biedrzycka and Radwan 2008; Kuduk et al. 2012; Moreno‐Cugnon et al. 2015). The co‐occurrence of highly divergent DRB1 alleles provides evidence that balancing selection has occurred in the course of prairie dog evolution, as expected for this class II MHC locus. In contrast to the high level of nucleotide diversity at DRB1, levels of H E were low in C. gunnisoni (0.27–0.53) compared to spotted suslik (0.63–0.89), eastern woodchuck (0.77–0.83), and alpine marmot (0.51–0.63) (Biedrzycka and Radwan 2008; Kuduk et al. 2012; Moreno‐Cugnon et al. 2015). DRB1 heterozygosity in C. gunnisoni was also lower than H E at microsatellite loci (0.50–0.59). These microsatellite H E estimates fell in the middle range of C. ludovicianus (0.53–0.63) (Roach et al. 2001; Sackett et al. 2012, 2013; Castellanos‐Morales et al. 2015) and other sciurids (0.33–0.75) (Stevens et al. 1997; Reid et al. 2010; Ricanova et al. 2011; Fitak et al. 2013) and suggests that C. gunnisoni is maintaining about the same level of genetic diversity as these phylogenetic relatives. The level of genetic structure among the three C. gunnisoni populations we studied was somewhat less than regional differentiation in C. gunnisoni zuniensis (0.38) (Sackett et al. 2014), C. ludovicianus (0.54) (Sackett et al. 2012, 2013; Castellanos‐Morales et al. 2015), and C. mexicanus (0.17) (Castellanos‐Morales et al. 2015). As such, the level of genetic differentiation and diversity at microsatellites in our study populations seems typical for prairie dogs and other ground squirrel species.

Selection appears to have played an important role in shaping patterns of genetic structure at the MHC DRB1 locus in C. gunnisoni. This has led to a reduced level of genetic differentiation (F ST) at the MHC locus compared to the level at neutral markers (Table 1). The variation in pairwise DRB1 F ST estimates (0.018–0.051) is also much narrower than that for 12 microsatellites (Fig. 4). However, two lines of evidence suggest that demography has also imparted a detectable signal on the DRB1 locus. First, all of the pairwise F ST values for DRB1 lie inside the 95% confidence interval generated by the fdist2 simulation based on microsatellite markers (Fig. 3). The same outcome has been observed in another ground squirrel, the spotted suslik (Spermophilus suslicus) (Biedrzycka and Radwan 2008). In both species, genetic drift appears to be more important than balancing selection in shaping MHC variation. This pattern has also been reported in other strongly bottlenecked or fragmented populations (Miller and Lambert 2004; Strand et al. 2012; Sutton et al. 2015). Second, we did not find evidence of positive selection acting on specific amino acid positions in C. gunnisoni, even though we did find this in C. ludovicianus. Genetic drift may explain the absence of this expected selective signature. Additional explanations may be that epizootic plague is not the only source of selection in these colonies and plague does not affect every prairie dog generation. Assuming a generation time of 1–2 years, C. gunnisoni has existed with endemic plague for a maximum of 85 generations and perhaps not enough time has elapsed to detect an evolutionary response at MHC loci. Another mechanism that can explain variation in F ST values at the DRB1 locus is that selection pressures on the MHC can fluctuate over time and space (Hill 1991; Hedrick 2002). In this case, determining the relative influence of plague may require sampling at a broader geographic scale. It is known that the average F ST among 34 rangewide samples of C. gunnisoni subspecies zuniensis is 0.38 for microsatellite markers (Sackett et al. 2014), which is considerably greater than the level we observed (microsatellite F ST 0.08–0.123 for AV‐SE‐ES). It is possible that MHC differentiation across these same 34 rangewide locations will be much lower.

In genotyping the 60 C. gunnisoni that were used in the previous experimental plague challenge (Busch et al. 2013), we found evidence for a potential MHC‐based survival advantage in individuals that carried at least one copy of allele MhcCygu‐DRB1*01. Despite a large number of challenge studies on laboratory rodents, only a handful of plague challenge studies have been conducted using wild rodents besides prairie dogs (Quan and Kartman 1962; Rust et al. 1972; Quan et al. 1985; Thomas et al. 1988, 1989; Tollenaere et al. 2010b). Although the presence of the 01 allele does not guarantee survival during a plague infection, the absence of this allele appears to be disadvantageous (60% vs. 20% survival). Individuals carrying allele 04 seemed particularly ill‐suited for fighting plague, as the three prairie dogs with this allele all succumbed to the infection, including a 0104 heterozygote. We recognize the potential problem of imbalanced sample size for these allelic comparisons, as 55 prairie dogs had at least one copy of the 01 allele but only five individuals lacked this allele. The pattern that we observed might suggest the common MHC allele is an important factor for prairie dog survival during plague infection, but it is likely not the only factor. It may also be possible that plague survival is influenced by a closely linked gene rather than this MHC class II locus. The six qPCR assays that we developed will facilitate rapid MHC genotyping in future studies of prairie dogs and plague.

We found evidence of trans‐species polymorphisms for two allelic lineages (02 and 04) that shared high sequence similarity between C. gunnisoni and C. ludovicianus (Fig. 3). These two species do not hybridize and are members of different subgenera that last shared a common ancestor about five million years ago (Goodwin 1995; Herron et al. 2004). Alleles 02 and 04 have been maintained over an evolutionary time frame by C. gunnisoni and C. ludovicianus with only a small number of species‐specific changes (Appendix 3). Therefore, we predict that very similar alleles will also be found in C. parvidens and C. leucurus. The maintenance of ancient polymorphisms through speciation is frequently observed in MHC loci (Klein 1987; Figueroa et al. 1988; Sena et al. 2003; Califf et al. 2013) and may be due to balancing selection over extended evolutionary timescales (Takahata and Nei 1990) or convergent evolution from similar pathogen pressures (Edwards and Hedrick 1998; Yeager and Hughes 1999). Some of the prairie dog alleles share past ancestry with major allele lineages found in the spotted suslik (Spermophilus suslicus), an Asian ground squirrel that naturally occurs in areas with native plague. Ancient polymorphisms shared by ground squirrels species would probably have arisen before the divergence of Y. pestis from its ancestor (Y. pseudotuberculosis) in the past 6000 years (Cui et al. 2013; Rasmussen et al. 2015). In a broad sense, this means that these MHC alleles could still be experiencing selection by this recently evolved pathogen. It remains to be seen whether all prairie dog species share MHC allele lineages that may provide a survival advantage during Y. pestis infection.

Most of what is currently known about the genetic basis of plague resistance is based on evidence obtained from laboratory experiments on model organisms. The existing studies suggest that plague resistance is a polygenic trait, probably because Y. pestis attacks the immune system with a complex array of virulence mechanisms. Our study is generally consistent with this idea, because the class II MHC locus is clearly not the only factor determining, which prairie dogs will survive plague. In the prairie dog challenge experiments, very few individuals of either species produced a strong antibody response, indicating that other immune mechanisms (e.g., early inflammation or innate pathways) may be more critical for surviving plague infection (Rocke et al. 2012; Busch et al. 2013). In laboratory mice, at least four plague resistance loci (prls) contribute to high resistance levels in Mus spretus (Blanchet et al. 2011; Chevallier et al. 2013) and another prl maps to the MHC region of M. musculus line BALB/cJ (Turner et al. 2008). Resistance in M. musculus line SEG is associated with rapid induction of early innate immune responses that occur in the first 2–3 days of infection (Demeure et al. 2011). The MHC class II apparatus is known to present the epitopes from several Y. pestis proteins (Parent et al. 2005; Shim et al. 2006), including the F1 and LcrV antigens that are used as a plague vaccine for prairie dogs (Mencher et al. 2004; Rocke et al. 2008, 2010). However, Y. pestis can use an effector protein (Yop M) to suppress MHC class II antigen presentation (Soundararajan et al. 2011; Bosio et al. 2012) that likely serves as a mechanism to interrupt the activation of T cells. The results of our challenge study raise the possibility that some class II alleles might be more easily suppressed with this type of mechanism than others. It could also be that certain prairie dog alleles are simply less effective at binding and presenting Y. pestis antigens. The application of genomewide scans to find new plague resistance markers (Tollenaere et al. 2010a, 2013), including those within the MHC (Tollenaere et al. 2012), will hopefully lead to a better understanding of the complex immune responses that mammals employ during exposure to plague.

Conclusions

Characterizing genes integral to the immune system is an important first step in identifying possible mechanisms of resistance to diseases in nonmodel species like prairie dogs. This is the first published report of an expressed MHC class II DRB1 locus in any prairie dog species. The small number of alleles at this DRB1 locus all share similarity to other ground squirrels, including C. ludovicianus and the spotted suslik (S. suslicus). We found weak evidence for positive selection acting on the DRB1 locus in the black‐tailed prairie dog, but not in Gunnison's prairie dog, which may be a result of recent population reductions or the short time frame during which plague has been acting as a selective force. Plague is probably not the only source of selective pressure in this system, and other diseases (e.g., Francisella tularensis) or parasites may play an important role in shaping variation at immune system loci in prairie dogs (Zeidner et al. 2004). In addition to disease, two other important factors have contributed to major declines in all prairie dog species in the last 100 years: habitat conversion and poisoning campaigns (Roemer and Forrest 1996; Wagner et al. 2006). In combination, these factors may have contributed to reductions in the genomic diversity found in C. gunnisoni at both neutral loci and genes under selection. We are currently investigating whether greater microsatellite and MHC diversity existed prior to the introduction of plague to the western United States by genotyping C. gunnisoni specimens collected prior to 1930. These studies will contribute to an understanding of the development of plague resistance in prairie dogs and other rodents, and help elucidate the influence of diseases in shaping natural genetic variation in wild populations.

Conflict of Interest

The authors declare no conflict of interest.

Data Archiving

Sequence data have been submitted to GenBank under accession numbers GenBank accessions KU695893–KU695896 and KR338362–KR338368.

Sample	Pop	Collection date	Zone 12N UTM_N	Zone 12N UTM_E	Sex	Age	Capture wt. (g)	Inf. dose (cfu)	Fate	DPI	DRB1 allele 1	DRB1 allele 2
AV01 V	AV	24AUG2005	3916926	320515	M	A	–	–	–	–	02	03
AV02	AV	24AUG2005	3917126	320555	F	A	–	–	–	–	02	06 m
AV03 V	AV	24AUG2005	3917823	320968	M	A	–	–	–	–	01	02
AV04	AV	24AUG2005	3918959	320552	M	A	–	–	–	–	01	01
AV05	AV	24AUG2005	3918994	320603	M	A	–	–	–	–	01	01
AV06	AV	24AUG2005	3919256	320725	F	A	–	–	–	–	01	01
AV09	AV	26AUG2006	3921136	314603	U	U	–	–	–	–	01	02
AV10	AV	26AUG2006	3922121	313441	F	U	–	–	–	–	01	01
AV11	AV	26AUG2006	3922420	313017	U	U	–	–	–	–	01	01
AV12	AV	26AUG2006	3922520	313161	F	U	–	–	–	–	01	03
AV13	AV	26AUG2006	3922748	312612	M	U	–	–	–	–	01	03
AV14	AV	26AUG2006	3927834	310982	M	U	–	–	–	–	01	01
AV15	AV	26AUG2006	3926003	308507	M	U	–	–	–	–	01	02
AV16	AV	26AUG2006	3922808	312411	M	U	–	–	–	–	01	03
AV17	AV	26AUG2006	3922615	312647	F	U	–	–	–	–	01	01
AV18	AV	26AUG2006	3922593	312522	U	U	–	–	–	–	01	01
SE07	SE	09SEP2006	3909725	327090	U	U	–	–	–	–	01	01
SE08	SE	09SEP2006	3908249	325340	F	U	–	–	–	–	01	06 m
SE19	SE	09SEP2006	3909656	327066	M	A	–	–	–	–	01	01
SE20	SE	09SEP2006	3907874	325037	M	U	–	–	–	–	01	01
SE21	SE	09SEP2006	3907805	325694	F	U	–	–	–	–	01	01
SE22	SE	09SEP2006	3907669	325116	M	U	–	–	–	–	01	01
SE23	SE	09SEP2006	3907338	325471	M	U	–	–	–	–	01	01
SE24	SE	09SEP2006	3906958	325301	M	U	–	–	–	–	01	01
SE25 V	SE	09SEP2006	3907279	325606	M	U	–	–	–	–	01	04
SE26	SE	09SEP2006	3907775	325014	M	U	–	–	–	–	01	06 m
SE27	SE	09SEP2006	3911917	327397	M	U	–	–	–	–	01	01
SE28	SE	09SEP2006	3907673	330195	M	U	–	–	–	–	01	01
SE29	SE	09SEP2006	3913203	332392	M	U	–	–	–	–	01	01
SE30	SE	09SEP2006	3914578	332705	M	U	–	–	–	–	01	01
SE31 V	SE	09SEP2006	3913326	327153	M	U	–	–	–	–	01	04
SE32	SE	09SEP2006	3913218	327220	F	U	–	–	–	–	01	01
SE33	SE	09SEP2006	3912854	327322	F	U	–	–	–	–	01	01
SE34	SE	09SEP2006	3911973	327342	F	U	–	–	–	–	01	04
60‐28	ESC	01SEP2009	3961380	356869	F	J	500	50000	SURV	30	01	01
61‐30	ESC	01SEP2009	3961020	357068	F	A	850	50000	DIED	9	01	01
63‐29 V	ESC	01SEP2009	3961430	356880	M	J	750	50	DIED	12	02	02
64‐19	ESC	31AUG2009	3967379	360787	M	J	700	50000	SURV	30	01	03
65‐20	ESC	01SEP2009	3967386	360711	M	J	550	50	DIED	8	01	02
66‐23	ESC	01SEP2009	3967438	360489	M	A	1100	50000	SURV	30	01	01
67‐24	ESC	01SEP2009	3967226	360624	M	A	1000	50	SURV	30	01	01
68‐33	ESC	01SEP2009	3961114	356884	M	A	1175	5000	DIED	10	01	01
69‐04	ESC	25AUG2009	3961396	356886	F	J	525	5000	SURV	30	01	01
70‐22	ESC	01SEP2009	3967386	360723	F	J	475	50	SURV	30	01	05 m
71‐17 V	ESC	29AUG2009	3961555	357242	M	J	700	50000	DIED	11	02	04
72‐12	ESC	29AUG2009	3967200	360621	F	J	575	5000	DIED	16	01	01
73‐31 V	ESC	01SEP2009	3961575	357274	F	A	675	5000	DIED	8	04	06 m
74‐35	ESC	01SEP2009	3967193	360751	M	J	500	5000	DIED	4	01	03
75‐32	ESC	01SEP2009	3961076	357066	F	J	625	5000	SURV	30	01	01
76‐21	ESC	01SEP2009	3967365	360690	F	J	500	50000	DIED	16	01	01
77‐01	ESC	29AUG2009	3967194	360872	M	J	750	5000	DIED	13	01	04
78‐08	ESC	25AUG2009	3967393	360730	M	A	850	50000	DIED	16	01	01
79‐16	ESC	29AUG2009	3967423	360681	M	J	650	50000	DIED	10	01	01
80‐34	ESC	01SEP2009	3961086	356990	F	A	675	50	SURV	30	01	06 m
81‐05	ESC	25AUG2009	3961380	356869	F	J	525	50	SURV	30	01	05 m
82‐03	ESC	25AUG2009	3961387	356881	M	A	925	50	DIED	12	01	01
83‐36	ESC	04SEP2009	3967525	360500	M	J	600	50	SURV	30	01	05 m
85‐18	ESC	29AUG2009	3961083	357887	F	J	600	50	SURV	30	01	06 m
86‐06	ESC	25AUG2009	3961362	356862	M	J	600	5000	SURV	30	01	02
87‐25	ESC	01SEP2009	3967535	360540	M	J	700	50000	SURV	30	01	01
88‐81	ESC	29AUG2009	3967373	360804	M	J	450	50	DIED	13	01	02
89‐02	ESC	25AUG2009	3961426	356866	F	A	700	5000	SURV	30	01	02
90‐07	ESC	25AUG2009	3961349	356873	F	A	750	50000	DIED	11	01	02
127‐14	ESC	29AUG2009	3967395	360730	M	J	575	5000	SURV	30	01	02
92‐49	AVC	15SEP2009	3928536	309811	M	J	650	50000	SURV	30	01	01
93‐56	AVC	15SEP2009	3922057	312544	F	A	700	50	SURV	30	01	02
94‐40	AVC	15SEP2009	3927905	308962	M	J	675	50	SURV	30	01	01
95‐65	AVC	18SEP2009	3921130	313790	M	J	625	5000	SURV	30	01	01
96‐53	AVC	15SEP2009	3922191	312523	M	A	1000	50000	SURV	30	01	01
97‐57 V	AVC	15SEP2009	3921451	313231	M	J	775	5000	DIED	6	02	02
98‐63	AVC	18SEP2009	3921070	313751	M	J	750	50000	SURV	30	01	03
99‐69	AVC	18SEP2009	3921837	313308	M	A	800	5000	DIED	11	01	01
100‐46	AVC	15SEP2009	3927956	309420	F	A	775	5000	SURV	30	01	01
101‐54	AVC	15SEP2009	3922177	312600	M	J	725	5000	SURV	30	01	01
102‐62	AVC	18SEP2009	3927972	309418	M	A	1150	50	SURV	30	01	01
103‐11	AVC	28AUG2009	3922013	313157	M	J	750	50000	DIED	7	01	01
104‐43	AVC	15SEP2009	3928034	309467	F	J	525	5000	SURV	30	01	01
105‐55 V	AVC	15SEP2009	3922100	312478	M	J	750	50	SURV	30	02	03
106‐09	AVC	28AUG2009	3922043	313167	F	A	700	50000	DIED	9	01	01
107‐67	AVC	18SEP2009	3921785	313179	F	J	550	50000	SURV	30	01	02
108‐66	AVC	18SEP2009	3921236	313762	F	A	775	5000	DIED	11	01	01
110‐48	AVC	15SEP2009	3928476	309863	F	J	550	5000	DIED	8	01	01
112‐45	AVC	15SEP2009	3928004	309439	M	J	600	50000	SURV	30	01	01
114‐51	AVC	15SEP2009	3922434	312012	M	A	875	50	SURV	30	01	01
115‐44	AVC	15SEP2009	3927901	309393	F	J	500	50	SURV	30	01	01
116‐61	AVC	18SEP2009	3928026	309279	F	A	650	50000	DIED	5	01	01
117‐41 V	AVC	15SEP2009	3927918	308938	M	J	725	5000	SURV	30	01	03
118‐68	AVC	18SEP2009	3921579	313264	M	A	850	50000	SURV	30	01	02
119‐47	AVC	15SEP2009	3928460	309703	M	J	750	50	SURV	30	01	02
120‐58	AVC	15SEP2009	3921514	313260	M	J	750	50	SURV	30	01	02
121‐10	AVC	28AUG2009	3922876	312444	F	J	625	50	DIED	7	01	01
122‐50	AVC	15SEP2009	3922376	311977	F	J	600	50	SURV	30	01	01
123‐42	AVC	15SEP2009	3928066	309474	F	J	625	50000	SURV	30	01	01
124‐52	AVC	15SEP2009	3922546	312087	M	A	1050	5000	DIED	17	01	01

m“Missing allele” that was not observed in sequences generated from the Sanger and 454 platforms.

V“Validation individual” used to confirm expression of the DRB1 gene using mRNA from spleen tissue.

Name	Sequence (5′–3′)	Target allele	T a (°C)	Size (bp)	Melt temp (°C)	Comments
PCR and sequencing primers
GH46F	CCGGATCCTTCGTGTCCCCACAGCACG	All	65	271	–	Erlich and Bugawan 1990
GH50R	CTCCCCAACCCCGTAGTTGTGTCTGCA	All	65	271	–	Erlich and Bugawan 1990
CyguDRB1_exon2_R1	CTCTCCTCTCCACAGTGAAGCTCTCA	All (exon 2) except Cygu‐DRB1*03	65	293		New reverse primer specific to C. gunnisoni.
CyguDRB1_exon2_R2	CTCTCCGCTCCACAGCGAAGCTCTCA	Cygu‐DRB1*03	65	293	–	New reverse primer specific to C. gunnisoni.
CyguDRB1exon3‐R2	AGACCAGGAGGTTGTGGTGCTG	All (exon 3)	65	334	–	For cDNA PCR and sequencing.
454 adapter‐A	cgtatcgcctccctcgcgccatcagNNNNNNNNNN	–	64	306	–	Tail for GH46F (with 10 bp indices).
454 adapter‐B	ctatgcgccttgccagcccgctcag	–	64	296	–	Tail for reverse primer (e.g., GH50R or CyguDRB1_exon2_R1 and R2).
Assays 1–4: Noncompetitive PCRs for all known CyguDRB1 alleles
PD1DRB_F	CGTTTCCTGGAGCAAGTTTCACATG	All	67	96	–	This forward primer can amplify all six alleles.
PD1DRB_R	CGAACTCCTCTCGGTTGTGG	Cygu‐DRB1*01 Cygu‐DRB1*05 Cygu‐DRB1*06	67	96	80–81	This reverse primer amplifies three highly similar alleles (01, 05, and 06).
PD2DRB_F2	GGAGCGGATACGGTTCCTAC	Cygu‐DRB1*02	67	48	73.5–74.5
PD2DRB_R	GACCTCCTCCCGGTTATAGAT	Cygu‐DRB1*02	67	48	73.5–74.5
PD3DRB_F2	TTTCCTGGAGCAAGGTTCACA	Cygu‐DRB1*03	67	102	80.75–81.75
PD3DRB_R2	AAAGCGCACGTACTCCTCCA	Cygu‐DRB1*03	67	102	80.75–81.75
PD4DRB_F2	GTTTCCTGGAGCAAGGTTTCAG	Cygu‐DRB1*04	67	183	86–87
PD4DRB_R2	GCGTCCTTCTGGCTGTTGA	Cygu‐DRB1*04	67	183	86–87
Assay 5: Noncompetitive PCRs to resolve allele 01
PD5DRB_F	CGTTTCCTGGAGCAAGTTTCACATG	All	65		–	Same sequence as primer PD1DRB_F above.
PD5DRB_R2_Anc	CCTCTCCACAGTGAAGCTCTCACCA	Cygu‐DRB1*01 (poor amp for 05, 06)	67	264	87–88	PCR with this reverse primer will amplify three highly similar alleles (01, 05, and 06). However, the C t value for allele 01 will occur 8–10 cycles earlier compared to a PCR using the alternate reverse primer (PD5DRB_R2_Der).
PD5DRB_R2_Der	TCTCCTCTCCACAGTGAAGCTCTCAAAC	Cygu‐DRB1*05 Cygu‐DRB1*06 (poor amp for 01)	67	267	87–88	PCR with this reverse primer will amplify three highly similar alleles (01, 05, and 06). However, the C t value for alleles 05 or 06 will occur 8–10 cycles earlier compared to a PCR using the alternate reverse primer (PD5DRB_R2_Anc).
Assay 6: Competitive PCR to resolve allele 05 versus 06
PD6DRB_164G_F_Anc	AGTACCGCGCGGTGAGAG	Cygu‐DRB1*01 Cygu‐DRB1*05	56	33	75.5–76.5 and 84–85	Double peaks observed in dissociation curve.
PD6DRB_164A_F_Der	gcccgcccgcccgcccAGTACCGCGCGGTGAGTA	Cygu‐DRB1*06	56	49	82–83
PD6DRB_R_Cons	CCGGCCGCCCCAGCT	All alleles (two SNPs in D allele)	56		–	Conserved reverse primer.

	Outcome of 6 qPCR assays						Comments
DRB1 alleles	1	2	3	4	5	6
By allele
Cygu‐DRB1*01	1	0	0	0	Anc	Anc	Most common allele in western AZ populations
Cygu‐DRB1*02	0	1	0	0	n/a	n/a
Cygu‐DRB1*03	0	0	1	0	n/a	n/a
Cygu‐DRB1*04	0	0	0	1	n/a	n/a
Cygu‐DRB1*05	1	0	0	0	Der	Anc	Very similar to allele Cygu‐DRB1*01
Cygu‐DRB1*06	1	0	0	0	Der	Der	Very similar to allele Cygu‐DRB1*01
By diploid genotype
0101	1	0	0	0	Anc	Anc	Most common diploid genotype
0102	1	1	0	0	Anc	Anc
0103	1	0	1	0	Anc	Anc
0104	1	0	0	1	Anc	Anc
0105	1	0	0	0	Anc + Der	Anc
0106	1	0	0	0	Anc + Der	Der
0202	0	1	0	0	n/a	n/a
0203	0	1	1	0	n/a	n/a
0204	0	1	0	1	n/a	n/a
0205	1	1	0	0	Der	Anc	Not observed
0206	1	1	0	0	Der	Der
0303	0	0	1	0	n/a	n/a	Not observed
0304	0	0	1	1	n/a	n/a	Not observed
0305	1	0	1	0	Der	Anc	Not observed
0306	1	0	1	0	Der	Der	Not observed
0404	0	0	0	1	n/a	n/a	Not observed
0405	1	0	0	1	Der	Anc	Not observed
0406	1	0	0	1	Der	Der
0505	1	0	0	0	Der	Anc	Not observed
0506	1	0	0	0	Der	Anc + Der	Not observed
0606	1	0	0	0	Der	Der	Not observed

Cynomys gunnisoni
Model†	ln L §	ΔAIC	Parametersφ
M0 – one ω	−675.361	37.50	0.290
M1a – nearly neutral (ω 0  < 1, ω 1  = 1)	−655.611	best *	p 0  = 0.733, ω 0  = 0.034
M2a – positive selection (ω 0 < 1, ω 1 = 1, ω 2 > 1)	−653.777	0.331	p0 = 0.805, p2 = 0.195,ω 0 = 0.073, ω 2 = 2.145
M7 – beta (p, q)	−656.058	0.892	p = 0.017; q = 0.034
M8 – beta and ω (p, q, ωs > 1)	−653.782	0.341	p0 = 0.806, (p1 = 0.194), p = 8.015, q = 99.0, ω = 2.155
Cynomys ludovicianus
Model†	ln L §	ΔAIC	Parametersφ
M0 – one ω	−552.710	26.682	0.325
M1a – nearly neutral (ω 0 < 1, ω 1 = 1)	−542.278	5.818	p0 = 0.654, ω 0 = 0.059
M2a – positive selection (ω 0  < 1, ω 1  = 1, ω 2  > 1)	−539.369	best *	p 0  = 0.631, p 2  = 0.091, ω 0  = 0.089, ω 2  = 13.813
M7 – beta (p, q)	−542.612	6.487	p = 0.111; q = 0.171
M8 – beta and ω (p, q, ωs > 1)	−539.275	−0.187	p0 = 0.905, (p1 = 0.095), p = 0.253, q = 0.502, ω = 12.803

Estimated proportions of sites (p x) evolving at corresponding estimated rates (ω x = d N/d S) are given in the parameters column.

†Alternative models of codon evolution from PAML (Yang 1997, 2007).

§Log likelihood score.

φProportion of sites (p) evolving at corresponding rate (ω).

*Model of best fit.

		# Alleles (w/richness)				H O (H E)				F IS
Locus	T a (°C)	AV	SE	ESC	AVC	AV	SE	ESC	AVC	AV	SE	ESC	AVC
aA2	55	4 (4)	3 (2.9)	5 (4.8)	4 (3.5)	0.56 (0.59)	0.61 (0.52)	0.63 (0.66)	0.47 (0.42)	−0.118	−0.133	−0.033	−0.082
aA104	51	6 (6)	6 (5.8)	6 (5.7)	6 (5.7)	0.75 (0.73)	0.77 (0.69)	0.83 (0.74)	0.90 (0.74)	0.016	−0.087	−0.093	−0.192
aA111	55	2 (2)	2 (2)	3 (2.9)	3 (2.8)	0.68 (0.48)	0.39 (0.46)	0.43 (0.54)	0.47 (0.52)	0.266	0.031	−0.017	0.071
aD2	51	5 (5)	5 (4.8)	5 (4.4)	4 (3.9)	0.50 (0.65)	0.66 (0.66)	0.73 (0.72)	0.57 (0.59)	0.167	−0.041	0.151	−0.322
aD12	59	4 (4)	3 (2.8)	4 (3.4)	4 (3.5)	0.43 (0.41)	0.11 (0.11)	0.30 (0.31)	0.30 (0.26)	0.211	0.227	0.151	−0.074
aD115	59	6 (6)	4 (3.9)	6 (5.9)	6 (5.8)	0.75 (0.77)	0.61 (0.68)	0.73 (0.74)	0.80 (0.78)	−0.034	−0.199	0.047	0.360
bIGS‐1	55	5 (5)	H5 (4.9)	5 (4.7)	5 (4.9)	0.75 (0.74)	0.44 (0.66)	0.57 (0.64)	0.80 (0.76)	0.060	0.130	0.033	−0.004
bIGS‐6	59	2 (2)	2 (2)	3 (3)	2 (1.9)	0.25 (0.31)	0.22 (0.27)	0.57 (0.65)	0.17 (0.15)	−0.040	−0.015	0.070	−0.101
bIGS‐9d	51	4 (4)	3 (3)	2 (2)	4 (3.6)	0.37 (0.32)	0.61 (0.52)	0.35 (0.32)	0.27 (0.24)	0.022	0.360	0.132	−0.035
bIGS‐24d	59	4 (4)	3 (3)	3 (3)	4 (3.6)	0.50 (0.57)	0.50 (0.46)	0.53 (0.61)	0.73 (0.54)	−0.398	0.185	0.227	0.124
bIGS‐BP1	59	2 (2)	2 (2)	2 (2)	2 (1.9)	0.12 (0.11)	0.61 (0.49)	0.45 (0.46)	0.10 (0.15)	0.094	−0.133	0.063	−0.084
cGS‐14	53	5 (5)	4 (3.9)	4 (3.9)	5 (4.9)	0.93 (0.75)	0.55 (0.58)	0.67 (0.64)	0.73 (0.73)	−0.203	0.079	−0.010	0.022
Mean		4.1 (4.1)	3.4 (3.3)	3.9 (3.8)	4.1 (3.8)	0.55 (0.54)	0.51 (0.51)	0.57 (0.59)	0.52 (0.50)	0.010	0.037	0.059	−0.044
± SE		0.4	0.4	0.4	0.4	0.07 (0.06)	0.05 (0.05)	0.05 (0.04)	0.08 (0.07)
MHC DRB1		4	3 (2.9)	6 (5.5)	3 (2.9)	0.50 (0.49)	0.28 (0.25)	0.53 (0.52)	0.27 (0.31)	0.004	−0.090	−0.007	0.164

aJones et al. (2005), bMay et al. (1997), cStevens et al. (1997), HOut of Hardy‐Weinberg equilibrium (P = 0.01).

Test	Model	AV	SE	ESC	AVC
H E excess	SMM	7:5	6:6	3:9	6:6
		0.733	0.733	0.092	0.518
	IAM	3:9	1:11	1:11	5:7
		0.176	0.034	0.002	0.233
Mode shift		L‐shaped	L‐shaped	L‐shaped	L‐shaped