No evidence for MHC‐based mate choice in wild giant pandas

Abstract Major histocompatibility complex genes (MHC), a gene cluster that controls the immune response to parasites, are regarded as an important determinant of mate choice. However, MHC‐based mate choice studies are especially rare for endangered animals. The giant panda (Ailuropoda melanoleuca), a flagship species, has suffered habitat loss and fragmentation. We investigated the genetic variation of three MHC class II loci, including DRB1, DQA1, and DQA2, for 19 mating‐pairs and 11 parent‐pairs of wild giant pandas based on long‐term field behavior observations and genetic samples. We tested four hypotheses of mate choice based on this MHC variation. We found no supporting evidence for the MHC‐based heterosis, genetic diversity, genetic compatibility and “good gene” hypotheses. These results suggest that giant pandas may not use MHC‐based signals to select mating partners, probably because limited mating opportunities or female‐biased natal dispersal restricts selection for MHC‐based mate choice, acknowledging the caveat of the small sample size often encountered in endangered animal studies. Our study provides insight into the mate choice mechanisms of wild giant pandas and highlights the need to increase the connectivity and facilitate dispersal among fragmented populations and habitats.

of mature sperm can be identified by females, which remains to be proved and is not our main concern here. Thus, the offspring attain additional fitness benefits in precopulatory mate choice through MHC-mediated selective mechanisms in terms of immune resistance to disease and adaptations to the changing environment, rendering sexual selection evolutionarily favorable.
Previous research on MHC-based mate choice has identified four main genetic mechanisms underlying precopulatory mate choice in vertebrates (Brown, 1997;Hamilton & Zuk, 1982;Kempenaers, 2007;Møller & Alatalo, 1999;Penn, Damjanovich, & Potts, 2002;Puurtinen, Ketola, & Kotiaho, 2005). First, mate choice is mediated through heterosis, which supposes that the heterozygotes of a targeted MHC locus are preferred over homozygotes. A direct benefit of heterozygosity in the MHC genes is the reduced risk of disease in offspring. Heterozygosity also provides indirect benefits by allowing the individual to be a carrier of rare alleles (Hughes & Nei, 1988).
Second, mate choice occurs based on MHC genetic diversity, measured by the number of variable sites at a targeted MHC locus (Lenz, Wells, Pfeiffer, & Sommer, 2009). Third, animals prefer to mate with MHC-dissimilar partners (i.e., genetic compatibility hypothesis), which is expected to produce offspring with higher fitness (Raveh et al., 2014;Sommer, 2005). In addition, if the choice of MHC gene is based on the immune effect, the partner's MHC allele is important. The optimizing hypothesis assumes that choosing a mate with the most appropriate number of MHC diversity or dissimilarity, but not the maximum, will benefit future generations, that is, preferred in some animals (Milinski, 2006;Wegner, Kalbe, Kurtz, Reusch, & Milinski, 2003), which is similar to the genetic compatibility hypothesis. Fourth, the "good gene" hypothesis emphasizes the choice for specific MHC alleles. For example, sires of wild and captive tucotucos were shown to carry distinctive alleles (Cutrera, Fanjul, & Zenuto, 2012).
The giant panda (Ailuropoda melanoleuca) is one of the most threatened animals, although the IUCN recently downlisted it from "endangered" to "vulnerable" based on population and habitat recovery (Swaisgood et al., 2016). The 4th national survey of giant pandas estimated the global population size at 1,864 individuals living in 33 fragmented populations (State Forestry Administration, 2015). A recent genetic study also found that the level of inbreeding in wild giant pandas was greater than expected (Hu et al., 2017). However, multiple studies have shown that giant pandas still harbor relatively high genetic diversity and evolutionary potential (Zhang et al., 2007), suggesting that some mechanisms may maintain genetic variation in the population. Studies have reported that MHC-based mate choice mechanisms function in shaping genetic diversity. For example, Santos, Michler, and Sommer (2017) found that the MHC-based mechanism in mate choice may allow raccoons to overcome inbreeding constraints during population expansion.
Wild giant pandas are rare and elusive and live in complex habitats, making direct field observations extremely difficult. Knowledge of wild giant panda reproduction is inferred from three long-term radio/GPS-telemetry studies (Nie, Swaisgood, Zhang, Liu, & Wei, 2012;Pan et al., 2014;Schaller, Hu, Pan, & Zhu, 1985). The courting and mating behaviors of giant pandas often occur from the end of February to the end of April. At first, males gather around an estrous female in a mating site. Each mating site comprises one female and at least one male. In most cases, males determine dominance through combat, and the winner receives access to the female. The female usually mates with the male winner, but occasionally mating opportunities will be granted to other males. A previous study on captive giant pandas showed that higher copulation and birth rates can be obtained if females are allowed to mate with their preferred partners, implying that mate choice may occur in wild pandas (Martin-Wintle et al., 2015). In addition, as a polygynous species, female giant panda lives by herself during the gestation, parturition and infant rearing period. Therefore, based on the parental investment theory (Trivers, 1972), females tend to be more active in mate choice to ensure maximum reproductive benefits.
The genetic diversities of MHC genes as adaptive genetic markers have been assessed among giant panda populations (Zeng, Yu, Pan, Wan, & Fang, 2005;Zhu, Ruan, Ge, Wan, & Fang, 2007;Zhu, Sun, et al., 2013;Zhu, Wan, Yu, Ge, & Fang, 2013). In particular, Zhang, Wu, Hu, Wu, and Wei (2015) analyzed the correlation between MHC genetic diversity and parasite infection in wild giant pandas and found a target MHC gene associated with parasite infection. However, the role of MHC genes in mate choice of giant pandas remains unexplored. In our study, we implemented a combination of long-term field observations of reproductive behavior F I G U R E 1 The Foping and Changqing Nature Reserves and the locations of 19 mating sites found from 2008 to 2016 and noninvasive genetic sampling (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) in the Foping and Changqing Nature Reserves to test whether MHC-based mechanisms play a role in male-male competition and mate choice in wild giant pandas.

