MHC-associated Baylisascaris schroederi load informs the giant panda reintroduction program.

Reintroducing captive giant pandas (Ailuropoda melanoleuca) to the wild is the ultimate goal of their ex situ conservation. Choosing higher fitness candidates to train prior to release is the first step in the giant panda reintroduction program. Disease resistance is one important index of individual fitness and presumed to be related to variation at major histocompatibility complex genes (MHC). Here, we used seven polymorphic functional MHC genes (Aime-C, Aime-I, Aime-L, Aime-DQA1, Aime-DQA2, Aime-DQB1 and Aime-DRB3) and estimate their relationship with Baylisascaris schroederi (Ascarididae) infection in giant panda. We found that DQA1 heterozygous pandas were less frequently infected than homozygotes. The presence of one MHC genotype and one MHC allele were also associated with B. schroederi infection: Aime-C*0203 and Aime-L*08 were both associated with B. schroederi resistance. Our results indicate that both heterozygosity and certain MHC variants are important for panda disease resistance, and should therefore be considered in future reintroduction programs for this species alongside conventional selection criteria (such as physical condition and pedigree-based information).

Introduction

Ex situ conservation has been conducted for the endangered giant panda for half a decade and great progress has been made: the number of giant pandas in captivity has increased to 600 (Jia, 2019). One important goal of ex situ conservation is to reintroduce captive individuals to the wild, ensure they survive and breed, and therefore reinforce the wild population and reduce the threat of extinction (Frankham et al., 2009). The reintroduction program for giant panda begun in 2003 and has proven to be labor intensive and costly (Zhang, 2013). In total, eleven pandas have been released to the wild to 2018 (https://www.pandasinternational.org). The typical procedure – employed since 2010 – is to identify pandas for release at age 2–3, and then train them with their mothers prior to reintroduction. Training consists of looking for bamboo, water and shelter, and avoiding animals which could cause them harm in a semi-wild environment (Zhang, 2013). Due to the high intensity of these activities, choosing appropriate (high quality) individuals to train and release is the first step in the panda reintroduction program, with individuals usually ranked on the basis of physical condition (weight, size, teeth, blood index, medical history) and pedigree information (Zhang, 2013). Parasitic disease is a major threat to wild pandas (Zhang et al., 2008), and so individual measures of disease resilience should also be considered alongside these factors.

The ascaridoid Baylisascaris schroederi (B. schroederi) is the most common intestinal parasite of the giant panda, and specific to the species; it has been reported in multiple wild and captive populations (Hu, 2001; Zhang and Wei, 2006). Infection with this parasite can cause intestinal obstruction, inflammation and death (Hu, 2001). Wild panda deaths from parasitic disease (mainly B. schroederi) have increased significantly from the 1970s–2000s; as many as 50% of deaths resulted from this parasite at the beginning of the 2000s (Zhang et al., 2008). Therefore, individuals with strong resistance (immunocompetence) to B. schroederi should be considered valuable options in reintroduction planning. According to parasite screening data from the Chinese Research and Conservation Center for Giant Panda (CRCCGP), wide inter-individual variation in parasite infestation was observed amongst pandas housed under the same conditions. We therefore hypothesize a relationship between immunogenetic background and B. schroederi infection in giant pandas.

Major histocompatibility complex genes (MHC) encode molecules that recognize and present foreign peptides to T cells and subsequently initiate an immune reaction (Klein, 1986). At least two non-mutually exclusive modes have been proposed to explain an MHC-mediated co-evolutionary arms race between hosts and pathogens: (i) Heterozygote advantage, whereby MHC heterozygotes recognize a broader range of antigens, and thus suffer a lower level of infection, compared to MHC homozygous individuals (Doherty and Zinkernagel, 1975). MHC heterozygote advantage has received convincing evidence from controlled experimental infection studies (Carrington et al., 1999; McClelland et al., 2003; Penn et al., 2002) and has also been supported in studies of natural populations, such as Chinook salmon Oncorhynchus tshawytscha (Evans and Neff, 2009), water voles Arvicola terrestris (Oliver et al., 2009), striped mouse Rhabdomys pumilio (Froeschke and Sommer, 2005). Other studies have failed to find support for this hypothesis (Deter et al., 2008; Langefors et al., 2001; Meyer-Lucht and Sommer, 2005; Schad et al., 2005, 2012; Tollenaere et al., 2008; Zhang et al., 2015). (ii) Negative frequency-dependent selection, where the dynamic relationship between the frequencies of MHC alleles, and the parasites targeted by them, results in high levels of host MHC diversity (Takahata and Nei, 1990). Evidence in line with MHC-based frequency-dependent selection has been reported for natural populations. In some cases, certain MHC alleles were correlated with increased resistance, such as Soay sheep Ovis aries (Paterson et al., 1998), hairy-footed gerbil Gerbillurus paeba (Harf and Sommer, 2005) and water vole Arvicola scherman (Tollenaere et al., 2008). Susceptible MHC alleles were observed in the Neotropical bat Noctilio albiventris (Schad et al., 2012). Both resistance and susceptibility MHC alleles were documented in striped mouse Rhabdomys pumilio (Froeschke and Sommer, 2005), yellow-necked mice Apodemus flavicollis (Meyer-Lucht and Sommer, 2005), mouse lemurs Microcebus murinus (Schad et al., 2005), bank voles myodes glareolus (Deter et al., 2008) and Atlantic salmon Salmo salar (Langefors et al., 2001).

MHC genes has been comprehensively studied in giant panda, in which six functional MHC class II genes (Aime-DRA, Aime-DRB3, Aime-DQA1, Aime-DQA2, Aime-DQB1 and Aime-DQB2) (Wan et al., 2009, 2011) and four classical MHC class I genes (Aime-C, Aime-F, Aime-I and Aime-L) (Zhu et al., 2012, 2013b) have been characterized, and mature genotyping methods developed. MHC-parasite studies in wild giant panda using Aime-DRB1, Aime-DQA1 and Aime-DQA2 found a putative susceptibility allele (Aime-DRB1*10), which was associated with parasite infection (Zhang et al., 2015). This previous study did not include the polymorphic Aime-DQB1 gene, nor any MHC class I genes, which might lead to incomplete conclusions as shown in our previous study (Zhu et al., 2019). Here, by genotyping multiple MHC genes, we aim to investigate heterozygote advantage and allele-based selection processes in respect of B. schroederi infection using all polymorphic functional MHC genes known in giant panda. The results will inform the selection of individuals for the giant panda reintroduction program.

Materials and methods

Sample collection

We obtained fecal and blood samples from 75 giant pandas housed in Dujiangyan (DJY, N = 21), Bifengxia (BFX, N = 22) and Wolong breeding bases (WL, N = 32) of CRCCGP respectively. Blood samples were collected during a routine medical examination and preserved in liquid nitrogen for DNA extraction. We obtained permission from the CRCCGP to collect all the samples. Blood samples were obtained with permission from the China Giant Panda Protection and Management Office during routine medical examinations.

Fecal samples were used for parasite analysis. Since anti-parasitic treatments are conducted every month for each panda in all three breeding facilities, but with different schedules, we collected fecal samples on the day before treatment. To avoid possible variation produced by sampling date, we collected all samples across a three-day period in August. We collected 6 samples for each panda (N = 450 fecal samples in total), which were combined prior to parasite detection to avoid possible variation in parasite distribution among samples. Fresh fecal samples (<12 h) were collected in the morning and transported to the laboratory under 4 °C, within 2 h. We obtained permission from the CRCCGP to collect fecal samples and confirmed that we did not impact the animal during the sampling. Biological samples were obtained according to the guidelines and approval of the Animal Ethics Committee of Sichuan Provincial Academy of Natural Resource Sciences (171009–1).

Parasite detection

We quantified B. schroederi infection followed methods described previously with minor modification (Zhang, 2015). Briefly, 40 g of fecal material was transferred to a beaker containing 300 ml sterile water. The mixture was filtered through 425 and 180 μm mesh sequentially to remove residual bamboo. After a 30 min incubation at room temperature, the supernatant was removed and sediment transferred to a 50 ml Eppendorf tube. The sediment was centrifuged at 5000 rpm/min for 20 min and the supernatant immediately removed. We added 30 ml saturated NaCl solution to the pellet and mixed thoroughly. Three 1-ml aliquots of each sample were examined on a compound microscope using a McMaster slide under 10 × magnification to detect helminths based on their morphological characteristics. If no B. schroederi were detected from three samples, another 5 samples were examined. The number of eggs per gram of feces (EPG) was calculated as the parameter of infection intensity.

DNA extraction and MHC genotyping

Genomic DNA was extracted from blood samples as described in Zhu et al. (2013). We examined seven MHC loci including three class I (Aime-C, Aime-I and Aime-L) and four class II (Aime-DQA1, Aime-DQA2, Aime-DQB1 and Aime-DRB3). Three loci that were previously found to be monomorphic (Aime-F, Aime-DQB2 and Aime-DRA (Chen et al., 2010; Zhu et al., 2013a);) were excluded from this study.

For MHC class I genes, we adopted our previous locus-specific amplification protocol with the following modifications (Zhu et al., 2013b). We amplified a fragment comprising exon2, inron2 and exon3 of each gene and cloned PCR products into DH5 competent cells (TaKaRa, Ltd, Dalina, China). Then we selected 6 positive clones per amplicon to sequence and determine the exon 2–3 genotypes. For MHC class II genes, PCR amplification and genotyping were performed according previous methods (Chen et al., 2013). The primer sets and PCR amplification conditions are presented in Table S1.

Data analysis

Associations between MHC genotypes and infection status were examined using linear regression in R 3.5.3 (R Development Core Team, 2011). The goal of our analyses was to test whether heterozygosity per se, particular genotypes, or particular alleles were associated with infection status and parasite load in giant pandas. Data exploration was first conducted to check data distribution and to investigate relationships among individual-level predictor variables.

Our parasite data showed that 38 pandas (nearly 50%) had 0 EPG and a few pandas had more than 100 EPG, a distribution suggesting excess zeros and overdispersion. We compared several count models using the rootgram function in the countreg package (Kleiber and Zeileis, 2016) and chose to fit hurdle negative binomial models to accommodate both excess zeros and overdispersion. A hurdle model is a type of linear regression that includes two parts, one for the zero count and one for the positive counts. The first part is a binary logistic model, which fits predictors of giant panda infection probability; the second part is a truncated negative binomial model, which fits predictors of a positive count, and thus the determinants of infection intensity. Hurdle models were fitted using the pscl package (Zeileis et al., 2008) for R.

We tested the effect of MHC heterozygosity at four levels, namely: (1) multilocus heterozygosity (MLH), i.e., the proportion of all genotyped MHC loci that were heterozygous; (2) multilocus heterozygosity at the three MHC class I loci (MLHI); (3) multilocus heterozygosity at the four MHC class II loci (MLHII) and (4) observed heterozygosity of each of the seven individual MHC loci (H), coded as 1 for heterozygote and 0 for homozygote. In addition to the genetic variable, each model also included three non-genetic predictors: sampling location (three levels), sex and age. We did not fit any interactions between predictor variables. In order to evaluate the effect of the three non-genetic predictors, we generated all submodels from the full model using the dredge function implemented in the package MuMIn (Barton and Barton, 2013). Only location appeared in the top model, which we interpreted as indicating little impact of age and sex (Table S2). We therefore excluded sex and age from subsequent analyses.

We hypothesized that, independent of the effect of heterozygosity, specific MHC genotypes or alleles could be associated with infection status and load. We therefore tested whether genotype and allele frequencies of each locus had an effect on infection with B. shroederi. We focused on those genotypes and alleles that were observed in at least 10 of the genotyped individuals (total N = 75; we had insufficient sample size to examine the effects of rarer genotypes or alleles). We modelled the effect of each genotype or allele on parasite infection by coding the data with 0/1 for representing absence/presence of the genotype or allele. In addition to the genetic variable, these models also included sampling location as a predictor.

Model selection was based on corrected Akaike Information Criteria for small sample size (AICC (Burnham et al., 2011),). For locus-level heterozygosity, genotype and allele effect models, we ranked AICC values of all models along with a base model that excluded genetic data (intercept and location only). Based on the ranked models, if the best model showed a ΔAICC ≥2 relative to other models and superior to the base model, we considered it evidence that the heterozygosity metric, genotype or allele in question influences B. shroederi in giant pandas (following (Grueber et al., 2013; Sepil et al., 2013)).

To visualize our data, we evaluated the 95% confidence interval (CI) of observed infection probability using the binom.confint function in the package binom (Dorai-Raj and Dorai-Raj, 2009). The 95% CI of observed EPG was produced by resampling 10, 000 × with replacement. The 95% CIs of predicted infection probability and EPG were obtained using parametric bootstrapping (10,000 × ) of the corresponding model coefficients.

Results

The overall prevalence of B. shroederi was 49.3% (37/75) and there was wide variation in abundance, from 0 to 238 EPG (Fig. S1). The prevalence and intensity of infection varied greatly across the three sampling locations, with the highest intensity at BFX (39.7 EPG with 63.6% prevalence) and lowest in WL (2.4 EPG with 25% prevalence; top model in Table S2). Of the three non-genetic predictors we tested (location, sex, age), location was the most compelling predictor of parasite load (ΔAICC < 2; Table S2).

Effects of heterozygosity on B. Shroederi infection

For overall heterozygosity, models that excluded MLH (base model) showed superior AICC support relative to a model that included MLH (Table S2), suggesting that there is no evidence for an effect of overall MHC heterozygosity on giant panda infection status and parasite load. For multilocus heterozygosity by MHC class, there was no difference between base model and models that included three MHC class I and four MHC class II genes (Table S3).

When examining the effect of heterozygosity at each of the seven individual MHC loci, only DQA1 showed substantially greater evidence for an effect of heterozygosity on infection status/parasite load, relative to heterozygosity models for other genes and relative to the base model (which excludes genetic data) (Table 1). The negative coefficient for zero hurdle model suggests that DQA1 heterozygotes had reduced parasite infection probability relative to homozygotes. (Table 1). The observed and predicted infection probability for DQA1 homozygotes was more than 1.6 times and 1.9 times that of DQA1 heterozygotes, respectively (Fig. 1a, observed: 85% vs 53%, predicted: 95% vs 49%), although our hurdle models had poor predictive power for this trait (wide standard error; Table 1). A strong, negative trend was also seen for observed parasite load. The observed parasite load for DQA1 homozygous individuals was more than 1.2 times that of DQA1 heterozygotes (Fig. 1b, observed: 38.7 vs 32.3 EPG). However, the model predicted the opposite trend for DQA1 homozygotes and heterozygotes (Fig. 1b, predicted: 13.1 vs 15.6 EPG). An absence of DQA1 heterozygote infected pandas at DJX may account for the differences between the observed and predicted patterns.

Table 1

Effect of heterozygosity (H) at seven MHC genes on B. schroederi infection in giant panda (N = 56).

Model1	βzero (±SEβzero)2	βcount (±SEβcount)3	AICC	ΔAICCC	wi	N4	infection %5	NE6	EPG7
Base + HDQA1	−2.887 (1.357)	0.178 (0.885)	309.317	0	0.603	49	53.1	26	32.3
Base + HI	0.692 (0.691)	0.785 (0.658)	313.075	3.758	0.092	34	61.8	21	44.0
Base + HDQB1	−0.336 (0.691)	−0.728 (0.642)	313.740	4.423	0.066	35	54.3	19	34.5
Base			313.753	4.436	0.066	56*	57.1*	32*	33.5*
Base + HDRB3	−1.016 (1.001)	0.380 (0.822)	314.079	4.762	0.056	47	53.2	25	33.4
Base + HC	−0.482 (0.801)	0.540 (0.690)	314.475	5.158	0.046	43	55.8	24	34.9
Base + HL	0.518 (0.675)	0.229 (0.667)	314.688	5.371	0.041	23	65.2	15	29.1
Base + HDQA2	0.132 (0.678)	−0.155 (0.638)	315.296	5.979	0.030	21	57.1	12	42.3

Note: 1 Hurdle model incorporating a binomial distribution for the zero part and a negative binomial distribution for the count part, where the base model includes only location as a predictor. Each heterozygosity model includes location plus a 0/1 binary predictor for the absence/presence of heterozygosity at the specified MHC gene.

2 Effect size coefficient (βzero) and its standard error (SEzero) for the zero component of the hurdle model.

3 Effect size coefficient (βcount) and its standard error (SEβcount) for the truncated count component of the hurdle model.

4 Number of giant pandas that were heterozygous at the specified locus.

5 Infection probability of giant pandas that were heterozygous at the specified locus.

6 Number of infected giant pandas that were heterozygous at the specified locus.

7 EPG of infected giant pandas that were heterozygous at the specified locus.

* Overall sample size, infection probability and EPG count for all loci and all pandas.

Fig. 1

Effect of DQA1 heterozygosity on (a) infection probability, (b) EPG of B. schroederi. The Bars show the 95% CI.

Effects of MHC genotype on B. Shroederi infection

When examining the effects of specific MHC genotypes on parasite infection, strong relationships were only observed for genotypes of the three MHC class I loci (Aime-C, Aime-I and Aime-L). None of the genotypes observed at the four MHC class II loci (Aime-DQA1, Aime-DQA2 and Aime-DQB1 and Aime-DRB3) were strong predictors of parasite infection (Table 2). Three genes showed one to two genotypes that were predictive of parasite outcomes: Aime-C (Aime-C*0305 and Aime-C*0203), Aime-I (Aime-I*0204) and Aime-L (Aime-L*0202), even accounting for any effect of sampling location. We found that pandas containing the genotypes Aime-C*0305, Aime-C*0203 and Aime-I*0204 were less likely to be infected (negative coefficient in Table 2, Fig. 2a and b). However, genotype Aime-L*0202 provided a disadvantage for parasite infection status, as indicated by a positive model coefficient (βzero ± SEβzero = 1.030 ± 0.807, Table 2, Fig. 2a and b). The observed parasite loads for infected individuals with these four genotypes were all lower than the loads of infected individuals with other genotypes (Fig. 2c). However, the predictors, Aime-C*0305, Aime-I*0204 and Aime-L*0202 showed the opposite trend compared to the observed data. There were no infected pandas with the above three MHC types at BFX, which may explain the differences between the observed and predicted patterns (Fig. 2d).

Table 2

Effects of genotypes at seven MHC genes on B. schroederi infection in giant panda. Only genotypes observed in ≥10 individuals were examined.

Model1	βzero (±SEβzero)2	βcount (±SEβcount)3	AICCC	ΔAICCC	wi	N4	infection %5	NE6	EPG7
Aime-C
Base + 0305	−0.203 (0.733)	1.814 (0.956)	360.773	0.000	0.553	11	0.455	5	16.8
Base + 0203	−0.938 (0.797)	−1.107 (1.083)	362.560	1.787	0.226	10	0.300	3	17.7
Base			363.674	2.902	0.130	75*	0.493*	37*	29.1*
Base + 0303	0.086 (0.728)	−0.631 (1.009)	364.382	3.609	0.091	11	0.455	5	39.8
Aime-I
Base + 0204	−0.719 (0.813)	3.726 (1.274)	343.503	0.000	0.969	10	0.300	3	24.3
Base + 0202	−0.177 (0.587)	−1.002 (0.761)	351.515	8.013	0.018	21	0.524	11	10.5
Base			352.047	8.545	0.014	68*	0.529*	36*	29.9*
Aime-L
Base + 0202	1.030 (0.807)	2.459 (1.300)	352.481	0.000	0.722	10	0.600	6	14.7
Base + 0303	−0.915 (0.652)	−0.525 (0.877)	355.240	2.759	0.182	16	0.375	6	34.5
Base			356.504	4.023	0.097	73*	0.493*	30*	29.9
Aime-DQA1
Base + 0304	0.627 (0.846)	−1.009 (0.909)	331.360	0.000	0.417	10	0.600	6	6.8
Base			331.695	0.335	0.353	61*	0.557*	34*	31.6*
Base + 0105	0.117 (0.806)	−0.647 (0.884)	332.544	1.183	0.231	11	0.455	5	36.0
Aime-DQA2
Base			334.096	0.000	0.460	62*	0.548*	34*	31.6*
Base + 0102	0.457 (0.657)	0.094 (0.680)	335.047	0.951	0.286	18	0.611	11	45.8
Base + 0101	−0.244 (0.612)	0.185 (0.655)	335.276	1.179	0.255	40	0.525	21	27.0
Aime-DQB1
Base + 0203	1.128 (0.905)	−0.616 (0.748)	326.417	0.000	0.396	10	0.700	7	11.6
Base			327.328	0.911	0.251	60*	0.550*	33*	32.5*
Base + 0204	−0.624 (0.755)	−0.584 (0.805)	327.584	1.167	0.221	13	0.385	5	28.8
Base +0404	−0.022 (0.690)	−0.162 (0.698)	328.627	2.210	0.131	16	0.500	8	34.0
Aime-DRB3
Base +0802	0.323 (0.785)	−0.858 (0.711)	316.982	0.000	0.504	13	0.538	7	26.6
Base			317.016	0.033	0.496	58*	0.552*	32*	33.5*

Note: 1 Hurdle model incorporating a binomial distribution for the zero part and a negative binomial distribution for the count part, where the base model includes only location as a predictor. Each genotype model includes location plus a 0/1 binary predictor for the absence/presence of the specified genotype at the specified MHC gene. Each genotype is specified by two alleles at a locus, i.e, 0305 at Aime-C consists of the alleles Aime-C*03 and allele Aime-C*05.

2 Effect size coefficient (βzero) and its standard error (SEzero) for the zero component of the hurdle model.

3 Effect size coefficient (βcount) and its standard error (SEβcount) for the truncated count component of the hurdle model.

4 Number of giant pandas containing the specified genotype. Genotypes observed in fewer than 10 individuals (low genotype frequency) were not included in this analysis.

5 Infection probability of giant pandas that had the specified genotype.

6 Number of infected giant pandas that had the specified genotype.

7 EPG of infected giant pandas that had the specified genotype.

* Overall sample size, infection probability and EPG count for all genotypes at each given locus.

Fig. 2

Effects of MHC genotypes on (a–b) infection probability, (c–d) EPG of B. schroederi using important predictors. The Bars show the 95%CI.

Effects of MHC alleles on B. Shroederi infection

When examining the effects of specific alleles at seven MHC loci, we found compelling associations with parasite loads for seven alleles at four loci (Aime-C*05, Aime-I*04, Aime-L*01, Aime-L*08, Aime-DQA1*03, Aime-DRB3*01 and Aime-DRB3*03, Table 3). Except for Aime-L*01, Aime-DQA1*03 and Aime-DRB3*03, all other five alleles showed a defensive advantage to pandas with these alleles (negative coefficient, Table 3). For most alleles, model predictions were consistent with observed data, except DQA1*03 infection probability (Fig. 3). Pandas carrying Aime-L*08 showed lower infection probability and reduced parasite load relative to pandas without Aime-L*08 (Fig. 3a and b). Although Aime-C*05, Aime-I*04 and Aime-DRB3*01 provided an apparent defense against infection, pandas that were infected may suffer greater loads than those without these alleles (Table 3). The alleles Aime-L*01, and Aime-DRB3*03 showed the opposite trend: pandas with Aime-L*01 and Aime-DRB3*03 were easily infected, but the infections showed reduced load compared to pandas without these alleles (Fig. 3c and d).

Table 3

Effects of allelesat seven MHC genes on B. schroederi infection.

Model1	βzero (±SEβzero)2	βcount (±SEβcount)3	AICCC	ΔAICCC	wi	N4	infection %5	NE6	EPG7
Aime-C
Base + 05	−0.515 (0.547)	1.484 (0.676)	358.391	0.000	0.761	25	0.400	10	41.9
Base + 02	−0.654 (0.568)	−0.72 (0.725)	362.549	4.158	0.095	23	0.435	10	27.1
Base			363.674	5.283	0.054	75*	0.493*	37*	29.1*
Base + 06	0.300 (0.533)	0.551 (0.875)	363.993	5.601	0.046	27	0.556	15	36.6
Base + 03	−0.104 (0.523)	−0.695 (0.868)	364.142	5.751	0.043	45	0.467	21	20.0
Aime-I
Base + 04	−0.277 (0.715)	2.980 (0.936)	339.529	0.000	0.969	12	0.417	5	50.8
Base + 01	1.004 (0.773)	−1.874 (0.758)	347.070	7.540	0.022	11	0.636	7	12.3
Base + 03	0.881 (0.704)	−0.750 (0.930)	350.912	11.383	0.003	14	0.714	10	44.2
Base + 07	−0.212 (0.725)	0.949 (0.904)	351.889	12.359	0.002	11	0.545	6	42.2
Base			352.047	12.518	0.002	68*	0.529*	36*	29.9*
Base + 02	0.156 (0.644)	−0.517 (0.837)	352.725	13.195	0.001	53	0.528	28	29.3
Aime-L
Base + 01	0.813 (0.730)	−1.354(0.826)	354.035	0.000	0.373	12	0.667	8	25.0
Base + 08	−0.198 (0.634)	−1.661(0.808)	354.241	0.205	0.336	16	0.438	7	17.7
Base			356.504	2.469	0.108	73*	0.493*	36*	29.9*
Base + 02	0.410 (0.574)	0.641 (1.005)	356.640	2.605	0.101	22	0.545	12	30.4
Base + 03	−0.365 (0.536)	−0.046 (0.843)	357.074	3.039	0.082	44	0.455	20	30.5
Aime-DQA1
Base + 03	−0.438 (0.629)	1.317 (0.682)	328.741	0.000	0.501	28	0.571	16	32.9
Base + 05	−0.040 (0.628)	−1.115 (0.776)	330.954	2.213	0.166	22	0.500	11	33.2
Base + 04	0.004 (0.642)	−0.222 (0.567)	331.139	2.398	0.151	21	0.571	12	27.4
Base			331.695	2.954	0.114	61*	0.557*	34*	31.6*
Base + 01	−0.309 (0.605)	0.193 (0.673)	332.728	3.987	0.068	31	0.548	17	37.5
Aime-DQB1
Base + 04	−0.907 (0.633)	−0.582 (0.640)	325.790	0.000	0.464	36	0.444	16	33.7
Base			327.328	1.538	0.215	60*	0.550*	11	32.5*
Base + 03	0.655 (0.700)	−0.390 (0.622)	327.415	1.625	0.206	17	0.647	33*	29.5
Base + 02	0.238 (0.611)	0.149 (0.702)	328.567	2.777	0.116	33	0.576	19	22.9
Aime-DRB3
Base + 01	−0.022 (0.641)	1.319 (0.634)	313.965	0.000	0.447	24	0.542	13	39.8
Base + 03	1.440 (0.840)	−0.458 (0.613)	314.586	0.621	0.328	15	0.733	11	29.6
Base			317.016	3.051	0.097	58*	0.552*	32*	33.5*
Base + 08	−0.212 (0.677)	−0.726 (0.706)	317.350	3.384	0.082	19	0.526	10	36.5
Base + 02	0.057 (0.650)	−0.068 (0.654)	318.509	4.544	0.046	34	0.559	19	33.8

Note: 1 Hurdle model incorporating a binomial distribution for the zero part and a negative binomial distribution for the count part, where the base model includes only location as a predictor. Each model includes location plus a 0/1 binary predictor for the absence/presence of each specified allele at each given locus. Data for DQA2 are not shown as the model did not converge.

2 Effect size coefficient (βzero) and its standard error (SEzero) for the zero component of the hurdle model.

3 Effect size coefficient (βcount) and its standard error (SEβcount) for the truncated count component of the hurdle model.

4 Number of giant pandas containing the specified allele. Alleles observed in fewer than 10 individuals (low allele frequency) were not included in this analysis.

5 Infection probability of giant pandas that had the specified allele.

6 Number of infected giant pandas that had the specified allele.

7 EPG of infected giant pandas that had the specified allele.

* Overall sample size, infection probability and EPG count for all alleles given a locus.

Fig. 3

Effects of MHC alleles on infection probability (a and b) and EPG (c and d) of B. schroeder.i. The bars show the 95%CI.

Discussion

The ultimate goal of giant panda ex situ conservation is to reintroduce captive pandas to ensure they live and breed in the wild, recovering the species from endangerment. To increase survival after release, it is wise to choose individuals with greater resilience to parasitic diseases commonly faced in the wild. Despite large resources invested in the reintroduction program (Zhang, 2013), this study represents the first complete examination into the effect of immunogenetic background on panda fitness. Overall, the results showed that both heterozygosity and certain MHC types are associated with infection status and/or parasite load in the giant panda.

Variation in infection probability and parasite load of B. schroederi were associated with DQA1 heterozygosity. The lowest infection probability was observed in DQA1-heterozygote individuals, which is consistent with our previous study of giant panda mate choice (Zhu et al., 2019). That previous work concluded that females preferred DQA1-heterzygote males over DQA1-homozygote males (Zhu et al., 2019). The current findings support the contention that female pandas prefer DQA1-heterozygoty males due to their superior disease resistance. In wild pandas however, DQA1 heterozygosity was not related to infection status nor parasite load; the difference between captive and wild pandas could be due to significant heterozygote deficiency observed at DQA1 in wild pandas (Zhang et al., 2015). In addition to effects of heterozygosity, our study also found genotypes and alleles that appear to predict infection status and parasite load of B. schroederi in pandas. Among these, Aime-C*0203, showed the most compelling effect, with both lower infection probability and reduced parasite load.

Deworming is common in captive giant panda populations, but environmental control of parasites is currently infeasible in wild populations. Therefore, captive animals may suffer more from B. schroederi when released to the wild, than they do in captivity. Thus, choosing pandas with higher disease resistance for reintroduction may benefit the health, welfare and survival of released individuals. Our finding that MHC variation is associated with infection status and/or parasite load in giant panda suggests that considering variation at this gene region may improve survival following reintroduction of giant panda individuals. Whether to target conservation actions (such as breeding or release) towards individuals carrying particular genotypes or alleles can be a contentious issue in conservation programs, and should be considered carefully in order to avoid losses of overall genetic diversity (Ralls et al., 2000). This issue was addressed by conservationists working to preserve the California condor, which suffers from chondrodystrophy, a recessive genetic disorder that severely affects survival in the conservation breeding program (Ralls et al., 2000). The occurrence of chondrodystrophy could be minimized at little cost to maintenance of genetic diversity if genetic management is carefully planned and implemented (Ralls et al., 2000). In the case of the panda breeding program, it is possible that individuals carrying resistant genotypes might be related to one another, therefore kinship should be simultaneously evaluated when selecting animals for breeding or reintroduction, to avoid inbreeding.

Within our dataset, we found some support for both the heterozygote advantage and frequency dependent hypotheses (selection favoring particular alleles). Interestingly, these hypotheses were supported with different MHC loci. For example, DQA1 showed evidence for heterozygosity advantage, but showed no compelling effects of particular alleles in the captive population, consistent with a study of wild pandas (Zhang et al., 2015). These results emphasize the complexity of evolutionary and population processes driving MHC diversity. Furthermore, individuals carrying the Aime-C*05 and Aime-I*04 alleles were only observed as heterozygotes, and the two genotypes Aime-C*0305 and Aime-I*0204 were detected as predictors of B. schroederi infection. Disease resistance associated with genotypes Aime-C*0305 and Aime-I *0204 may therefore be due to inheritance of the alleles Aime-C*05, Aime-I*04 respectively under the good genes hypothesis; our data do not allow us to differentiate heterozygote advantage from allele-based selection (such as frequency dependent selection). We have previously reported that the alleles Aime-C*02, Aime-C*03, and Aime-I*02 are at high frequencies in the wild population (Zhu, 2012). Releasing Aime-C*0305 or Aime-I *0204 individuals would therefore represent genotypic combinations that may be particularly beneficial, as they could inherit the fitness advantage of Aime-C*05 and Aime-I*04.

To our knowledge, this is the first study to incorporate all functional MHC loci into a single analysis of the role of these genes on infection status and intensity in a non-model species. Those genes that showed a relationship with parasite resistance (Aime-DQA1, Aime-C and Aime-I) may be more important for B. schroederi resistance than other MHC genes, but it is also possible that other MHC genes and alleles are essential for resilience to other, unstudied parasites or diseases. Furthermore, we note that two of these MHC classical I loci, Aime-C and Aime-I have been implicated in female panda mate choice (Zhu et al., 2019). Taken together, these findings reiterate the conclusions of others (e.g. (Feng et al., 2009; Kamiya et al., 2014; Zhu et al., 2019) that studying only a few MHC loci may not give the full picture of the effect of these genes on fitness outcomes, and that all polymorphic, functional MHC loci should be addressed in future studies in other species.

We acknowledge that some of the effects we report were estimated with poor precision since the standard errors of relevant coefficients were large (Tables 1–3, Fig. 1–3). We also note that breeding pandas with particular MHC genotypes might also be difficult as it would depend on which breeding individuals are available, and the stochastic nature of Mendelian inheritance cannot ensure particular offspring genotypes. Nevertheless, our data do suggest that incorporating MHC genotype information into the giant panda's reintroduction program, alongside typical parameters such as kinship and age, could decrease the probability of B. schroederi infection. As the program progresses over time, and more data are collected, we hope that the patterns we observed will become increasingly refined.

Declaration of competing interests

All the authors declare that they have no known competing interests.