Disentangling the mechanisms of mate choice in a captive koala population

Study samples

A total of 70 koala DNA samples were extracted from banked whole blood and tissue samples (previously collected under San Diego Zoo Global IACUC protocols 10-008, 10-009, 11-029, 14-034) using either standard phenol-chloroform extraction, or the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) following manufacturer’s directions. Detailed pairing records, and studbook data, were provided by the Association of Zoos and Aquariums North American Regional Studbook Keeper (Chris Hamlin, personal communication, April 2017). These pairing records spanned 1984–2012 and contained mate choice data and breeding outcomes for every pairing (n = 964) at the zoo throughout this period. Breeding recommendations are reviewed annually, and predominately follow a minimize kinship strategy (Ballou & Lacy, 1995). During pairing, oestrous females are placed with males in enclosed cubicles for 5–10 min and mating behavior is monitored throughout this time (Bercovitch et al., 2006). If the pair does not copulate, the female and male may be paired with other conspecifics, or the same pairing may be trialed again at a later time. This pairing process means both female and male koalas are often exposed to multiple individuals of the opposite sex within and between seasons. We are therefore able to determine whether mate choice is occurring at the pre-copulatory level (i.e., copulation success as a proportion of pairing attempts) or post-copulatory level (i.e., breeding success (production of offspring) as a proportion of copulations), and whether it affects offspring success (offspring that survive more than 1 year as a proportion of offspring produced).

DNA samples were available for individuals across 28 years of pairing data (Fig. S1), accounting for 51% of the individuals and 49% of the pairings in the complete dataset (Table S1). A total of 22 sampled koalas were never paired for breeding and were excluded from the final mate choice analyses. These koalas were genotyped during the study and included in marker diversity statistics (see “Marker diversity” in the Supplemental Information) to maximize sample sizes where possible. All non-genetic factors were analyzed using the complete pairing dataset (including unsampled koalas).

Non-genetic determinants of koala mate choice

We used generalized linear models (GLMs) in R v 3.4.0 (R Core Team, 2017) to test for non-genetic effects, including year, age and familiarity, on mating success (copulation, breeding and offspring success; see below). GLMs were performed with binomial distribution as follows: for each pairing event (n = 964) we modelled whether the pair successfully copulated (1) or did not copulate (0), with predictor variables including year of the pairing (to account for changes in the breeding program over time), age of the female, age of the male and the number of years the male and female had previously been paired together (as a measure of familiarity). Age² was also included, as the relationship between age and mating success was not predicted to be linear (Rose, 1991). The dataset was then subset into only those pairs that successfully copulated (n = 304) and the same predictor variables were modelled against whether each of these pairs successfully bred (produced offspring) (1) or did not breed (0). The dataset was then further subdivided into only those pairs that successfully bred (n = 134) and the same predictor variables were modelled against whether those pairs produced offspring that survived more than 1 year (1) or did not survive more than 1 year (0).

Male body mass is a strong predictor of male koala breeding success in the wild (William & Bercovitch, 2011) and correlates strongly with male age (Tobey et al., 2006). We did not have male body mass available for our study animals. As age difference between the male and female was also highly correlated with both female and male age (ρ = −0.58, p = < 0.001 and ρ = 0.76, p = < 0.001, respectively, Table S2) it was not included in the model. As a result, any effects of male age that we observe may also be reflecting the effect of male body mass and/or the age difference of the pair; female age effects may also be driven by age difference.

Although some pairs were repeated in multiple years, Pair ID was not included as a random factor due to the majority of pairs (60%) only being represented in 1 year of the dataset (Fig. S2) (models with Pair ID fitted as a random intercept did not converge). Variance inflation factors (VIFs; Belsley, Kuh & Welsch, 1980) were calculated for the remaining predictor variables to ensure there were no adverse effects of multicollinearity. All VIFs were <2 and so year, female age, male age and familiarity were included in the models (Belsley, Kuh & Welsch, 1980). Model predictors were standardized by subtracting the mean and dividing by two standard deviations (following Gelman, 2008) to facilitate inference of regression coefficients within and between models (Schielzeth, 2010). Model fitted values were back-transformed onto the natural scale for plotting and interpretation.

MHC genotyping