| Study area, behavior observation, and sample collection
We implemented this study in the Foping nature reserve (Foping) and the adjacent Changqing nature reserve (Changqing) in Shaanxi Province (Figure 1) than other nature reserves, we chose them to ensure that we found as many mating sites (i.e., the aggregating sites of estrous males and female) as possible.
Two methods, tracking GPS-collared pandas  and following the unique vocalizations of courting or mating pandas, were used to search for mating sites from 2008 to 2016 in Foping and from 2008 to 2010 in Changqing (see details in Hu et al., 2017).
We identified the estrous female, dominant male (i.e., the male who won the combat), and the subordinate males (i.e., other males who participated in combat) by observing behavioral and spatial position cues (Pan et al., 2014;Schaller et al., 1985). After the reproductive activity was over, we collected fresh feces and shed hair samples using sterile polyethylene gloves following previous methods (Hu, Zhan, Qi, & Wei, 2010;Zhan et al., 2006). We also collected samples of mother-cub pairs through opportunistic encounters in the field, the tracking of collared females, and the monitoring of natal dens in winter. Feces samples of mother-cub pairs were attributed to either mothers or cubs based on fecal size. If the cub was not old enough to defecate, we collected the cub's hair sample when conditions allowed (see details in Hu et al., 2017). It was not possible to record data blind because our study involved focal animals in the field.

| Identification and sex determination of individuals
We extracted total genomic DNA from feces using the QIAamp DNA Stool Mini kit (QIAGEN, Hilden, Germany) and from hair using proteinase k in a PCR-compatible buffer (Allen et al., 1998). Blank controls were performed for both extractions and downstream amplifications. We used 14 giant panda-specific microsatellite loci to genotype the DNA samples (see details in Hu et al., 2017). To obtain reliable genotypes, we used a multi-tube amplification approach (Taberlet et al., 1996). First, we amplified each DNA sample three times, and if the genotype could not be determined, we performed two additional amplifications. PCR volumes and conditions as well as microsatellite genotyping are detailed in Hu et al. (2017). After obtaining multi-locus combined genotypes, we identified individuals following Zhan et al. (2006). We used MICRO-CHECKER (Van Oosterhout, Hutchinson, Wills, & Shipley, 2004) to detect the presence of genotyping errors, such as null alleles, large allele dropouts, and stuttering. No genotyping errors were detected in our final dataset. We used a Y-linked sexing marker (ZX1, 210-bp) in combination with an X/Y-linked amplification control (ZFX/ZFY, 130-bp) to determine the sex of each sample (Hu et al., 2017).

