Genetic variability and ontogeny predict microbiome structure in a disease-challenged montane amphibian

Amphibian populations worldwide are at risk of extinction from infectious diseases, including chytridiomycosis caused by the fungal pathogen Batrachochytrium dendrobatidis (Bd). Amphibian cutaneous microbiomes interact with Bd and can confer protective benefits to the host. The composition of the microbiome itself is influenced by many environment- and host-related factors. However, little is known about the interacting effects of host population structure, genetic variation and developmental stage on microbiome composition and Bd prevalence across multiple sites. Here we explore these questions in Amietia hymenopus, a disease-affected frog in southern Africa. We use microsatellite genotyping and 16S amplicon sequencing to show that the microbiome associated with tadpole mouthparts is structured spatially, and is influenced by host genotype and developmental stage. We observed strong genetic structure in host populations based on rivers and geographic distances, but this did not correspond to spatial patterns in microbiome composition. These results indicate that demographic and host genetic factors affect microbiome composition within sites, but different factors are responsible for host population structure and microbiome structure at the between-site level. Our results help to elucidate complex within- and among- population drivers of microbiome structure in amphibian populations. That there is a genetic basis to microbiome composition in amphibians could help to inform amphibian conservation efforts against infectious diseases.


Introduction
Multicellular organisms act as hosts to a diverse suite of bacterial communities, collectively referred to as the microbiome. Recent research across multiple taxa has highlighted that these bacterial communities perform beneficial functions for the host, including protection from infectious pathogens [1][2][3]. Despite the importance of symbiotic relationships between microbes and their hosts, a comprehensive understanding of factors that influence the diversity and composition of the microbiome, particularly for non-human animals, is lacking. Though previous studies have revealed landscape-scale variation among populations in microbiome structure [4][5][6], and fine-scale variation among individuals within populations due to life stage [6], diet [7], gender [7] and genotype [7,8], few have investigated both within and among population drivers of microbiome structure in a single framework [6,8]. Studying factors that determine variation in the distribution of bacterial symbionts at both the individual and landscape scale is crucial for understanding how individual susceptibility to pathogens varies within and among populations.
The amphibian skin microbiome represents a model system for understanding the ecological drivers of host microbiome structure and the role of symbiotic bacteria in protecting the host from pathogens. The lethal and globally distributed amphibian chytrid fungus (Batrachochytrium dendrobatidis; Bd) is causing mass mortalities and rapid population declines of amphibians around the world and is a major driver of the current amphibian extinction crisis [9]. Certain host-associated bacteria have been shown to influence susceptibility to Bd, with particular bacterial genera conferring increased or decreased protection against Bd in experimental and field studies [2,[10][11][12]. Host genotype is one factor that may underpin variation in the presence and persistence of these protective microbes, both within and among host populations [7,8,[13][14][15][16]. For example, polymorphism of the major histocompatibility complex (MHC) Class IIb gene has been shown to co-vary with gut microbiome structure in vertebrates [7]. The expression of antimicrobial peptides (AMPs) on frog skin is governed by innate immunity [17], and so variation in immune genes may drive structural changes in the microbiome through differential AMP expression. The population genetic structure could therefore yield differences in microbiome at the landscape level if it is representative of the non-random distribution of functional genetic variation at immune loci. Furthermore, shifts in the environmental reservoir of bacteria across large geographic scales may also cause differences in the relative abundances of bacteria able to colonise the host [4]. Within populations, fine-scale variation in life-history traits, such as age, may further modify the ability of bacteria to coexist on the skin, or favour differential selection by the host from the environment. Though microbiome differences have been shown among larvae, juveniles and adults in some amphibian species [18,19], changes in microbiome structure throughout larval ontogeny have not been investigated. This may be especially pertinent, given that fine-scale variation in Bd infection loads has been found across developmental stages of tadpoles [20,21]. Research into the links between host genetics, developmental stage, microbiome composition and pathogen susceptibility may enable a better understanding of the factors governing disease in vulnerable populations.
Here we use a disease-challenged high-altitude amphibian species in the Drakensberg Mountains of southern Africa to investigate within and among population predictors of microbiome structure that may have implications for disease susceptibility. Using microsatellite genotyping and 16S amplicon sequencing, we examine how population genetic structure, individual genotype, and developmental stage determine differences in microbiome structure among individuals. Specifically, we aim to: (i) quantify Bd infection across Amietia hymenopus frog populations; (ii) characterise the microbiome of A. hymenopus; (iii) examine the genetic population structure of A. hymenopus on the Drakensberg plateau; (iv) investigate predictors of microbiome structure among sites (specifically river, geographic distance and genetic differentiation among host populations); (v) investigate within-site predictors of microbiome structure (specifically host genotype and ontogenetic stage); and (vi) examine the roles of host microbiome and host genetics in influencing Bd infection loads.

Study site
The Drakensberg Mountains of southern Africa are a vast mountain range reaching up to 3482 m above sea level (asl). Here, the Mont-aux-Sources plateau crosses north-east Lesotho and South Africa. On the South African side, four streams (Vemvane, Tugela, Bilanjil and Ribbon Falls) flow 1-3 km across the plateau from multiple origins in the high peaks on the South Africa-Lesotho border, before dropping off a~1000 m precipice into South Africa (Fig. 1). Two other streams on the plateau flow into Lesotho; Khubela and Kubedu. The Phofung river frog (Amietia hymenopus) is the only amphibian found in most of the Mont-aux-Sources stream system, with the exception of the Kubedu river, where A. hymenopus is absent and A. umbraculata is present [22]. Long-term (10 year) sampling in this region indicates A. hymenopus has a consistent history of Bd infection with associated mass mortalities [23].

Sample collection
We collected tadpoles from Mont-aux-Sources by permission of Ezemvelo KwaZulu Natal Wildlife and the Kingdom of Lesotho Ministry of Tourism, Environment and Culture. We sampled eight sites across five rivers over the plateau; 'Nampolice' and 'Vemvane' from the Vemvane River (South Africa); 'Tugela', 'Tugela 2' and 'Tukelahed' from the Tugela River (South Africa); 'Khubnam' from the Khubela River (Lesotho); and 'Bilanjil' and 'Ribbon Falls' from their respective rivers (South Africa; Fig. 1). Highresolution mapping was unavailable for the study area, so we used QGIS (version 2.18) to combine Open Source data sets and then digitise the water courses within the study area (Supplementary Information).
We collected 20-28 tadpoles per site (nine at Khubnam) along a~150 m transect within an altitudinal range of 3010-3090 m asl to minimize variation in Bd infection. We conducted two sampling trips 4-weeks apart; in March 2015, we collected tadpoles from six sites (Nampolice, Tugela, Tugela 2, Tukelahed, Bilajil and Ribbon Falls; Fig. 1) and analysed these for Bd infection, host genotype and microbiome composition. In April 2015, we collected tadpoles from Vemvane and Khubnam to enhance the population genetics analyses, but did not carry out microbiome analyses on these samples in order to avoid spurious variation caused by temporal changes (Table 1). We stored tadpoles in 95% ethanol prior to Gosner staging [24] and exported carcasses for molecular analyses in the UK by permission of the Department of Agriculture, Forestry and Fisheries, South Africa, and the Department for Environment, Food and Rural Affairs, UK. We excised full mouthparts for DNA extraction as due to their keratinised nature, Bd primarily affects this area in tadpoles. DNA was extracted using the DNeasy® Blood and Tissue Kit (Qiagen) and was used for all subsequent analyses.

Bd infection
We constructed serial dilutions of DNA and conducted quantitative PCRs on a Qiagen Rotor Gene according to Boyle et al. [25]. We quantified Bd infection intensity in triplicate using a 1:10 dilution of DNA with Bd standards ranging from 0.025 to 250 zoospores/µl. We calculated infection prevalence and intensity for each site.

16S amplicon sequencing of bacterial communities
We sequenced the microbial communities present on tadpole mouthparts using a modified version of the protocol in Kozich et al. and Bates et al. [26,27]. Transient bacteria most likely rinsed off during the ethanol preservation and so microbial analyses reflect the resident microbiome in the mouthparts of tadpoles. Briefly, we amplified a V4 region of the 16S rRNA gene in triplicate for 139 samples, as well as two negative controls (one per PCR plate). Triplicates were pooled and normalized using a titration run on an Illumina MiSeq using v2 nano chemistry, followed by a full run using 250 bp paired-end v2 chemistry. We also sequenced a mock community containing 20 bacteria isolated from amphibian skin (genera Acinetobacter, Citrobacter, Enterobacter, Flavobacterium, Plantibacter, Pseudomonas and Serratia).

Analyses
We performed sequence processing in DADA2 v1. 4 [31]. Following Longo and Zamudio [32], we excluded SVs that contained fewer than 100 reads across the entire dataset (n = 19,395), yielding 2384 SVs. Results using all SVs are similar and not shown here. After filtering, mean library size across individuals was 20,160 reads (range 13,156-83,063). We rarefied all libraries to 13,156 reads per sample for alpha and beta diversity analyses. We calculated alpha diversity metrics (Shannon index) using the 'estimate_richness' command in phyloseq. We visualised differences in Bray-Curtis distance among samples using non-metric multi-dimensional scaling (NMDS) ordination of bacterial communities using the R package vegan [33], setting k to 3 to yield a stress value of <0.17. We also used the vegan function 'vegdist' to calculate Bray-Curtis distance between bacterial communities to compare to genetic distances among samples, and the 'adonis' function to compare bacterial community structure among populations using permutational MANOVA. We used the R package gplots [34] to visualize the number of shared SVs among populations with a Venn diagram. To simplify the Venn diagram, we collapsed the two 'Tugela' sampling points into one group.

Analyses
Details of quality control and summary statistics carried out can be found in Supplementary Information. To assess levels of population differentiation, we calculated global F ST and pairwise F ST values between sites with corrections for null alleles using the ENA method of Chapuis and Estoup [37] with 1000 bootstrap replicates, and pairwise D [38] in Genodive v2.0b23. [39] We used ADZE v1.0 [40] to calculate site-level rarefied allelic richness (maximum g = 20) to study the effects of intraspecific genetic diversity at the population level on microbiome variation. We excluded Khubnam due to the low sample size (n = 9) compared with other sites (n = 20-28) to ensure sufficient power to detect differences among sites was retained. We created a pairwise matrix of Bruvo genetic distance [41] between all individuals in GenoDive. We also calculated individual-level heterozygosity (proportion of heterozygous loci in all loci successfully genotyped) using Genhet [42] in R to study its effect on associated microbial communities, as individual heterozygosity can be correlated with fitness.
We conducted principal coordinates analysis (PCoA) in GenAlEx v6.502 [49,50] using pairwise D values and the standardised covariance method. We ran an analysis of molecular variance (AMOVA) in Genodive using the Infinite Allele Model. We tested for the presence of isolation by distance among all locations and within population clusters. We conducted a Mantel test in ade4 in R [51] in R with 999 permutations to test associations between matrices of pairwise linearised F ST values (F ST /1-F ST ) and the logarithm of geographic distances. We used GENECLASS2 v2.0 [52] to identify putative first-generation migrants among sites [53,54].
We used INEST v2.1 [55] to identify the presence of two genetic signatures of bottleneck events; heterozygosity excesses in respect to allelic richness [56] and deficiency in M-Ratio values (mean ratio of the number of alleles to the range in allele size; [57]).

Individual level data
To examine the factors affecting Shannon diversity, we fitted a Gaussian Mixed effects model in the R package lme4 [58] with Shannon index as the response, Gosner stage and heterozygosity (specified as proportion of heterozygous loci, [59]) as predictors, and population ID as a random intercept (n = 126 data points for which there were both microbial (Shannon) and host genetic (proportion of heterozygous loci) data). Gosner stage and proportion of heterozygous loci were z-transformed prior to model fitting. The most complex model contained a proportion of heterozygous loci × Gosner stage interaction. We assessed significance of terms using likelihood ratio tests between models estimated using maximum likelihood. To quantify differences in mean richness among populations, we extracted the posterior modes of the population random effects from a Bayesian version of the best-supported model, specified in the R package MCMCglmm [60].
To examine the factors predicting microbial community structure (beta diversity), we calculated the withinpopulation 'divergence' metric of beta diversity using the microbiome package (Lahti et al. 2017) to be used as a response in a Gaussian mixed effects model. The community structure model contained Gosner stage, proportion of heterozygous loci and their interaction as fixed effects, and population ID as a random effect. We calculated the r 2 values of minimal models using the [61] method for mixed effects models, as implemented in the R package MuMIn [62]. We used a partial Mantel test to examine the correlation between genetic distances among A. hymenopus individuals (calculated from microsatellite data) and Bray-Curtis dissimilarities among their microbial communities (calculated from 16S amplicon data), while controlling for the effect of geographic distance. We used the 'mantel.partial' function in the vegan package, specifying the Pearson correlation statistic.

Population level data
We fitted population-level regressions (n = 6 sampling sites) to test the influence of allelic richness and Bd infection prevalence on alpha diversity. We did not fit both Bd prevalence and allelic richness in the same model to maximize residual degrees of freedom in the model. Significance of terms was assessed by likelihood ratio test.

Bd infection
The average Bd infection intensity across all individuals was 0.19 (±0.02) zoospores, with an average prevalence of 10.91 (±3.92)% across populations. There was relatively low infection prevalence and intensity at all sites, with highest infection rates at Tugela 2 and a complete absence of infection at Ribbon Falls (Table 2).
We observed significant differences in Shannon diversity index among sampling sites (Fig. 3). Tukelahed had the highest alpha diversity, whilst the values for the remaining five sampling locations were significantly lower (95% credible intervals of difference all >0; Fig. 3). There was no significant difference in the alpha diversity of the two Tugela sampling locations (95% credible interval of difference −0.92 to 0.05). NMDS ordination revealed distinct separation of bacterial community centroids corresponding to sampling site (Fig. 4). PERMANOVA of community distances supported this pattern, with significant separation among sampling sites (F 5,127 = 18.64, r 2 = 42.32%; p < 0.001). Shannon diversity of tadpole mouthparts was not affected by either Gosner stage of tadpoles (χ 2 = 2.61, p = 0.1), proportion of heterozygous loci (χ 2 = 2.19, p = 0.14), or their interaction (χ 2 = 0.19, p = 0.67).
Gosner stage significantly influenced within-population divergence in beta diversity (χ 2 = 5.01, p = 0.025; Fig. 5). There was no evidence of an effect of an interaction between Gosner stage and proportion of heterozygous loci on beta diversity, or proportion of heterozygous loci as a main effect (all p > 0.47). The best model explaining differences in beta diversity contained Gosner stage as the only fixed effect and had a marginal r 2 of 4%.
There was a positive correlation between host genetic distance calculated from microsatellites and microbial community distance when controlling for geographic distance (Partial Mantel Test r = 0.145, p < 0.001). There was no evidence that the Shannon diversity at a sampling site was affected by site-level allelic richness (F 1,4 = 0.95, p = 0.38) or Bd infection prevalence (F 1,4 = 1.17, p = 0.34).

Population genetics
Strong genetic structure was found in A. hymenopus, with F ST levels reaching 0.292 between Ribbon Falls and Nampolice (Table S2 and Supplementary Information). Ribbon Falls was highly differentiated from all sites, and in general sites from different rivers showed moderate levels of genetic differentiation, while sites within the same river showed low differentiation. STRUCTURE showed the optimum K to be two (Fig. 6), with Ribbon Falls forming a distinct population cluster and the remaining sites another. At K = 3, Vemvane and Nampolice also formed a separate cluster (Fig. 6). Hierarchical clustering analysis conducted on the population cluster of the K = 2 analysis (i.e., repeating the STRUCTURE run without Ribbon Falls) supported the presence of a Nampolice and Vemvane cluster. The PCoA and AMOVAs showed results consistent   with the STRUCTURE analysis (Fig. 7, Table S4), and together indicate that Ribbon Falls is a distinct population cluster, and Vemvane and Nampolice form a subpopulation cluster within the cluster of remaining sites.
Bottlenecks were detected at three sites: Heterozygosity excess in respect to allelic richness was found at Nampolice (p = 0.003), and deficiencies in M-ratios were found at Nampolice (p = 0.003), Vemvane (p < 0.001) and Ribbon Falls (p = 0.042).

Discussion
Population structure in A. hymenopus was strongly influenced by river network and geographic distance. However, at the landscape level, the microbiome of A. hymenopus is not structured according to host population genetic structure or river system. For example, Tukelahed did not have a similar microbial community structure to downstream sites Tugela and Tugela 2 (Fig. 4). Despite a lack of similarity between host population genetic structure and site-level microbiome structure, genetic distance among hosts is still significantly correlated with associated microbial community dissimilarity when controlling for geographic distance. This has not been previously shown in amphibians, but host genotype influences the gut microbiota of chickens [16], mice [63] and fish [7], and inbreeding and relatedness (although not population genetic structure) influence gut microbiomes of gopher tortoises (Gopherus polyphemus; [8]). Although not tested here, these patterns may be due to genetic variation in microsatellites reflecting variation in phenotypic traits that affect microbiome composition.  K h u b n a m N a m p o li c e V e m v a n e T u g e la 2 T u k e la h e d T u g e la B il a n ji l R ib b o n F a ll s K = 2 K = 3 Fig. 6 Bayesian population assignment analysis STRUCTURE plots showing membership coefficients (Q) for K = 2 and K = 3 genetic clusters in Amietia hymenopus. Vertical bars represent individuals, and colours represent a genetic cluster; the height of the bar represents the proportion of the genotype assigned to that cluster Microsatellites can be linked to genes coding for functional traits [64][65][66], including immune response genes in a number of species [67][68][69]. Variation in these genes could influence the composition of microbiomes by controlling the horizontal transmission of microbes from the environment. In A. hymenopus, the microbiome of the tadpole mouthparts is most likely transmitted from the environment as they are feeding or behaviourally selecting areas in the river. Behaviours can also influence symbiont transmission [70], and therefore there may be a link between genetic variation, behavioural variation, and microbe acquisition. That the influence of host genotype does not extend to landscape level microbiome structure may indicate that sitelevel environmental differences have more influence in structuring tadpole microbiomes, and host genotype simply exerts some selection pressure or behavioural influence on horizontal transmission. We detected a link between Gosner stage and mouthpart microbiome composition. Studies to date have predominantly used skin swabs to determine host microbiome of amphibians, as this is the region that becomes infected by Bd on postmetamorphic individuals. However, Bd primarily infects the mouthparts of larvae as keratin is localised to this area until the onset of metamorphosis when skin is remodelled to contain keratin. Although Bd infection intensities tend to increase as larvae develop due to increasing size [71], we show that alpha diversity (i.e., number of SVs) of the larval mouthpart microbiome does not increase as Gosner stage increases. This indicates that ontogenetic variation in microbiome is not simply an artefact of increasing mouthpart size, but rather it is the specific composition of the community that is changing as tadpoles develop. These microbiome changes may be related to phenotypic changes in mouthparts over development [24,72], which occur in other montane river species in the region during development as an adaptation to fast-flowing water [73]. Ontogenetic variation may also arise from changes in behaviour throughout larval development that means individuals are exposed to different environmental conditions and thus, microbial influences; however, this has not been studied in this species.
Despite these genotypic and ontogenetic influences on microbiome composition, populations shared a considerable core microbiome (Fig. 2), which has also been shown for other amphibian species [6,74], including among genetically isolated populations [75]. Our population genetics analyses detected low numbers of first-generation migrants between populations, with high-genetic differentiation given the relatively small geographic area studied. Together, these results indicate very low levels of dispersal among rivers, and in some instances, low dispersal among tributaries in the same river. This implies that core microbiome observed across the metapopulation of A. hymenopus has not arisen from contact between populations. A number of studies have now indicated that amphibians are able to select their microbiome based on functional traits [76,77]. This may be related to conserved host traits, such as anti-microbial peptide composition, or adaptive genetic loci. In addition, the river system and associated vegetation is remarkably homogenous on the uninhabited and unspoiled ecosystem of the Mont-aux-Sources plateau, with very low diversity of plants [78]. This may also in part explain the shared microbiome exhibited across populations. That said, there is likely to be fine-scale variation between environmental conditions in the rivers, which may be driving some variation in microbiome composition [6,79,80]. Microbial communities have been shown to vary between rivers even when colonisation substrate is controlled for [81]. Individuals from Nampolice and Tukelahed had the most distinct microbiome (Fig. 4), and the highest alpha diversity (Fig. 3), and this may be related to proximity from the source of the river. For example, Górniak et al [82] found bacterial diversity generally decreased as they sampled downstream from the source (glacial lake), and linked this to changes in environmental conditions and nutrient availability.
Historically there is wide annual and geographical variation in Bd infection prevalence in the Mont-aux-Sources region, ranging from~0 to 100% [23]. Although, we found uncharacteristically low infection loads in this study, our data correspond to historical patterns for the three sites that have undergone long-term monitoring. Ribbon Falls has historically low infection prevalence, followed by Bilanjil, and then Tugela, which shows the highest infection prevalence [23]. This may be driven by a range of environmental factors similar to those that drive microbiome diversity, including water flow, temperature, pH and depth [71]. Particularly low infection intensities such as those found here are, however, difficult to assess accurately using quantitative PCR. Given these low infection intensities, it is not possible to identify meaningful relationships between Bd infection and microbiome diversity, although this has been demonstrated in other studies ( [2,12,27,76,83,84]).
Ribbon Falls, Nampolice and Vemvane showed evidence of recent population bottlenecks, and these sites also have some of the lowest allelic richness. Considerable mortalities of A. hymenopus were seen at Ribbon Falls and Vemvane in September 2006 during routine monitoring (every 2-4 months) of Vemvane, Tugela, Bilanjil and Ribbon Falls ( [23]; C. Weldon, pers. obvs.). Diagnostic analyses identified Bd infections in 100% of the dead metamorphic A. hymenopus (n = 10) collected at each site. Further dead post-metamorphic individuals (n = 1-12) were encountered on most of the subsequent surveys, all infected with Bd [23]. Low genetic diversity in mortality-affected sites is unsurprising given the low level of migration from other sites. This presents a potential danger in bottlenecked locations struggling to recover evolutionary potential, impacting their ability to respond to stressors such as climate change, disease or local climatic events. Local chytridinduced mass mortality events may eliminate genetic diversity from the species as a whole due to the distinct population genetic structure and low migration. Interestingly, Bd infection prevalence at Ribbon Falls is generally low compared to other rivers in the region (this study and [23]), and Ribbon Falls has the lowest genetic diversity, which may be because of selection of Bd-resistant individuals. Addis et al. [85] showed montane boreal toads (Bufo boreas) with higher heterozygosity had higher Bd infection intensity, which they attribute to increased dispersal of heterozygous individuals leading to increased spread of Bd. Therefore, it is possible that the lack of immigrants to Ribbon Falls is also limiting the spread and maintenance of Bd at this site.

Conclusion
Here we have shown that the microbiome composition of amphibians is determined by genetic and ontogenetic variation, as well as geographic site, indicating within and among population predictors of microbiome diversity. Bd prevalence was low in the year we sampled, but we found lower genetic diversity and signals of genetic bottlenecks in sites where mass mortalities have been recorded in previous years, demonstrating the vulnerability to and persisting effects of disease in amphibian populations. Together these results demonstrate the key role of host genotype on microbial colonisation, which may be important in pathogen susceptibility in the amphibianchytridiomycosis system. Identifying individuals based on this characteristic may aid amphibian conservation efforts against this lethal pathogen.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.