Major histocompatibility complex-linked primers were designed within 10 kb of MHC Class II genes that had been annotated by the Koala Genome Consortium (Johnson et al., 2018). RepeatMasker (Smit, Hubley & Green, 2013–2015) was used to identify microsatellite sequences <10 kb away from the MHC genes (Cheng et al., 2009b). Candidate microsatellite sequences (PhciDBB001M3, PhciDCBM1 and MHCIIDAB001M1) were selected based on minimal interruptions to the repeat sequence and low proximity to other repeat regions. These microsatellites were linked to genes of the DB, DC and DA families respectively (Abts, Ivy & DeWoody, 2015; Lau et al., 2013, 2014a; Lau, Griffith & Higgins, 2014b), allowing us to incorporate a representative for each classical marsupial MHC class II family (Belov et al., 2006; Belov, Lam & Colgan, 2004). The repeat motifs for each microsatellite were (TG)₁₃, (GA)₂₈ and (AC)₂₉ respectively. We extracted these microsatellite sequences with 300 bp of flanking sequence and designed PCR primers using Oligo 7 (Rychlik, 2007). Primer sequences were then used in a BLAST search (Priyam et al., 2015) against the koala genome (Johnson et al., 2018) to ensure specificity and prevent amplification of non-target sequences. Primer sequences used to amplify the three microsatellites were as follows (5′ to 3′): PhciDBB001M3 F:TTCTCTTGTCCTTCTTGTGTC, R:TTCTCCCTACAAAGATGATCC; PhciDCBM1 F:AGTCTGGTGTCATTAGCAATAGG, R:CTGAATGAGGCAAGGGAGAG; MHCIIDAB001M1 F:ACACTACTTCCCTGAATCTGAC, R:TACAGTGTTACTTCATGCAGAG.

All loci were initially screened for polymorphism using DNA samples from koalas previously found to be polymorphic at MHC loci (Cheng et al., 2017) (see “Initial primer screening and optimization methods” in the Supplemental Information) before genotyping the study population at these markers. Since all three markers had similar product lengths, they were not multiplexed for typing. PCRs were carried out for each locus using the Type-it Microsatellite PCR Kit (Qiagen, Hilden, Germany) with a modified total reaction size of 10 μL and the following modified primer concentrations: 0.06 μM tagged primer, 0.6 μM untagged primer and 0.6 μM 6-FAM labelled CAG tag (Schable, Fischer & Glenn, 2002) (see Table S3 for thermocycling conditions). Capillary electrophoresis was undertaken at the Australian Genome Research Facility using MCLAB DSMO-100 Orange Size Standard. Alleles were manually called using GeneMarker (Hulce et al., 2011). Controls included a negative control using water, and a positive control using DNA from a koala that was successfully genotyped during the initial primer development.

Genotyping of non-MHC-linked microsatellites

Koalas were genotyped at a further 15 microsatellite loci, not known to be linked to MHC, using primers from previous studies (Cristescu et al., 2009; Dennison et al., 2017; Houlden, England & Sherwin, 1996) (Table S4). The genomic location of each microsatellite was confirmed using the NCBI koala assembly browser (Kitts et al., 2015). All microsatellites were located on scaffolds not containing any MHC genes, and seven were >10 kb away from any genes. Seven markers (Pcin05, Pcin08, Pcin11, Pcin20, Pcin21, Pcin22 and Pcin23) were split into three multiplexes and typed using a fluorescently labelled (6-FAM) CAG-tag (Schable, Fischer & Glenn, 2002) (Table S4), with PCRs carried out using the Type-it Microsatellite PCR Kit (Qiagen, Hilden, Germany) (as above). The remaining eight markers were amplified separately using a fluorescently labelled (HEX or 6-FAM) M13 tail (0.6 μM) (Schuelke, 2000) or forward primer (Table S4). PCRs were carried out in a 10 μL reaction volume containing 1× PCR buffer, 2.5 mM MgCl, 0.2 mM dNTPs (Apex; Genesee Scientific, San Deigo, CA, USA), 0.2 μM forward primer, 0.6 μM reverse primer and 0.5 U AmpliTaq Gold (Applied Biosystems, Waltham, MA, USA) (see Table S3 for thermocycling conditions). Samples were genotyped on an ABI 3130xl using an internal GeneScan 500 ROX size standard and alleles were automatically called then manually checked using GeneMapper (ABI).

Microsatellite diversity