| MHC gene genotyping
The second exons of the three MHC class II loci (DRB1, DQA1, and DQA2) were selected for mate choice study of giant pandas because of three reasons: Firstly, MHC class II loci are the most reported loci involved in immune response, especially against extra-cellular parasites (Sommer, 2005); secondly, based on previous studies, the second exons of the three MHC class II loci (DRB1, DQA1, and DQA2) can be correctly amplified without the risk of falsely amplifying other paralogous loci (Chen et al., 2010;Zhang et al., 2015;Zhu et al., 2007), which is very important for reliable biological analysis; and thirdly, these loci were shown to have normal expression levels (Chen et al., 2010). Exons were amplified from DNA extracted from feces or hair. At the MHC class II DRB1 locus, all samples were genotyped using the primer pairs To ascribe the DQA alleles to a specific locus, for the MHC class II DQA1 locus, samples were genotyped using the primer pairs DQA1UP3B (5′-GTTTAGTAATCATCTTTTCTCCC-3′) and DQA1DN2 (5′-AGAGGCAGAGCATTGGACACATAC-3′; Chen et al., 2010). We modified two nucleotide bases in the forward primer DQA1UP3B according to the giant panda genome sequence downloaded from Ensembl (http://www.ensembl.org/info/data/ftp/index.html). For the MHC class II DQA2 locus, samples were genotyped using the primer pairs DQA2UP1 (5′-GTTTCTTCCGTCACTTGGCTTAATAAG G-3′) and DQA2DN1 (5′-AGGCAGAGCATTGGACACATACCAT-3′; Chen et al., 2010). All amplified alleles were reconstructed using the PHASE function (Stephens & Donnelly, 2003) in DnaSP 5.10 (Librado & Rozas, 2009) with the "recombination" model (−MR0) and 1,000 iterations after 100 burn-ins (Berggren & Seddon, 2008;Galaverni, Caniglia, Fabbri, Lapalombella, & Randi, 2013). The DQA alleles from the clone sequencing were used as the background for the PHASE analysis. Alleles were matched to sequences available in GenBank via Blastn (Johnson et al., 2008).
When the clone sequencing and PHASE data were combined, most of the individuals from the mating sites and mother-father-cub trios were successfully genotyped for the DQA1 and DQA2 loci.
Individuals with missing genotypes were excluded in the following data analysis.

| Testing the MHC-based mate choice hypotheses
We tested four mate choice hypotheses based on the MHC data from wild giant pandas, with different estimates of MHC variation for each hypothesis. To determine the MHC-based heterosis mechanism, we calculated the distribution of homozygotes and heterozygotes at each locus for dominant and subordinate males, respectively. To determine whether MHC genetic diversity affects panda mate choice, two measures were used: (a) the difference in coding amino acids between two haplotypes of an individual (Setchell, Charpentier, Abbott, Wickings, & Knapp, 2010); the difference in coding amino acids within the antigen-binding sites (ABS) between two haplotypes of an individual (Setchell et al., 2010). To investigate whether the MHC genetic compatibility hypothesis plays a role in mate choice, three methods were used to calculate the MHC genetic dissimilarity between males and females: (a) the mean number of differences in coding amino acids between female and male alleles using the formula (P CD + P Cd + P cD + P cd )/4 (Landry, Garant, Duchesne, & Bernatchez, 2001), where P CD , P Cd , P cD , and P cd are the number of amino acid differences between two individuals whose genotypes are Cc and Dd; (b) the mean number of differences in coding amino acids within the ABS between female and male alleles, with the same formula as the first method; and (c) twice the number of the proteins female and male individuals share divided by the sum of proteins of both individuals using D = 2Fab/(Fa + Fb), where Fa and Fb are the number of proteins in individuals A and B, and Fab is the number of proteins shared by both individuals (Wetton, Carter, Parkin, & Walters, 1987). At last, to determine whether the MHCbased "good gene" hypothesis was supported, the allele distribution of each locus between dominant and subordinate males was compared.

| Correlation test
To detect the impacts of genomic background on MHC variation, we tested the Pearson correlation among female-male relatedness based on microsatellite loci and female-male genetic compatibility parameters based on MHC variation, and tested the correlation among male's genetic heterozygosity based on microsatellite loci and male's MHC genetic diversity parameters, using psych package in RSTUDIO. Based on these correlation test results, we removed one of pairwise significantly correlated parameters for 10 female-male genetic similarity measures and 8 male's genetic diversity measures separately, and retained no significantly correlated parameters for further GLMM analysis.

| GLMM analysis
We performed generalized linear mixed model (GLMM) analysis (McCullagh & Nelder, 1989) to assess the roles of MHC heterosis, diversity and compatibility of three MHC class II loci in matingpair formation. We also included female-male relatedness and two measures for male's genetic heterozygosity (standardized individual heterozygosity (SH), Coltman, Pilkington, Smith, & Pemberton, 1999; internal relatedness (IR), Amos et al., 2001) based on microsatellite data as the explanatory variables (Hu et al., 2017). The response variable (mating-pair formation) was defined as binary ("1" denotes the pairing of dominant male and estrous female, and "0" represents the pairing of subordinate male and estrous female), given random effects for different sampling years, mating sites, individual identity or their combinations. The model was fit using the lme4 package in RSTUDIO, with a logit-link and binomial error distribution.

| Chi-square test
We used the chi-square test to investigate the good gene hypothesis, and the significance of differences in allele distribution between dominant and subordinate males was tested.

| Monte Carlo randomization tests
To detect whether females or males choose partners in a significant, nonrandom way, we compared the partners' mean MHC genetic dissimilarity with the individuals' randomized mean obtained by Monte Carlo simulations (Sin et al., 2015). We first calculated the mean MHC genetic dissimilarity for real parent pairs, that is, the mother and father of sampled cubs, then calculated the mean MHC genetic dissimilarity by randomly shuffling all female and male individuals (excluding cubs) 10,000 times. We calculated all sites of coding amino acid sequences and the ABS, respectively.

| Genotyping of the DRB1, DQA1, and DQA2 loci
For the DRB1 locus, we found five alleles in total, all of which had been reported in previous genetic diversity studies (Zhang et al., 2015; Figure 2a). For the DQA1 locus, we found eight alleles in total, including six alleles that had been reported previously  and two novel alleles that we deposited in GenBank (Aime-DQA*17-18, accession numbers: MG050735 and MG050736; Figure 2b). For the DQA2 locus, we found 3 alleles that had been reported previously ; Figure 2c). The corresponding NCBI allele names are listed in Supporting Information Table S1.

| Correlation estimate for MHC-based and microsatellite-based variables
Several individual samples failed to be amplified successfully due to low quality/quantity of fecal and hair DNA. As a result, for the DQA1 loci, 15 of 19 mating sites were included for the data analysis, while for the DQA2 and DRB1, 16 mating sites were included.  Table S2), suggesting that MHC genetic diversity is not dependent on genome-scale background diversity.
F I G U R E 2 Coding amino acid sequence alignment of the second exon of (a) 5 DRB1 alleles, (b) 8 DQA1 alleles, and (c) 3 DQA2 alleles. Note: "*" above sequences represents putative antigenbinding sites (ABS) according to Reche and Reinherz (2003). Among them, DQA1*07 and DQA1*08 are two novel alleles found in our study  Information Table S3). Here, compatibility1 denotes the method from Landry et al. (2001), and compatibility2 denotes the method from Wetton et al. (1987). The results showed that for each of three MHC class II loci, there were significant correlations among compatibility1, ABS_compatibility1, and com-patibility2. Among the nine genetic compatibility measures, only DQA1_compatibility2 and DRB1_compatibility2 had significant middle correlation coefficients with female-male relatedness. These results suggest that on the whole, MHC genetic compatibility is not dependent on genome-scale background relatedness.

| Testing mate choice hypotheses using mating site data
Based on the correlation test results, we selected one of significant correlated parameters for further GLMM analysis. As a result, we included as DQA1_heterosis, DQA1_diversity, DQA1_compatibil-ity1, DQA2_compatibility1, DRB1_heterosis, DRB1_diversity, and DRB1_compatibility1 as the explanatory variables, together with microsatellite-based SH and relatedness. The GLMM analysis (13 mating sites, n = 38) showed that neither of these genetic variation measures was significantly associated with mate-pair formation (Table 1). In addition, to increase sample size and detection power, we performed GLMM analysis for DQA1 (15 mating sites, n = 43) and DRB1 loci (16 mating sites, n = 43) separately (Supporting Information Tables S4 and S5). However, the two analysis results also detected no significant association between mate-pair formation and MHC variation measures. These findings suggest that matepair formation could not be explained by MHC heterosis, genetic diversity, and genetic compatibility hypotheses, acknowledging the caveat of the small sample size.