Approximately 20% of the koalas were re-genotyped to determine genotyping error rate. Tests for evidence of null alleles, deviation from Hardy–Weinberg equilibrium, and linkage disequilibrium were performed to ensure all of the non-MHC and MHC markers were suitable for use in the final statistical analysis (Methods and Results for these analyses are provided in “Marker diversity” in the Supplemental Information and Table S5). Standardized heterozygosity (H_s) was calculated as a measure of individual multilocus heterozygosity at both the non-MHC and MHC markers using the Rhh package (Alho, Välimäki & Merilä, 2010) in R. We chose this method, as H_s gives equal weighting to all loci examined despite variation in the number and frequency of alleles present across the markers used (Aparicio, Ortego & Cordero, 2006; Coltman et al., 1999). A Spearman’s rank correlation between standardized heterozygosity at MHC-linked markers and standardized heterozygosity at non-MHC markers was also performed to test whether MHC diversity and genome-wide diversity were correlated. This was necessary to determine whether any MHC-associated mate choice findings were by-products of genome-wide variation (Ferrandiz-Rovira et al., 2016).

Statistical analysis

GLMs were used to test the three mate choice hypotheses at pre-copulatory, post-copulatory and offspring survival stages using the three binomial response variables: copulation success, breeding success and offspring success respectively (defined above). In each model, the predictor variables (male heterozygosity, pair similarity and allele presence/absence) were as described in the sections that follow. Multicollinearity of predictor variables in all models was checked by calculating VIFs, and model predictors were standardized to facilitate inference across predictors, as described above.

Non-genetic determinants of koala mate choice

Male age had a significant effect on copulation success (Table 2). Expected copulation success rates increased from ∼20% when males were 2-years-old, to 40% when males reached 12 years of age, and decreased to below 35% when males were 17 years or older (fitted values are taken from the regression model in Table 2). Copulation success increased significantly with increasing familiarity between pairs (Table 2). Dyads that had previously been paired together for 5 years or more had expected copulation success rates above 50% (95% CI [0.36–0.65]), compared to the 34% success rate of dyads that had never been paired (95% CI [0.28–0.41]; fitted values are taken from the regression model in Table 2). No association was found between female age and mating success (Table 2). None of the factors we tested influenced breeding or offspring success and could not explain the strong declines in copulation and breeding success across time (Table 2).

Microsatellite genotyping

All koalas were genotyped at the three MHC markers and >75% of the study population was genotyped at 13 or more of the non-MHC loci (Table S5). Genotyping error rate was very low (0.53%). For the MHC-linked markers, the mean number of alleles per locus (N_a) was 9.67 (range 7–13; n = 3 loci, 70 koalas), mean observed heterozygosity (H_o) was 0.771 ± 0.1 (SE), and mean expected heterozygosity (H_e) was 0.755 ± 0.083 (Table S5). Standardized heterozygosity (H_s) for the MHC markers ranged from 0.432 (more homozygous) to 1.3 (more heterozygous) (Table S5). Non-MHC marker H_s ranged from 0.347 (more homozygous) to 1.62 (more heterozygous). For the non-MHC markers, N_a = 6.53 (range 2–10; n = 15 loci, 70 koalas), H_o = 0.673 ± 0.05 and H_e = 0.662 ± 0.046 (Table S5). Standardized heterozygosity at MHC-linked loci was not correlated with standardized heterozygosity at non-MHC loci (n = 70 koalas, ρ = −0.06, 95% CI [−0.30–0.19], p = 0.614).

Genetic determinants of koala mate choice

Non-genetic determinants of koala mate choice

Our analysis showed that koala copulation success is significantly influenced by male age and/or the age difference between males and females, in line with previous studies (Bercovitch et al., 2006). Studies in other species suggest that females may prefer to mate with older males, likely due to older males being of a higher genetic quality through viability selection (Manning, 1985; Trivers, 1972). In koalas, male size, bellowing and sternal scent secretions have been found to convey age-related information (Charlton et al., 2012; Salamon & Davies, 1998; Tobey et al., 2006), and so females may use visual, auditory and chemical cues to select mates based on age (Bercovitch et al., 2006). Here, we provide additional evidence that male age influences captive koala mate choice, although the precise chemical and auditory mechanisms by which females receive and utilize this information remain unclear (Ellis et al., 2015; Tobey, Nute & Bercovitch, 2009).