| Testing mate choice hypothesis using parentpair data
Of 11 parental pairs, nine pairs were successfully genotyped for the MHC class II DRB1 locus. Due to the low quality/quantity of fecal or hair DNA and the difficulty in amplifying the DQA1 and DQA2 loci, the majority of parent-pair samples were not successfully genotyped at the DQA1 and DQA2 loci. Thus, we did not perform parent pair data analysis for the DQA1 and DQA2 loci.
We calculated the MHC genetic compatibility parameters at the DRB1 locus. Using a Monte Carlo randomization test, we found that MHC genetic compatibility measures did not differ significantly between the real parent pairs (n = 9) and random female-male pairs (n = 40) for the coding amino acid sites (3.9722 vs. 4.125, p = 0.89) and for the coding amino acids within the ABS (2.52 vs. 2.8126, p = 0.62). Indeed, the sample size was very small, which may affect the reliability of our conclusion.

| D ISCUSS I ON
We implemented a long-term ecological, behavioral and genetic study of wild giant pandas in the Foping and Changqing Nature Reserves by tracking GPS-collared giant pandas (Hu et al., 2017;Nie et al., 2012;Wei et al., 2014;Zhang et al., 2014). One aim of the long-term field project is to explore the genetic mechanisms to 2016, we collected behavioral and genetic data from 19 mating sites and 13 mother-cub pairs. Despite our long-term attempts to find additional mating sites and mother-cub pairs, the sample sizes were yet small, which is common in endangered species studies.
Based on these mating sites and mother-cub data, Hu et al. (2017) assessed the inbreeding level and inbreeding avoidance mechanisms of giant pandas and found that female-male relatedness and male  (Nie et al., 2012;Pan et al., 2014;Schaller et al., 1985). An estrous female is fertile for approximately 3 days every mating season. Within one mating site, an estrous female is often detected by 1-6 males, which greatly limits access of males. These circumstances are not ideal for male mate choice. Males must invest considerable time and energy in the pursuit of and competition for females, and the number of available estrous females is small.
Considering these constraints, giant pandas may use a "better than nothing" mating strategy, that is, a less optimal mate is better than no mate. As Hu et al. (2017) found, giant pandas also did not use female-male relatedness or male genetic heterozygosity to locate a mate. In a similar manner, a study on snub-nosed monkeys failed to find evidence supporting MHC-based mate choice mechanisms.  (Brennan, 2004). It is important that body size may impact mate choice. In the recent past, Martin-Wintle et al. (2015) found that in captive giant pandas, male body mass significantly affected cub production in both female and male panda reproductive performance experiments and significantly affected copulation success in females. Nie et al. (2012) also found that body size may be the primary factor in male dominance status by visually ranking the body size of males involved in the same mating site.
These findings are consistent with our field behavior observations that the dominant males often appeared to be stronger than the subordinate males.
At last, our failure to find evidence of mate choice may have resulted from a small sample size and the low diversity in the selected MHC loci. As discussed in Hu et al. (2017), although we performed long-term fieldwork, the number of mating sites and mother-cub pairs were yet small, which may affect our conclusions. In addition, we found 5 DRB1 alleles, 8 DQA1 alleles and 3 DQA2 alleles Nature Reserves for assistance with sample collection, and Lei Zhang for help with experiments.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
FWW and YBH designed the study. YBH and YGN performed the field study and collected the genetic samples. YBH, LJY, and LY conducted the genetic analysis and data analysis. LJY, YBH, and FWW wrote the manuscript.

DATA ACCE SS I B I LIT Y
The two MHC DQA1 gene sequences new identified are deposited in GenBank (accessions numbers MG050735 and MG050736).