In addition to age, we also found that familiarity may promote copulation success in captive koalas. Other mammal mate choice studies show that females have a preference for more familiar males (Roberts & Gosling, 2004) and mating with familiar males leads to increased reproductive success (Martin & Shepherdson, 2012). This familiar male preference often arises in territorial scent-marking species, as females encounter scent marks of locally territorial males and select these males due to their ability to defend a territory (Rich & Hurst, 1998). Although koalas are a territorial scent-marking species (Allen et al., 2010), female koalas do not show a preference for locally territorial males in the wild (Ellis, Hale & Carrick, 2002). The familiarity trend in our study may be driven by pairing previously successful pairs together in subsequent years, although most (60%) of the pairings in our dataset were from first-time pairings. Further research, which directly examines the role of familiarity, is needed to confirm its influence in koala mate choice.

Genetic determinants of koala mate choice

Previous research suggests that females of many species are often more attracted to heterozygous males, and heterozygosity has been linked to numerous advantages such as greater sexual ornamentation, mating success and overall reproductive success (reviewed in Kempenaers, 2007). Despite these advantages, genome-wide heterozygosity did not influence mating success in our captive koalas. Some species also display a preference for dissimilar individuals, which may reduce inbreeding and increase genetic diversity of offspring (Ferrandiz-Rovira et al., 2016; Kempenaers, 2007). In contrast, we found a positive association between genome-wide similarity and breeding success, suggesting female koalas are more likely to produce offspring with males that are more genetically similar overall. A recent review showed mating with similar individuals can allow populations to adapt more quickly to virulent diseases and parasites (Campbell et al., 2017). Assortative mate preferences may help protect koala populations from threatening infectious diseases such as chlamydia (Polkinghorne, Hanger & Timms, 2013) and koala retrovirus (Denner & Young, 2013), and should be examined further. Although similar numbers of neutral microsatellite markers have been used in recent studies to examine genome-wide mating preferences (Ferrandiz-Rovira et al., 2016; Huchard et al., 2013), estimates of genome-wide diversity based on 15 microsatellites may not be sufficient, and larger numbers of markers should be employed in future studies to provide more accurate measures of genome-wide diversity (Miller et al., 2014).

In addition to genome-wide mating preferences, many species select mates based on MHC (Table 1). Consistent with these studies, we found that koala mating success showed a significant association with male heterozygosity and pair similarity at MHC loci, as well as the presence or absence of particular MHC alleles. In contrast to the quantity of alleles hypothesis, males that were less heterozygous at MHC-linked loci showed a greater rate of copulation success. This indicates that female koalas prefer to copulate with males that have fewer alleles at MHC loci, particularly at DAB loci (we note that individual locus heterozygosity and overall MHC heterozygosity are related and do not provide multiple lines of evidence). Interestingly, among those males that did copulate, breeding success was higher for more heterozygous males. This implies that females are more likely to produce offspring when breeding with males of higher heterozygosity, than when breeding with males of lower heterozygosity. The standardized slopes of the trends at each stage were of similar magnitude (Table 3), suggesting that the effect of heterozygosity is similar at both mate choice stages. Assessed together, these results reflect differences in the pre-copulatory and post-copulatory MHC-dependent choice mechanisms in koalas.

Vertebrate females are known to select sperm based on heterozygosity or diversity at MHC loci (Wedekind, 2002; Winternitz et al., 2013). Males heterozygous at MHC loci show significantly greater fertilization success relative to homozygous males (Skarstein et al., 2005). In koalas, males with low heterozygosity at MHC loci overall (particularly at DAB loci) have a higher probability of copulating; however, more heterozygous males experience a fertilization advantage, so that their copulations are more likely to result in the production of offspring. While the benefits of breeding with heterozygous males can be explained by the increased antigenic peptide repertoire and immunocompetence of heterozygotes (Kamiya et al., 2014; Landry et al., 2001), we are unaware of any other reports where less-heterozygous males have a higher copulation success rate. Future work should investigate whether this unexpected relationship may be driven by an unmeasured male trait that is correlated with MHC heterozygosity.

Contrary to many previous findings under the genetic compatibility hypothesis (Table 1), captive koala pairs that were more similar at the MHCII DAB-linked locus had greater breeding success than less similar pairs. Female koalas were more likely to produce offspring with males that shared alleles at DAB loci, which is consistent with a greater reliance on post-copulatory MHC-dependent mechanisms of mate choice. Numerous studies have shown that females select sperm based on the genetic dissimilarity of mates (Olsson et al., 1996; Thuman & Griffith, 2005), particularly at MHC loci (Schwensow, Eberle & Sommer, 2008; Yeates et al., 2009). A preference for mates that are more similar at MHC loci may evolve in response to disadvantages associated with mating with individuals that are too dissimilar, including increased risk of autoimmune disorders due to suboptimal T-cell selection (Kaufman, Völk & Wallny, 1995), reduced recognition of foreign peptides due to T-cell loss (Vidovi & Matzinger, 1988), and disruption of co-adapted gene complexes (Hendry et al., 2000).

Studies have shown that, in some circumstances, carrying multiple copies of the same MHC allele allows for higher disease resistance (Grimholt et al., 2003; Nuismer, Otto & Blanquart, 2008). However, MHC assortative mating may make populations more vulnerable to future disease outbreaks or other stochastic events (Campbell et al., 2017). We suggest that female koalas might not solely choose more-similar mates but may rather optimize the quantity and combination of MHC alleles in the offspring (see also Milinski, 2006). A similarly complex mate choice mechanism has been demonstrated in sticklebacks (Gasterosteus aculeatus), whereby females prefer to mate with males with genotypes that, when combined with their own MHC alleles, will produce offspring with an optimal number of alleles and provide the highest possible resistance against common parasites (Milinski, 2003). A study of the numbers and combinations of MHC alleles in offspring with reference to the parents’ MHC is needed to confirm whether this mechanism exists in koalas.

In line with the advantage of particular alleles hypothesis, we found that the presence of certain MHC-linked microsatellite alleles was associated with increased or decreased mating success, whilst other alleles showed no association. Since the MHC-linked microsatellites are found in non-coding regions of the genome (and likely have no direct functional consequences), this finding suggests that females are selecting for and/or against males that carry the respective MHC alleles. Although some of these models appear to provide strong evidence for the influence of certain alleles on mating success, we note cautious interpretation of these findings is warranted as small sample sizes and subject: predictor ratios limit the reliability of the models (Table 5). Even so, by using a multimodel inference framework, i.e., ranking models under information theory, we provide support for the hypothesis that some MHC-linked microsatellite alleles are associated with mating success. This is probably because they are in linkage disequilibrium with important functional MHC sequence variants that potentially play a role in disease resistance and susceptibility to common infectious diseases such as Chlamydia (Lau, Canfield & Higgins, 2012).

Mate choice influences offspring viability in a variety of species, particularly when mating preferences are MHC-dependent (Agbali et al., 2010; Von Schantz et al., 1996). In koalas there was no association between offspring survival and mate choice preferences, however, data on offspring survival was only present for koala pairings that produced offspring (n = 13 males, 26 pairs, 28 offspring in total). Investigation into early joey loss will ascertain if observed mating preferences produce offspring with optimal MHC, as successful matings resulting in offspring with optimal immunity will increase offspring viability whilst suboptimal MHC will, in theory, produce early offspring deaths.

By using MHC-linked microsatellites we were able to examine three families of MHC loci to look for general patterns as well as locus-specific patterns. The majority of MHC-dependent mate choice studies in the current literature often only examine a single locus due to the limitations of MHC typing techniques (Kamiya et al., 2014). Our data indicate that some loci may play a larger role in mate choice than others, and different loci may act in different ways, further indicating the importance of examining multiple MHC loci. Huchard et al. (2013) found that female grey mouse lemurs (Microcebus murinus) chose males based on a particular MHCII locus that was under stronger diversifying selection than other MHCII loci. Similarly, DAB loci in the koala have previously been found to be under stronger selection than other MHCII loci (Lau et al., 2013), which may explain the strong effect of DAB on koala mate choice in the current study. Multi-locus approaches are vital in gaining a holistic understanding of MHC-dependent mate choice mechanisms (Kamiya et al., 2014) and can be easily achieved using MHC-linked microsatellite markers. Although numerous studies have confirmed MHC-linked microsatellite markers as a good proxy for MHC diversity in other species (Cheng & Belov, 2012; Cheng et al., 2009a; Crouau-Roy et al., 1996), this association needs to be confirmed in the koala.