Skip to main content

Assessing the causes and consequences of gut mycobiome variation in a wild population of the Seychelles warbler

Abstract

Background

Considerable research has focussed on the importance of bacterial communities within the vertebrate gut microbiome (GM). However, studies investigating the significance of other microbial kingdoms, such as fungi, are notably lacking, despite their potential to influence host processes. Here, we characterise the fungal GM of individuals living in a natural population of Seychelles warblers (Acrocephalus sechellensis). We evaluate the extent to which fungal GM structure is shaped by environment and host factors, including genome-wide heterozygosity and variation at key immune genes (major histocompatibility complex (MHC) and Toll-like receptor (TLR)). Importantly, we also explore the relationship between fungal GM differences and subsequent host survival. To our knowledge, this is the first time that the genetic drivers and fitness consequences of fungal GM variation have been characterised for a wild vertebrate population.

Results

Environmental factors, including season and territory quality, explain the largest proportion of variance in the fungal GM. In contrast, neither host age, sex, genome-wide heterozygosity, nor TLR3 genotype was associated with fungal GM differences in Seychelles warblers. However, the presence of four MHC-I alleles and one MHC-II allele was associated with changes in fungal GM alpha diversity. Changes in fungal richness ranged from between 1 and 10 sequencing variants lost or gained; in some cases, this accounted for 20% of the fungal variants carried by an individual. In addition to this, overall MHC-I allelic diversity was associated with small, but potentially important, changes in fungal GM composition. This is evidenced by the fact that fungal GM composition differed between individuals that survived or died within 7 months of being sampled.

Conclusions

Our results suggest that environmental factors play a primary role in shaping the fungal GM, but that components of the host immune system—specifically the MHC—may also contribute to the variation in fungal communities across individuals within wild populations. Furthermore, variation in the fungal GM can be associated with differential survival in the wild. Further work is needed to establish the causality of such relationships and, thus, the extent to which components of the GM may impact host evolution.

Video Abstract

Background

The vertebrate gut microbiome (GM) is a diverse microbial ecosystem of fundamental importance to the host. Species within the GM make critical contributions to many aspects of host biology, from nutrient acquisition to influencing host behaviour, immunity, and development [1,2,3]. GM structure can vary substantially amongst individuals living within the same natural population [4,5,6]. Given the importance of microbial species to host function, such variation may have significant consequences for traits associated with host fitness, including disease susceptibility [7, 8], survival [9, 10], and reproductive success [11, 12]. However, most GM research is biassed towards studying bacteria, the dominant taxonomic group (according to the proportion of sequencing reads) in the vertebrate gut [13]. In contrast, few studies have investigated the extent to which other microbial kingdoms in the GM, such as viruses, archaea, and microbial eukaryotes, vary across individuals living in wild vertebrate populations, and the potential drivers of that variation [14,15,16]. Furthermore, none of these studies has investigated the fitness consequences associated with variation in non-bacterial GM communities.

Fungi are an important but understudied component of the vertebrate GM, being present at lower abundance than bacterial taxa [13]. Despite this, evidence is accumulating that host-associated fungal communities, known as “the mycobiome”, can impact host health by contributing to important host processes. In the gut, fungal species can play an important role in the digestion of complex dietary components [17, 18], as well as the development and activity of the host immune system [19,20,21]. For example, laboratory experiments have demonstrated that increased fungal diversity in the GM can reduce the severity of diseases such as colitis in captive mice (Mus musculus) [22, 23], and certain fungal species can functionally replace intestinal bacteria in mitigating mucosal tissue injury [19].

Although numerous studies on captive animals have indicated the potential importance of the gut mycobiome to host health, very few studies have investigated gut mycobiome variation across individuals living in wild populations and the possible drivers of this variation [14, 24, 25]. Indeed, although a handful of studies have shown that environmental factors, such as season and habitat differences, influence fungal GM composition [14, 24, 25], there are no studies investigating the extent to which host genetic differences are associated with fungal GM variation across individuals living in wild populations. Moreover, we are not aware of any studies that have looked at the potential consequences of this mycobiome variation for host fitness components (e.g. differential survival) in the wild. This is despite growing evidence that captivity can significantly alter, and in many cases simplify, the GM [26, 27], including the mycobiome [21, 28]. Wild animals are exposed to highly complex and dynamic environments which are very difficult to replicate in captivity, but which shape the GM and interact with host fitness [29]. Furthermore, inbred captive hosts harbour reduced levels of genetic variation and, thus, are not fully representative of wild populations [29, 30]. As such, studying the extent to which host-associated mycobiome communities vary under natural conditions, the possible drivers of this variation, and the consequences for host fitness will be crucial if we are to fully understand the dynamics of the GM and the functional significance of the mycobiome to the host.

Variation in immune receptor genes, which recognise specific microbial antigens and activate an immune response, may play a key role in driving inter-individual differences in GM structure. A number of innate immune receptors respond to commensal and pathogenic fungi at mucosal surfaces, including C-type lectins, interleukins, and Toll-like receptors [19, 31]. Similarly, components of the adaptive immune system also interact with fungi and may play a key role in sculpting the microbiome. One such component is the major histocompatibility complex (MHC), a gene family present in virtually all vertebrate species [32,33,34]. MHC genes encode cell-surface glycoproteins which bind to protein antigens, including those derived from intestinal microbes. The resulting MHC-peptide complex is then presented to T-lymphocytes which, depending on the nature of the antigen, trigger an appropriate immune or tolerogenic response [35,36,37]. There are two major classes of MHC genes which process antigens via distinct pathways [33, 38]. MHC class I (MHC-I) genes are expressed on almost all somatic cells and predominantly process intracellular peptides (e.g. those derived from viruses). In contrast, MHC class II (MHC-II) genes are expressed on the surface of antigen-presenting cells and primarily bind to extracellular peptides, such as those derived from bacterial and fungal cells in the intestinal tract.

MHC genes exhibit extreme polymorphism, with many allelic variants existing in natural populations [32, 34, 39]. These polymorphisms can determine the range of antigen peptides (and therefore microbial species) that are recognised by the immune system [32]. As such, several studies have demonstrated that variation in MHC genotype may drive inter-individual differences in the bacterial GM [4, 40, 41]. For example, a previous study on the Seychelles warbler (Acrocephalus sechellensis) [4] demonstrated that specific MHC-I and MHC-II alleles are associated with shifts in the diversity and composition of the bacterial GM. In contrast, very little is known about the extent to which commensal fungal communities are shaped by variation at the MHC. Most research has focussed on specific fungal pathogen species, such as susceptibility to Batrachochytrium dendrobatidis infection in amphibians [42,43,44]. Comparatively few studies have looked at the influence of MHC variation on wider host-associated fungal communities [45, 46], and as far as we are aware, there are no studies investigating the relationship between MHC variation and fungi in the GM.

In this study, we investigate the environmental and genetic drivers of fungal GM variation in a natural population of the Seychelles warbler. We also investigate the relationship between fungal GM differences and host survival; to our knowledge, this is the first time this has been characterised in a natural vertebrate species. Importantly, the Seychelles warbler population on Cousin Island is closed, with virtually no inter-island dispersal, thus enabling accurate measures of survival to be achieved [47]. Although the population harbours reduced genetic variation as a result of a population bottleneck [48], variation still exists at the MHC-I and MHC-II loci [4, 49] and has been linked to survival differences [50]. In the Seychelles warbler, inter-individual variation in bacterial GM structure has been associated with age, sex, sampling season, genome-wide heterozygosity, and, importantly, MHC genotype [4]. Furthermore, differences in bacterial GM composition have been associated with survival [9]. Here, we expand our investigation of the GM to the mycobiome, by sequencing the fungal internal transcribed spacer subregion two (ITS2) ribosomal DNA.

We first test whether differences in fungal GM structure are associated with environmental variables, including season, territory quality, and time of day. We also test for an association with host variables, including sex, age, and genome-wide heterozygosity. Following this, we restrict our analysis to those individuals with complete immunogenetic data to test whether MHC diversity and/or specific MHC alleles are associated with variation in the diversity or composition of the fungal GM. In addition, we test for an association between mycobiome structure and TLR3 genotype, one of the few TLR loci for which functional variation has been maintained in the Seychelles warbler [51, 52]. The TLR3 genotype is correlated with host survival and reproductive success in this population [51]. Thus, although the TLR3 gene is thought to play a role in viral sensing [53], it may have an indirect impact on fungal communities, by influencing host susceptibility to viral infection or other aspects of host health. Finally, using the full dataset, we test for an association between variation in fungal GM structure and the probability of host survival to the next breeding season.

Methods

Study species and sample collection

The study was carried out on the population of Seychelles warblers inhabiting Cousin Island (29 ha; 04° 20′ S, 55° 40′ E), which has been monitored as part of a long-term research project since 1985 [51, 54, 55]. The population on Cousin consists of ca. 320 adult individuals which are distributed across ca. 115 year-round territories [55, 56]. Almost all individuals (> 96%) have been individually marked with a unique combination of a British Trust for Ornithology (BTO) metal ring and three plastic colour rings [57]. A population census is carried out twice a year in the minor (January to March) and major (June to September) breeding seasons [58]. As the annual resighting probability of adult individuals is high (98% ± 1% SE) [50] and inter-island dispersal is virtually absent [47], individuals not seen during a breeding season can be confidently assumed dead. Seychelles warblers have a median life expectancy at fledging of 5.5 years and a maximum recorded lifespan of 19 years [55, 59]. Since they lack natural predators and experience limited human disturbance, extrinsic mortality is lower than in many passerine species living in temperate regions [60]. Indeed, the average annual survival probability is exceptionally high in adults (0.84 ± 0.04 SE) and juveniles (0.61 ± 0.09) [60].

Individuals were caught in mist nests. Age was estimated based on a combination of lay, hatch, or fledge date or, if these dates were unknown, estimated from eye colour which changes from grey in fledglings to red-brown in adult individuals [54]. The mean age of individuals was 2.46 ± 2.88 SD (range 0.12–15.78 years). Blood samples were taken via brachial venipuncture and stored in absolute ethanol at 4 °C for later molecular analysis. As Seychelles warblers are insectivorous, an index of territory quality was calculated based on the number of insect prey available, the territory size, and foliage cover in an individual’s territory during that breeding season [54]. For territories with missing quality measures in a season, quality was calculated as the mean from the preceding and following breeding periods [as in 60].

Faecal samples were collected as described previously [see 4]. Briefly, captured birds were placed into a disposable, flat-bottomed paper bag containing a sterilised weigh boat protected by a metal grate; this set-up allows faecal matter to fall into the tray whilst minimising the possibility of contamination from the bird’s surface. Each bird was removed from the bag after defecation (or 30 min in the case of no defecation). The faeces were then collected into a sterile microcentrifuge tube containing 1 ml of absolute ethanol. Control samples were collected from empty collection bags and the fieldworker’s hands to capture possible sources of contamination (n = 7). All samples were stored at 4 °C for the remainder of the field season, before transferring to − 80 °C for long-term storage. The time kept at 4 °C was controlled for in all downstream analyses. A total of 277 faecal samples (one per individual) were taken across five sampling seasons (the minor breeding seasons of 2018 and 2019 and the major breeding seasons of 2017, 2018, and 2019) spanning 3 years. The survival of sampled individuals was assessed during the following breeding season (i.e. up to the minor breeding season of 2020).

Molecular methods

DNA was extracted from blood samples using the DNeasy Blood and Tissue kit (Qiagen, Crawley, UK) according to the manufacturer’s protocol. Molecular sexing was carried out using an established PCR-based method [59, 61]. Individuals were genotyped at 30 polymorphic microsatellite loci [57, 59], and standardised individual microsatellite heterozygosity (Hs) was calculated using genhet 3.1 [62] in R 4.1.1 [63]. Variation at a single non-synonymous SNP within the leucine-rich repeat domain of TLR3 exon 4, determined as part of a separate study, resulted in three genotypes (TLR3AA, TLR3AC, or TLR3CC) (see [51] for details). MHC-I exon 3 and MHC-II exon 2 were amplified and sequenced using Illumina MiSeq technology as part of a previous study (see [4] for methods). Complete MHC genotype data were available for 202 out of the 277 individuals with GM samples in the present study. Here, the term “allele” refers to different variants amplified for each class of MHC. Genotyping identified ten MHC-1 alleles that were present in > 5% but < 95% of individuals, and a further 10 alleles had a frequency of < 5%. Based on all alleles, individuals carried a mean of five MHC-I alleles (range 2–7 alleles). Three MHC-II alleles were present in > 5% but < 95% of individuals. A further two alleles were present in almost all individuals (> 97%), and nine alleles had a frequency of < 5%. Individuals carried a mean of 2.9 MHC-II alleles (range 1–5 alleles). All alleles were included in MHC diversity calculations but not in analyses of allele presence/absence, for which only alleles present in > 5% but < 95% of individuals were included.

Microbiome extraction and sequencing of the fungal community

Total genomic DNA was extracted from faecal samples and collection controls using the DNeasy PowerSoil Kit (Qiagen), following an optimised version of the manufacturer’s instructions (see [4] for detailed methods). As part of this protocol, samples were lysed using a combination of chemical lysis and bead beading. Chemical lysis was performed by means of the C1 lysis buffer (included in the kit). For bead beating, each sample was added to a PowerBead tube (provided) and secured to a vortex adapter; samples were then vortexed at maximum speed for 10 min. Modifications to the kit instructions included heating samples in the C1 buffer (65 °C for 10 min) prior to bead beating to improve lysis and eluting DNA in a final volume of 60 µl elution buffer. Samples were randomised across extractions. The bacterial component of the same extractions was sequenced as part of previous studies [4, 9]. Five randomly selected faecal samples were extracted twice so that the consistency of the extraction method for isolating fungal DNA could be assessed. One negative control was carried out per extraction kit using a blank sterile swab (n = 5 extraction blanks in total). Two types of positive control were also extracted to verify the completeness of the extraction protocol and reproducibility of sequencing methods. The first was extracted from a D6300 Microbial Community Standard (ZymoBIOMICS) which contains two fungal species (Saccharomyces cerevisiae and Cryptococcus neoformans) that are typically hard to lyse. The second was extracted from a synthetic fungal mock community containing 12 unique sequences [64]. The DNA concentration in each extraction was quantified using a Qubit dsDNA High Sensitivity Assay kit (Invitrogen).

Following DNA extraction, the ITS2 subregion of the fungal ITS rRNA was amplified and sequenced using a universal tail tag dual index barcoding approach [65]. For the first-round PCR, amplification was carried out in 96-well plates; 25-μl reactions were used for each sample. The primers gITS7_fw and ITS4_rev (Additional file 1: Table S1) were used for amplification [66, 67]. Primers were adapted to include recognition sequences at the 5′ end so that a secondary nested PCR could be carried out to incorporate Illumina adapter sequences and barcodes. The PCR mix contained 12.5 μl of 2 × Kapa HiFi Hotstart ReadyMix (Roche), 0.75 μl each of the forward and reverse primers (10 μM concentration), 8 μl of sterile water, and 3 μl of DNA template per reaction. Thermocycler conditions were as follows: 95 °C for 3 min; 20 cycles of 98 °C for 20 s, 55 °C for 15 s, and 72 °C for 60 s; and 72 °C for 10 min. One negative control (using 3 µl of sterile water instead of DNA template) was carried out per 96-well plate (n = 4 in total). First-round PCR products were then submitted to the Centre for Genomic Research (CGR), Liverpool, where they were cleaned using AMPure beads and entered into a second-round PCR (as described in [65]). Briefly, an Illumina index primer set (i7 and i5) comprising a unique 8 nucleotide index sequence and bases complementary to the sequence introduced in the first-round PCR was used to tag the library [65]. Amplification was carried out using the same reagents and conditions as in the first-round PCR, with 20 cycles. Amplicons were quantified via Qubit, and fragment distributions were assessed using an Agilent DNA high-sensitivity kit (2100 Bioanalyzer, Agilent). Amplicons were then pooled on an equimolar basis. The amplicon pool was then purified, and size was selected within the 300–700 bp range using PippinPrep. The fragment distribution was further checked via qPCR. Libraries underwent 2 × 300 bp paired-end sequencing on an Illumina MiSeq platform (Illumina, San Diego). Each library was run twice to increase sequencing depth; data from both runs were combined at the data pre-processing stage. A total of 300 samples were sequenced. This included 277 faecal samples and 5 extraction duplicates (thus, 282 faecal extractions in total), along with 7 collection controls, 2 positive controls (the D3600 standard and the synthetic mock community), 5 blank extraction controls, and 4 PCR negative controls.

Bioinformatic processing of fungal reads

The raw reads were trimmed at the CGR using Cutadapt 1.2.1 [68] to remove Illumina adapter sequences. Reads were further trimmed using Sickle 1.200 [69] with a minimum window quality score of 20. Data for the same sample (across the two MiSeq runs) were concatenated at this point. Upon receipt at UEA, sequences were imported into QIIME2 2019.10 [70] for further processing. The DADA2 plugin [71] was used to trim 6 base pairs from the 5′ end of reads to remove low-quality base calls. Reads were not truncated further as the ITS2 subregion can vary substantially in length across different fungal taxa [72]. Amplicon sequencing variants (ASVs) were inferred for each sample, followed by dereplication, paired-end joining, and the removal of chimeras and singletons (ASVs represented by a single read in the entire dataset). Two faecal samples which had very low read counts prior to DADA2 filtering contained zero reads after this step and so were removed. A total of 21,868,562 reads (mean per sample of 72,895.207 ± 88,523.704 SD) remained across 298 samples. Following this, ASVs were taxonomically classified by training a naïve-Bayes classifier on the developer version of the QIIME-compatible release of the UNITE database (version 8.3) [73]. A total of two ASVs were detected in the D6300 ZymoBIOMICs positive control; these corresponded to the two constituent fungal species Saccharomyces cerevisiae and Cryptococcus neoformans. Additionally, 12 ASVs with > 50 reads were detected in the synthetic mock community; these corresponded to the 12 synthetic sequences, suggesting that the extraction and sequencing methodologies had accurately captured this community. ASV and taxonomy tables were exported and further processed using phyloseq 1.36.0 [74] in R 4.1.1 [63].

Before conducting downstream analysis, sequences were filtered to remove ASVs that were unassigned at the phylum level. Additionally, potential contaminants were identified and filtered using the prevalence method in DECONTAM 1.12.0 [75]. First, laboratory contaminants were identified using PCR and extraction blanks as a reference; 7 ASVs were removed. Several of these were species of Candida which are known to contaminate reagents [76]. The second filtering step aimed to identify potential contaminants introduced at the sampling stage; collection controls taken from observer hands or from sample bags were used as a reference for this step. A total of 104 ASVs were filtered as possible contaminants—many of these were at a higher prevalence in low biomass samples. Following filtering, a total of 3441 fungal ASVs were detected across 279 faecal samples (a further low biomass faecal sample had zero reads following filtering). Sample completeness and rarefaction curves were generated using the R package iNEXT 2.0.20, with 50 bootstrap replicates per sample [77]. Sample completeness plateaued at approximately 5000 reads (Additional file 1: Fig. S1); all faecal samples with fewer than 5000 reads were therefore excluded from downstream analyses (14 samples). As the final filtering step, ASVs with fewer than 50 reads in total across all samples were also removed prior to downstream analysis, as these may represent possible sequencing errors. A total of 2540 ASVs were retained across 265 faecal sample extractions (mean ASVs per faecal sample = 51.69 ± 27.90 SD).

Statistical analyses

Alpha diversity analysis

Samples were rarefied to a depth of 5000 reads prior to calculating alpha diversity metrics to control for variation in library size across samples. A total of 2530 fungal ASVs remained across 265 samples following rarefaction. The observed ASV richness in each sample and the Shannon diversity index (which accounts for ASV abundance in a sample) were calculated using phyloseq 1.36.0 [74]. To assess whether the DNA extraction method was consistent, pairwise Euclidean distances were calculated based on the observed ASV richness in samples that had been extracted twice. The Kruskal–Wallis test was used to compare how pairwise distances varied according to the type of sample comparison, i.e. comparisons of duplicate extractions from the same sample, versus comparisons of extractions from different samples. Following this analysis, duplicates were filtered such that only the extraction with the highest read count was retained in downstream analyses, excepting one set of duplicates in which both were filtered to avoid pseudoreplication as this bird had been sampled again on a separate occasion (259 samples remained after filtering).

Following the calculation of alpha diversity metrics, correlations between fungal and bacterial alpha diversity metrics for the same samples were compared. We then analysed different subsets of our data using an information-theoretic (model averaging) approach [78, 79]. We first analysed the full dataset (n = 259 samples/individuals). Four “floater” individuals were removed prior to analysis as these had no assigned territory (n = 255 samples/individuals remaining). To establish whether fungal alpha diversity varied according to host or environmental factors (not including immunogenetic variables), a global linear model with a Gaussian distribution was constructed using stats 4.1.1 [63]. Shannon diversity was used as the response variable. Host sex, age, genome-wide heterozygosity (Hs), breeding season (major or minor), territory quality, and the amount of time the sample was stored at 4 °C prior to freezing were included as independent variables. Squared terms were also included for age and Hs but were subsequently removed due to a lack of significance (P > 0.05) and to enable the interpretation of first-order effects. Continuous variables were centred and scaled using arm 1.12.2 [80]. Collinearity between independent variables was checked using the vif function in car 3.0.12 [81]; vifs were < 2 for all variables in the model. Further model diagnostics were carried out using DHARMa 0.4.5 [82]. The dredge function in MuMIn 1.44.3 [83] was used to build all plausible models and rank them by AICc scores. All models with AICc scores within 7 of the top model were included in the final model set and used to obtain the conditional averaged estimates [78]. To assess whether results were consistent across different alpha diversity metrics, a second model was run with observed ASV richness as the response variable; a generalised linear model (GLM) was constructed with a negative binomial distribution using glm.nb in MASS 7.3.54 [84] to account for ASV counts being over-dispersed. The same independent variables were included in this model.

For the second set of analyses, the dataset was restricted to only include individuals with complete MHC and TLR3 genotype data (n = 189 samples/individuals). Two model types were generated for each alpha diversity metric; this was to establish whether fungal alpha diversity varied according to (1) MHC diversity and (2) the presence/absence of specific MHC alleles. The first model contained MHC-1 and MHC-II diversities as independent variables, as well as their squared terms, since optimal diversity could be more beneficial [85]. The squared terms were subsequently removed due to a lack of significance (P > 0.05) and to enable the interpretation of first-order effects. The second model included the presence/absence of each MHC allele. For MHC-I, these were Ase-ua1, 3–8, and 11. Ase-ua10, and Ase-ua1 were highly correlated (98% co-occurrence), so only Ase-ua1 was included in the analysis. The MHC-II alleles were Ase-dab3, 4, and 5. Both model types also included heterozygosity (Hs) and TLR3 genotype (TLR3AA, TLR3AC, or TLR3CC). Host sex, age, breeding season (major or minor), territory quality, time of day (AM or PM), and time stored at 4 °C were also controlled for in each model. Since squared terms for age and Hs were not significant predictors in the larger model set, these terms were not included to aid interpretation of the first-order effects. In both models, vifs were < 5 for all variables. Model averaging was carried out as above.

Finally, to assess whether fungal alpha diversity was associated with differences in host survival, a generalised linear model with a binomial error structure and logit link function was constructed using stats 4.1.1 [63]. Survival to the next breeding season (1, yes; 0, no) was included as the dependent variable. Fungal alpha diversity, host sex, age, and territory quality were included as independent variables. Sample year (2017, 2018, or 2019) was included to control for the differences in survival probabilities between years [60]. Squared terms for age and alpha diversity were tested but were subsequently removed due to a lack of significance. This analysis was conducted using the full dataset (n = 259 samples; 37 individuals died, 218 survived to the following season). The significance of terms was assessed as above. Although the exact date of death is not known, for 25 of the individuals that died by the next breeding season, the point of GM sampling was the last time they were observed in the population. The remaining 12 individuals that died were observed (but not sampled) again in the same breeding season as GM sampling took place; in some cases, this was up to 7 weeks after their last GM sample was taken. There was also a median period of 4 months (range 3–7) between the point when GM samples were taken and when the population was next censused to assess survival. Thus, it is possible that some of the individuals were sampled several months before their point of death.

Beta diversity

The unrarefied reads (filtered to remove samples with < 5000 reads) were used in beta diversity analyses. ASVs were further filtered from the dataset if they were present in less than 1% of samples (approximately 3 samples) to remove exceptionally rare taxa that may disproportionately influence beta diversity metrics. There were 794 ASVs remaining across 259 samples after filtering. ASV abundances were transformed using the centred log ratio (CLR) transform function in microbiome 1.14.0 [86]. This produces values that are scale invariant (i.e. not influenced by differences in library size across samples) and controls for the compositional nature of microbiome datasets [87]. To assess whether the DNA extraction method produced consistent results for beta diversity measurements, pairwise Euclidean distances were calculated, using the CLR-transformed ASV abundances, between duplicate extractions of the same sample or between extractions of different samples. Distances were compared using the Kruskal–Wallis test. Extraction duplicates were then removed, as described for alpha diversity analysis. The four “floater” individuals were also removed, as described above.

To test whether fungal GM composition varied according to host and environmental variables, a permutational analysis of variance (PERMANOVA) was performed on a Euclidean distance matrix of CLR-transformed ASV abundances, using the adonis2 function in vegan 2.5.7 [88, 89], with 9999 permutations. We first analysed the full dataset (n = 255 individuals) to maximise power. The host’s age, sex, heterozygosity (Hs), breeding season (major or minor), territory quality, and time of sampling (AM or PM) were included in the analysis. Survival to the next breeding season (yes or no) was included as a factor in the model. The amount of time each sample was stored at 4 °C in the field was also controlled for in the analysis. The function betadisper was used to check for the homogeneity of group dispersion values [88, 89]. Differences in fungal GM composition between the groups were visualised using a principal components analysis (PCA). We further tested for associations between variables and shifts across specific principal component axes using the envfit function in vegan 2.5.7; this uses linear model permutations to map variables onto an ordination.

In the second analysis, the dataset was restricted to include individuals with complete MHC and TLR3 genotype data, as described for alpha diversity analysis (n = 189 individuals). Two different models were run to test whether fungal GM composition was associated with immunogenetic variation. Firstly, a PERMANOVA was performed as described above for the full dataset, but with MHC-I and MHC-II diversities as additional predictor variables. The second PERMANOVA model instead included the presence/absence of each MHC allele as predictor variables, as described for alpha diversity analysis. TLR3 genotype (TLR3AA, TLR3AC, or TLR3CC) was also included in both models. All other variables included in the full analysis were controlled for. Neither MHC diversity, nor the presence/absence of specific alleles correlated with survival to the next breeding season (P > 0.1 in GLMs), and so these terms could be included in the same PERMANOVA analysis.

To establish whether specific ASVs were differentially abundant across particular groups of individuals, an analysis of compositions of microbiomes with bias correction (ANCOM-BC) was carried out, using ANCOMBC 1.1.5 [90]. As part of ANCOM-BC, the Benjamini and Hochberg method was used to correct P values for multiple testing [91]. A cutoff of Padj < 0.05 was used to assess significance.

Results

Characterising the fungal GM of the Seychelles warbler

Following quality filtering and the removal of samples with < 5000 reads, the number of high-quality fungal reads ranged from 5013 to 1,204,831 across the remaining 265 faecal samples. A total of 2540 fungal ASVs were detected across all samples, with a mean of 51.69 ± 27.90 (SD) ASVs inferred per sample. Comparisons of samples that had been extracted twice showed that the extraction method was consistent; pairwise distances were significantly lower between repeat extractions of the same sample, than between extractions from different samples (P < 0.05 in Kruskal–Wallis tests, Additional file 1: Fig. S2). This was the case for both observed ASV richness and fungal GM composition (Additional file 1: Fig. S2).

Following the removal of duplicate extractions and rarefaction to 5000 reads, 2530 fungal ASVs remained across 259 samples (mean ASV richness per sample = 49.49 ± 26.02 SD). The mean fungal Shannon diversity was 2.29 ± 0.66 across samples. On average, a greater number of bacterial ASVs were detected per sample than fungal ASVs (mean of 278.53 ± 167.77 SD per sample). The mean bacterial Shannon diversity was 3.86 ± 1.15. There was a significant, positive correlation between the fungal and bacterial ASV richness identified in each sample (r = 0.50, P < 0.001; Additional file 1: Fig S3). A significant, positive correlation was also identified between fungal and bacterial Shannon diversity measures (r = 0.18, P = 0.003; Additional file1: Fig S3).

A total of seven fungal phyla were detected across samples (Fig. 1a). The phylum Ascomycota was the most dominant and was present in 100% of samples with a mean relative abundance of 89.29% ± 11.69% (SD). This was followed by Basidiomycota which was present in 97% of samples with a mean abundance of 10.01% ± 11.07% (SD) across samples. All other phyla were present in fewer than 10% of samples. At lower taxonomic levels, the 2530 detected ASVs grouped into a total of 39 fungal classes, 249 families, and 475 genera.

Fig. 1
figure 1

The relative abundance (%) of fungal A phyla and B families in Seychelles warbler faecal samples (N = 259 individuals). Each vertical bar represents a separate faecal sample. Samples are ordered according to the abundance of (a) Ascomycota and (b) Cladosporiaceae. Families with a median relative abundance of < 1% across samples are collapsed into the category “other”

Core fungal taxa were identified as those that were present in at least 50% of samples at a minimum relative abundance of 0.1%. A total of 13 core families were identified, four of which had a prevalence of over 80% (Additional file 1: Table S2). These were Cladosporiaceae (mean relative abundance across all samples 34.11% ± 23.29% SD), Dothideales incertae sedis (10.79% ± 11.96%), Aspergillaceae (7.43% ± 12.72%), and Bionectriaceae (3.48% ± 8.56%) (Additional file 1: Table S2). A total of 13 core genera were also identified, three of which were present in more than 80% of samples: Cladosporium (34.03% ± 23.27%), Hortaea (10.79% ± 11.95%), and Penicillium (3.79% ± 8.78%) (Additional file 1: Table S2). Further to this, a total of 7 core ASVs were identified and classified as the following at species level: Cladosporium dominicanum (28.57% ± 20.77%), Hortaea werneckii (28.57% ± 20.77%), Cladosporium coloradense (3.06% ± 5.9%), Penicillium citrinum (2.13% ± 6.31%), an unidentified species in the family Didymellaceae (2.19% ± 5.38%), and two ASVs identified as Sympodiomycopsis kandeliae (1.96% ± 3.03% and 1.23% ± 2.55%) (Additional file 1: Table S2). Despite the presence of shared taxa, the composition and diversity of the fungal GM varied substantially across individuals in the population (Fig. 1).

The impact of host and environmental variables on fungal GM alpha diversity

We first analysed the full dataset (n = 255 individuals) to establish whether fungal GM alpha diversity varied according to host and/or environmental factors, before testing for the effects of immunogenetic variables using a reduced dataset. Fungal alpha diversity was significantly higher in the minor versus the major breeding season (Table 1, Fig. 2). However, neither territory quality nor time of sampling was significantly associated with fungal alpha diversity (Table 1). Host age, sex, and genome-wide heterozygosity were also not associated with fungal alpha diversity (Table 1). The length of time that samples were stored at 4 °C had no effect on fungal Shannon diversity (Table 1), but there was a weak, negative association between storage time at 4 °C and the fungal ASV richness detected in samples (Table 1). As such, storage time at 4 °C was controlled for in all downstream analyses.

Table 1 The impact of host and environmental factors on fungal alpha diversity in the gut microbiome of the Seychelles warbler (n = 255 individuals). (A) Shannon diversity and (B) observed ASV richness were used as the response variable in two separate models. Estimates and standard errors (SE) are based on linear conditional model-averaged estimates. Significant (P < 0.05) predictors are shown in bold. The reference categories for categorical variables were as follows: female (sex), major (season), and AM (time of day)
Fig. 2
figure 2

The impact of season on fungal A Shannon diversity and B ASV richness in the gut microbiome of the Seychelles warbler (N = 255 individuals). Major season = 182 samples, minor season = 73 samples. Black points represent the raw data. Boxes span the interquartile (25–75%) range. The median is marked by a horizontal line, and the mean is marked by a red diamond. Whiskers extend to 1.5 times the interquartile range. **P < 0.001, *P < 0.01 from a linear model

The effect of immunogenetic variation on fungal GM alpha diversity

There was no relationship between individual MHC diversity and fungal GM alpha diversity (when measured as either Shannon diversity or observed ASV richness; Additional file 1: Table S3). TLR3 genotype also had no significant effect on GM alpha diversity (Additional file 1: Table S3). In contrast, the presence of three MHC-1 alleles and one MHC-II allele was significantly associated with changes in ASV richness (Fig. 3). Individuals that carried the MHC-I allele Ase-ua4 or Ase-ua7 had significantly lower ASV richness compared to individuals without these alleles (Fig. 3, Additional file 1: Table S4 and Fig. S4a,b). In contrast, individuals that carried the MHC-I allele Ase-ua8 had greater fungal ASV richness compared to individuals without this allele (Fig. 3, Additional file 1: Table S4 and Fig. S4c). This was also the case for individuals possessing the MHC-II allele Ase-dab3 (Fig. 3, Additional file 1: Table S4 and Fig. S4d). The mean number of ASVs (± SD) detected in samples when Ase-ua4 was absent was 50.34 ± 24.73, compared to 48.91 ± 24.52 when the allele was present (Additional file 1: Fig. S4). In contrast, changes were considerably greater for the allele Ase-ua7; a mean of 57.96 ± 26.81 ASVs was detected when the allele was absent compared to 47.61 ± 24.51 when present (c.a. a 20% loss in ASV richness). The same was also true for the alleles Ase-ua8 and Ase-dab3 (Ase-ua8: 48.59 ± 24.81 absent, 55.97 ± 27.15 present; Ase-dab3 47.22 ± 21.35 absent, 57.80 ± 33.06 present). However, none of these alleles was significantly associated with shifts in Shannon diversity (Fig. 3, Additional file 1: Table S4). Instead, there was a marginally significant, negative relationship between the presence of Ase-ua11 and Shannon diversity (Fig. 3, Additional file 1: Table S4).

Fig. 3
figure 3

The effect of MHC alleles and TLR3 genotype on fungal alpha diversity in the gut microbiome of the Seychelles warbler, controlling for other host and environmental variables. Diversity was measured as A observed ASV richness and B Shannon diversity. N = 189 individuals. Estimates and standard errors are based on linear conditional model-averaged estimates. An estimate of greater or less than zero indicates increased or reduced alpha diversity, respectively. Significant terms (P < 0.05) are highlighted in yellow. Reference categories for categorical variables were as follows: absent (for all MHC alleles), TLR3.AA (TLR3 genotype), female (sex), major (season), and AM (time of day)

The impact of host and environmental variables on fungal GM composition

A PERMANOVA analysis using the full dataset revealed that fungal GM composition differed significantly between the major and minor breeding seasons (Table 2, Fig. 4); fungal GM composition was also more variable in the minor breeding season (betadisper test: F1,253 = 19.654, P < 0.001, Additional file 1: Fig. S5). Shifts in GM composition were additionally associated with differences in territory quality (Table 2, Fig. 4) and, weakly, with the time of day that an individual was sampled (Table 2, Fig. 4). GM composition also varied significantly according to the amount of time samples were stored at 4 °C (Table 2) indicating the importance of controlling for sample storage methods. Differences in host age, sex, and genome-wide heterozygosity were not associated with shifts in GM composition (Table 2).

Table 2 PERMANOVA analysis of fungal gut microbiome communities in the Seychelles warbler. Euclidean distances were calculated based on centred log ratio (CLR)-transformed amplicon sequencing variant (ASV) abundances. Significant predictors (P < 0.05) are shown in bold. The analysis included 255 individuals
Fig. 4
figure 4

Shifts in fungal gut microbiome composition according to A time of day, B sampling season, and C territory quality in the Seychelles warbler. PCA ordination was carried out using Euclidean distances calculated using centred log ratio (CLR)-transformed amplicon sequencing variant (ASV) abundances. Each point represents a unique gut microbiome sample (N = 255). Large circles represent the group centroids. For clarity, samples were grouped into discrete categories for territory quality (high or low depending on whether it was above or below the median territory quality of 12,678, respectively). Principal components 1, 2, 3, and 4 explained 4.31%, 3.54%, 2.90%, and 2.48% of the variation in gut microbiome structure, respectively

Points clustered on a PCA ordination according to significant environmental variables (Fig. 4); clustering along the first two principal component axes was significantly associated with the time of day (envfit model: R2 = 0.018, P = 0.010, Fig. 4a) and sampling season (R2 = 0.064, P < 0.001), although points clustered more strongly according to season along axes three and four (R2 = 0.106, P < 0.001, Fig. 4b). Placement of points along the third and fourth axes was also influenced by differences in territory quality (R2 = 0.048, P = 0.003, Fig. 4c).

Differential abundance tests revealed that 31 ASVs were significantly differentially abundant between the major and minor breeding seasons (P < 0.05 in an ANCOM-BC test; Additional File 1: Fig. S6 and Additional file 2). Of these, seven were more abundant in the major season. These were in the fungal orders Capnodiales (two ASVs in the genus Cladosporium and one ASV in the genus Recurvomyces) and Eurotiales (three ASVs in the genus Penicillium, one in the genus Aspergillus). The remaining 24 ASVs, belonging to ten identified fungal orders, were significantly more abundant in the minor breeding season (Additional File 1: Fig. S6 and Additional file 2). The fungal orders Pleosporales (10 ASVs) and the Hypocreales (three ASVs) contributed over 50% of these ASVs.

There were seven ASVs that were differentially abundant between high and low territory qualities (see Additional file 2 for a full taxonomic breakdown). Quality scores above or below the median territory quality (12,678) were classified as high or low, respectively, in this analysis. Six of the ASVs were more abundant in high territory qualities (Additional file 2); five of these were in the fungal phylum Basidiomycota and belonged to two orders, namely, Tremallales (one ASV in the genus Dioszegia and one in an unidentified genus) and Microstromatales (two ASVs in the genus Sympodiomycopsis and one in the genus Jaminaea, respectively). The other ASV was in the phylum Ascomycota (an unidentified member of the family Didymellaceae). One ASV in the phylum Ascomycota was more abundant in low-quality territories (Additional file 2); this was in the order Capnodiales (unidentified below order level). Differential abundance testing did not identify any ASVs that were significantly differentially abundant between samples taken in the morning or afternoon.

The impact of immunogenetic variation on fungal GM composition

A PERMANOVA analysis using the immunogenetic dataset demonstrated that there was a marginally significant association between MHC-I diversity and fungal GM composition (Table 3), but MHC diversity was not associated with differences in GM variability (betadisper test: F1,187 = 1.018, P = 0.405). MHC-I diversity was associated with the placement of points along the first two axes of a PCA ordination (envfit model: R2 = 0.044, P = 0.015; Fig. 5). Testing identified two ASVs that were differentially abundant across individuals with low (2–4 alleles) or high (5–7 alleles) MHC-I diversity (P < 0.05 in an ANCOM-BC test; Additional file 2). One ASV was more abundant in individuals that had low diversity; this was in the order Botryosphaeriales (genus Lasiodiplodia). The second ASV was more abundant in individuals that had high MHC-I diversity; this was in the order Saccharomycetales (genus Sympodiomyces). MHC-II diversity and TLR3 genotype had no effect on GM composition (Table 3). A separate PERMANOVA analysis showed that the presence/absence of specific MHC-I or MHC-II alleles was also not associated with variation in fungal GM composition (Additional file 1: Table S5).

Table 3 The influence of MHC diversity on fungal gut microbiome beta diversity in the Seychelles warbler. Based on a PERMANOVA analysis of Euclidean distances, calculated using centred log ratio (CLR)-transformed amplicon sequencing variant (ASV) abundances. Significant predictors (P < 0.05) are shown in bold. The analysis included 189 individuals
Fig. 5
figure 5

Shifts in fungal gut microbiome (GM) beta diversity according to MHC-I diversity in the Seychelles warbler. Principal components analysis (PCA) was carried out using Euclidean distances based on CLR-transformed abundances of amplicon sequencing variants (ASVs). Principal components 1 and 2 explained 4.32% and 3.66% of the variation in GM community structure, respectively. N = 189 individuals were included in the analysis. For clarity, individuals were grouped into discrete categories, according to whether they had 2–4 (< 5, green points) or 5–7 (5 + , yellow points) alleles. Large circles represent group centroids

Fungal GM structure and host survival

There was no significant association between fungal GM alpha diversity and the probability that an individual survived to the next breeding season (Additional file 1: Table S6) for either the Shannon diversity or observed ASV richness metric. However, survival to the next breeding season was significantly associated with shifts in fungal GM composition in a PERMANOVA analysis (Table 2, Fig. 6). Points predominantly separated along the PC3 and PC4 axis of a PCA ordination according to survival (envfit model: R2 = 0.025, P = 0.002, Fig. 6). A betadisper test indicated that there was no difference in fungal GM variability between individuals that died and those that survived by the following breeding season (F1,253 = 0.013, P = 0.910).

Fig. 6
figure 6

The association between fungal gut microbiome (GM) structure and differential survival in the Seychelles warbler. Principal components analysis (PCA) was carried out using Euclidean distances based on CLR-transformed abundances of amplicon sequencing variants (ASVs). Points are coloured according to whether individuals survived (blue) or died (red) by the next breeding season. Principal components 3 and 4 explained 2.90% and 2.48% of the variation in GM community structure, respectively. N = 255 individuals were included in the analysis (218 individuals survived, 37 individuals died). Large circles represent group centroids

Differential abundance testing indicated that six fungal ASVs were more abundant in individuals that survived to the following breeding season, compared to those that died (P < 0.05 in an ANCOM-BC test, Additional file 2). Three of these belonged to the fungal class Sodariomycetes; these were members of the genera Gliomastix, Chordomyces, and Acrostalagmus. A further two ASVs belonged to the order Microstromatales (both in the genus Sympodiomycopsis). The remaining ASV belonged to the order Venturiales (genus Ochroconis). No ASVs were identified that were significantly more abundant in individuals that died.

Discussion

In this study, we characterised variation in the fungal GM in a natural population of the Seychelles warbler. We then investigated the extent to which this variation was associated with environmental and host factors and, importantly, subsequent host survival. Variation in fungal GM composition was associated with breeding season, territory quality, and the time of day at which an individual was sampled. In contrast, there were few host variables associated with fungal GM variation. Neither host age, sex, genome-wide heterozygosity, or TLR3 genotype were associated with differences in fungal GM structure. However, the presence of three MHC-I alleles and one MHC-II allele was associated with changes in fungal ASV richness, whilst the presence of a different MHC-I allele was associated with a reduction in Shannon diversity. Despite these associations, differences in diversity based on allele presence were reasonably small and did not scale up to dissimilarities in GM composition. However, overall MHC-I allelic diversity was (to a limited extent) associated with fungal GM composition. Finally, variation in fungal GM composition was associated with differential survival, and the abundances of specific fungal ASVs differed between individuals that survived versus those that died by the following breeding season.

The fungal GM of the Seychelles warbler

The composition of the Seychelles warbler fungal GM was similar to that of other vertebrate taxa, including passerine species, in that it was dominated by the phyla Ascomycota and Basidiomycota [92]. Despite substantial variation in fungal GM structure across individuals in the Cousin Island population, we identified several core taxa that were shared across the majority (> 50%) of individuals. Several of these taxa are common members of the vertebrate GM. For example, members of the family Aspergillaceae (including Aspergillus and Penicillium species) are often identified in GM samples taken from human and non-human primates [14, 93,94,95]. The genus Cladosporium is one of the most dominant filamentous fungal taxa in the human GM [93, 96, 97] and is abundant in faecal samples taken from other vertebrate species [94, 98]. This genus is also more abundant in the GM of insectivorous, versus phytophagous, bat species suggesting that there may be some association between its abundance and having an insect-based diet [98].

It is likely that many of the core taxa identified in the Seychelles warbler GM are environmentally acquired. For example, Cladosporium species are widely distributed in the environment and have been isolated from a diverse range of habitats including soils, plant material and insects [99]. Similarly, species in the genus Sympodiomycopsis, have previously been isolated from insects [100] and mangrove plants [101]. Thus, Seychelles warblers may acquire such species by ingesting their insect prey. Since birds possess relatively short intestinal tracts as an adaptation to improve flight efficiency, they may be particularly prone to acquiring microbial strains from their environment. Thus, these species may only transiently pass through the gut ecosystem before being excreted in faecal matter [102,103,104].

Differentiating between resident gut commensals and transient colonisers is challenging given our limited understanding of host-microbe relationships within the intestinal tract and the fact that many fungal strains remain poorly characterised [13, 105, 106]. Furthermore, the extent to which transient strains can elicit changes in the resident gut microbiota and perform functions relevant to the host remains unclear [105, 107, 108]. In future studies, it will be important to establish the degree to which different fungal species consistently colonise their host and the functionality of these strains within the gut ecosystem [105, 106]. Although such studies may require experimental manipulation, longitudinal sampling of wild individuals could help to determine which microbial taxa are stably acquired, and metagenomic sequencing would provide further information on their functionality.

The selection of primers for fungal community analysis remains a subject of much debate since primer choice can significantly influence community composition [72, 109, 110]. The specific ITS2 primers used in this study were chosen as they provide superior coverage compared to other ITS2 primers (by amplifying > 90% of fungal groups) and those targeting the ITS1 region which, in comparative studies, have been shown to introduce greater taxonomic bias and underestimate levels of fungal diversity [72, 111]. The ITS2 region can also provide greater taxonomic resolution than rRNA markers (such as the 18S subunit) [72, 110]; this is likely to be important when assessing inter-individual variation at a fine scale. However, we acknowledge that 18S amplicons can be useful for calculating measures of phylogenetic diversity given the possibility of sequence alignment across a broad range of fungal taxa [72]. The primers used in this study were shown to accurately capture two mock community controls of varying complexity. However, biases could still exist, and combining different primer sets may allow fungal communities to be captured with further accuracy in the future.

Environmental drivers of fungal GM variation

Season explained the largest proportion of variance in fungal community composition across sampled individuals. Fungal GM alpha diversity and variation in GM composition were also greater in the minor, compared to the major, breeding season. Climatic variables can vary considerably between seasons on Cousin Island, for example, rainfall tends to be greater in the minor season [58]. Such variation likely impacts microbial communities in the external environment [112] meaning that hosts are exposed to different fungal species either directly or indirectly via their insect prey. Differing weather conditions could also have an impact on island-wide insect availability which, by altering the type or quantity of diet available, could influence host condition and/or conditions within the intestinal tract and, in turn, the fungal GM [113]. Indeed, previous research has shown that island-wide insect availability tends to be significantly greater during the major season on Cousin Island [114]. A similar effect of season on fungal GM communities has been observed in wild Tibetan Macaques (Macaca thibetana), in which alpha diversity and certain fungal taxa increased in winter months; these differences were thought to be driven by changing food availability and dietary components [14].

The reproductive behaviour of the warblers also changes notably across the year, with the majority of birds breeding in the major season [58]. As such, changes in fungal GM structure between seasons could also be driven, in part, by physiological or behavioural changes related to reproduction. Relationships between breeding activity (e.g. hormonal changes, nest building, and copulation) and the microbiota have been identified when studying the bacterial GM in other bird species [11, 115, 116]. Several of the fungal taxa that increased in abundance in the GM of Seychelles warblers sampled in the major season are known to be plant-associated and thus could have been acquired from nesting material. For example, Cladosporium species have been isolated from mangrove plants in tropical regions [117], and species in the genus Penicillium are known to live on, and decompose, plant material [118, 119].

Fungal GM composition also differed according to the quality of the host’s territory. As with seasonal differences, this relationship could arise due to environmental filtering, whereby the environmental conditions within a particular territory determine the range of fungal species that can proliferate and colonise the host [120, 121]. Thus, birds inhabiting territories that are exposed to differing levels of foliage cover, humidity and/or salt spray, may come into contact with different fungal taxa. However, variation in territory quality has also been associated with differences in host condition in the Seychelles warbler [122, 123]. Experimental manipulations and studies of wild animals have shown that prolonged stress, differing food availability, and host infection status can all have an impact on the GM, including the composition of fungal communities [124,125,126]. Thus, physiological differences associated with varying territory qualities may have consequences for fungi in the intestinal tract.

Host genetics and fungal GM variation

Compared to environmental variables, host genetic factors appeared to have a weaker influence on fungal GM structure. Genome-wide heterozygosity was not associated with shifts in GM composition across individuals. This is surprising given that heterozygosity is a measure of inbreeding. Inbreeding occurs frequently in the Seychelles warbler, and previous studies have identified an effect of inbreeding depression on fitness in this population [127,128,129]. Thus, it might be expected that, via an effect on host physiology and condition, heterozygosity would also be associated with GM differences. However, some of the costs of inbreeding have been found to be environmentally dependent in this species, whereby they are only observed under very poor conditions [127,128,129]. Thus, further years of data collection, which incorporate temporal variation in environmental conditions, may be needed to identify the effects of inbreeding on the GM.

It is perhaps not surprising that TLR3 variation had no influence on the fungal GM, since it encodes receptor molecules that bind to viral dsRNA [53]. It is possible that polymorphisms at other TLR loci may have a greater influence on fungal GM variation in wild populations. For example, molecules encoded by the TLR4 gene have been shown to specifically detect fungal peptides [19, 130]. However, the majority of TLR loci (including TLR4) lack variation in the Seychelles warbler [52].

The presence of two MHC alleles was associated with a significant reduction in fungal ASV richness (Ase-ua7 and Ase-ua4), whilst one MHC-I allele (Ase-ua8) and one MHC-II allele (Ase-dab3) were associated with an increase in ASV richness. However, these effects were limited; the mean number of ASVs that were lost or gained ranged from between two and ten ASVs, though in the latter case, this accounted for c.a. 20% of fungal ASVs carried by an individual. There was also no association between the presence of these alleles and Shannon diversity, suggesting that only a small number of rare ASVs may have been gained or selectively eliminated, rather than these alleles having a wider influence on taxon abundances and community evenness. Only one different allele (Ase-ua11) had a marginally significant influence on Shannon diversity. However, the presence/absence of these individual alleles did not have a significant influence on overall fungal GM composition. Such small effects are not unprecedented; a study on the bacterial GM of wild mouse lemurs (Microcebus griseorufus) found that the majority of MHC-I and MHC-II supertypes influenced the abundance of a very small number of ASVs (between one and three taxa), rather than having an overarching effect on the entire GM community [41]. This also concurs with an experimental study on laboratory mice (Mus musculus) which found that MHC genotype had a very limited effect on the GM; this study concluded that the immune activity of the MHC may be reserved for a small number of specific colonisers such as pathogens, or intrusive GM commensals [131]. More generally, small effect sizes may be expected when investigating relationships between MHC variation and immunocompetence, given that the MHC plays only a small part in a complex network of immune pathways [132].

The diversity of MHC-I alleles (but not MHC-II alleles) was associated with shifts in GM community composition. This is despite that fact that none of the MHC-I alleles was individually associated with significant compositional shifts. However, the association between MHC-I diversity and GM composition was relatively weak, explaining 0.6% of the variance in GM community structure. Furthermore, only two ASVs were identified as being differentially abundant between individuals with high and low MHC-I diversity; again, this is consistent with the idea that the MHC may act to target specific pathogens or tolerate specific commensal strains, rather than having pervasive effects on the whole GM [41, 131]. High MHC diversity was associated with a reduction in the abundance of an ASV in the fungal family Botryosphaeriaceae. Members of the Botryosphaeriaceae have predominantly been found in association with plants, often as pathogens [133], but certain genera can cause opportunistic infections in animal species [134]. Thus, it is possible that this species was more likely to be detected and selectively eliminated from the GM of individuals carrying a greater repertoire of MHC-I motifs. The second differentially abundant ASV—a member of the genus Sympodiomyces (order Saccharomycetales)—was more abundant in the GM of individuals with high MHC-I diversity. Although little is known about the Sympodiomyces, members of the Saccharomycetales are common in adult human faeces [135] and may play a beneficial role in maintaining host metabolic health [136, 137]. However, further research is needed to characterise the function of these ASVs in the avian GM to fully understand the observed relationships.

Similar to other studies on the bacterial GM, MHC-I alleles appeared to have a larger influence on the fungal GM than MHC-II alleles [4, 41]. The possible mechanism underlying this pattern remains unclear, particularly as MHC-I receptors primarily interact with intracellularly derived peptides. However, it is possible that MHC-I diversity affects the fungal GM indirectly by influencing host conditions and/or altering susceptibility to other infectious agents. The MHC-I allele Ase-ua4, which was associated with a reduction in fungal ASV richness in this study, has previously been associated with differential survival in the Seychelles warbler [50]. Although the agents underlying this interaction are unknown, it is possible that they also influence the fungal GM. Another important consideration is that linkage disequilibrium between the MHC and genes that regulate other aspects of immunity and/or host physiology could also generate correlations between MHC-I alleles and the fungal GM. This could be particularly important in the Seychelles warbler, where considerable linkage disequilibrium may exist because of a recent bottleneck [138].

One caveat of our study is that MHC diversity was calculated as the number of distinct sequences detected at MHC loci. However, this number can be a poor predictor of functional diversity [139]. Indeed, several studies on wild populations have shown that the degree of sequence divergence between MHC alleles is a stronger predictor of bacterial GM dissimilarity than the number of MHC motifs [40, 41]. Thus, sequence divergence between MHC-I alleles identified in the Seychelles warbler may be more strongly associated with shifts in GM composition and should be an avenue for future research.

Overall, the effect sizes linking immunogenetic variation to GM differences were smaller than those observed in a recent investigation of the bacterial component of the Seychelles warbler GM [4]. In the latter study, three MHC-I alleles and one MHC-II allele were associated with a reduction in bacterial GM alpha diversity, and four MHC-I alleles were associated with shifts in bacterial GM composition [4]. Out of those associated with changes in bacterial GM structure, the presence of three of the same alleles (Ase-ua4, Ase-ua7, and Ase-ua11) were also found to be associated with significant, but small, reductions in fungal alpha diversity in the present study. It is possible that at least some of the observed changes in fungal GM diversity were driven by cross-kingdom microbial interactions. It is well documented that bacterial taxa can promote or antagonise the growth of fungal species (and vice versa) and so changes in one community may influence both the structure and stability of the other [108, 140]. Thus, any impact of the MHC on bacterial communities may have downstream consequences for fungal components of the GM. Indeed, a positive correlation was found between the diversity of bacterial and fungal communities in Seychelles warbler GM samples. Although establishing the mechanisms underlying these relationships is beyond the scope of the current paper, investigating the interactions between different components of the GM will be a key area for future research.

Other host factors (age and sex) were not associated with differences in fungal GM structure. This is despite the fact that age was found to be a predictor of the bacterial component of the Seychelles warbler GM [4, 9] and has also been shown to influence the fungal GM of Tibetan macaques [14]. Our study did not incorporate chicks because of a low sample size, and so, it is possible that fungal GM communities differ in very early development, before becoming largely determined by the environment. However, future studies with longitudinal data for individuals will be needed to investigate the relationship between age and fungal GM structure in more detail. Host sex was also not associated with variation in the fungal GM. As male and female Seychelles warblers exhibit low levels of behavioural and morphological dimorphism, it may be that they are exposed to, and maintain, very similar fungal communities.

Fungal GM and survival

In our study, variation in fungal GM composition (although not alpha diversity) was associated with differential survival. Six fungal ASVs were more abundant in individuals that survived to the next breeding season; however, no ASVs were significantly more abundant in those that died before the next season, suggesting that mortality may not be driven by the proliferation of pathogenic fungal species. It is likely that some of the differences in fungal composition between individuals that survived and those that died are correlated with environmental factors resulting in differences in host condition. Several of the ASVs that were more abundant in individuals that survived (e.g. those in the genus Sympodiomycopsis) were also found to be more abundant in individuals inhabiting better-quality territories in which insect food is more abundant. Thus, differences in fungal ASV abundances could be the result of differing host condition relating to factors such as food availability, stress, and/or because of differential exposure to fungal taxa in territories associated with higher or lower survival. Another possibility is that relationships between fungal taxa and survival arise due to changing microbial interactions within the GM. A previous study on the bacterial component of the Seychelles warbler showed that GM composition and the abundances of some pathogenic bacteria (e.g. Mycobacterium species) differed between adult individuals that survived or died by the next breeding season [9]. Thus, changes in fungal abundances could also represent a response to changes in bacterial communities, either due to direct microbe-microbe interactions or because abiotic conditions are altered in the intestinal tract.

Subsequent survival explained a small proportion of the overall variance in fungal GM composition (0.5% of the variation in the larger dataset). One limitation of our dataset is that some birds may have been sampled up to several months prior to death, when only small differences in the GM may have been detectable. For the majority of birds that died (25 out of 37 individuals) the point of GM sampling was the last time they were observed in the population. However, 12 individuals were seen again following GM sampling; in some cases, this was up to 7 weeks later, suggesting individuals had remained alive for a significant period post-sampling. Furthermore, there was a median period of 4 months between the last GM sample from these individuals and the next population census; individuals could have died at any point during this period. A study on captive juvenile ostriches (Struthio camelus) found that relationships between bacterial GM diversity and survival probability became stronger in the weeks directly leading up to death [10]. Thus, the fungal GM of Seychelles warblers may become more distinct closer to the point of death, as pathogenic strains proliferate and/or host condition declines. Further work, including the collection of longitudinal samples, is needed to understand if this is the case and to unpick the extent to which changes in the GM are the cause of death or if they simply represent a response to declining host health. Such studies will be crucial to understanding the extent to which the GM influences host fitness and evolutionary processes in wild populations.

Conclusion

Very few studies have investigated variation in fungal GM structure across individuals in wild animal populations, despite the known impact of fungal species on host health in captivity. In this study, we identified substantial individual variation in the fungal GM of Seychelles warblers on Cousin Island. Understanding the relative importance of environmental factors and host genetics in driving such variation is crucial if we are to understand the dynamics of the GM and its potential influence on host evolution. We show that the fungal GM of the Seychelles warbler is primarily sculpted by environmental factors. Fungal communities were particularly distinct between the major and minor breeding seasons, although it is not yet clear whether variation in climatic variables (e.g. rainfall) or reproductive activity underlie these differences. In contrast, host genetic variation had a smaller impact on the fungal GM. The presence of certain MHC alleles was associated with differences in fungal GM alpha diversity, and there was some evidence that MHC-I diversity had an impact on fungal GM composition. However, further work is needed to determine the mechanisms underlying these associations. Fungal GM composition also differed between individuals that survived and those that died by the next breeding season suggesting some link between the fungal GM and host fitness. However, further work is needed to establish the causality of such relationships, or whether changes in the fungal GM are correlated with differences in host habitat and/or condition. Defining the function of specific fungal taxa in the GM and the importance of transient strains will also play a role in establishing this. Furthermore, a holistic approach that integrates information on other microbial kingdoms in the GM and quantifies covariation between the abundances of bacterial and fungal strains will be essential if we are to fully understand GM dynamics and the importance of different GM components to host fitness.

Availability of data and materials

All ITS2 gene amplicon sequences have been submitted to the European Nucleotide Archive (ENA) database under the study accession number PRJEB54641. All 16S rRNA gene amplicon sequences have been submitted to the ENA database under the study accession numbers PRJEB45408 (samples taken in 2017 and 2018) and PRJEB47095 (samples taken in 2019). The GenBank accession number for MHC class I alleles is MZ509455-74 and for MHC class II alleles is MZ509475-98.

The scripts and metadata to reproduce all analyses and figures can be accessed via the GitHub repository, https://github.com/Seychelle-Warbler-Project.

References

  1. Sommer F, Bäckhed F. The gut microbiota–masters of host development and physiology. Nat Rev Microbiol. 2013;11:227–38.

    Article  CAS  Google Scholar 

  2. Davidson GL, Raulo A, Knowles SCL. Identifying microbiome-mediated behaviour in wild vertebrates. Trends Ecol Evol. 2020;35:972–80.

    Article  Google Scholar 

  3. Round JL, Mazmanian SK. The gut microbiota shapes intestinal immune responses during health and disease. Nat Rev Immunol. 2009;9:313–23.

    Article  CAS  Google Scholar 

  4. Davies CS, Worsley SF, Maher KH, Komdeur J, Burke T, Dugdale HL, et al. Immunogenetic variation shapes the gut microbiome in a natural vertebrate population. Microbiome. 2022;10:41.

    Article  Google Scholar 

  5. Ren T, Boutin S, Humphries MM, Dantzer B, Gorrell JC, Coltman DW, et al. Seasonal, spatial, and maternal effects on gut microbiome in wild red squirrels. Microbiome. 2017;5:163.

    Article  Google Scholar 

  6. Ren T, Grieneisen LE, Alberts SC, Archie EA, Wu M. Development, diet and dynamism: longitudinal and cross-sectional predictors of gut microbial communities in wild baboons: Gut microbiota in wild baboons. Environ Microbiol. 2016;18:1312–25.

    Article  Google Scholar 

  7. Rosshart SP, Vassallo BG, Angeletti D, Hutchinson DS, Morgan AP, Takeda K, et al. Wild mouse gut microbiota promotes host fitness and improves disease resistance. Cell. 2017;171:1015-1028.e13.

    Article  CAS  Google Scholar 

  8. Weldon L, Abolins S, Lenzi L, Bourne C, Riley EM, Viney M. The gut microbiota of wild mice. PLoS ONE. 2015;10:e0134643.

    Article  Google Scholar 

  9. Worsley SF, Davies CS, Mannarelli M-E, Hutchings MI, Komdeur J, Burke T, et al. Gut microbiome composition, not alpha diversity, is associated with survival in a natural vertebrate population. Anim Microbiome. 2021;3:84.

    Article  CAS  Google Scholar 

  10. Videvall E, Song SJ, Bensch HM, Strandh M, Engelbrecht A, Serfontein N, et al. Early-life gut dysbiosis linked to juvenile mortality in ostriches. Microbiome. 2020;8:147.

    Article  CAS  Google Scholar 

  11. Leclaire S, Pineaux M, Blanchard P, White J, Hatch SA. Microbiota composition and diversity of multiple body sites vary according to reproductive performance in a seabird. Mol Ecol. 2022;00:1–18.

  12. Comizzoli P, Power ML, Bornbusch SL, Muletz-Wolz CR. Interactions between reproductive biology and microbiomes in wild animal species. Anim Microbiome. 2021;3:87.

    Article  CAS  Google Scholar 

  13. Huffnagle GB, Noverr MC. The emerging world of the fungal microbiome. Trends in Microbiol. 2013;21:334–41.

    Article  CAS  Google Scholar 

  14. Sun B, Gu Z, Wang X, Huffman MA, Garber PA, Sheeran LK, et al. Season, age, and sex affect the fecal mycobiota of free-ranging Tibetan macaques (Macaca thibetana). Am J Primatol. 2018;80:e22880.

    Article  Google Scholar 

  15. Bergner LM, Orton RJ, Benavides JA, Becker DJ, Tello C, Biek R, et al. Demographic and environmental drivers of metagenomic viral diversity in vampire bats. Mol Ecol. 2020;29:26–39.

    Article  Google Scholar 

  16. Hicks AL, Lee KJ, Couto-Rodriguez M, Patel J, Sinha R, Guo C, et al. Gut microbiomes of wild great apes fluctuate seasonally in response to diet. Nat Commun. 2018;9:1786.

    Article  Google Scholar 

  17. Yang S, Gao X, Meng J, Zhang A, Zhou Y, Long M, et al. Metagenomic analysis of bacteria, fungi, bacteriophages, and helminths in the gut of giant pandas. Front Microbiol. 2018;9:1717.

    Article  Google Scholar 

  18. Li J, Heath IB. Chytridiomycetous gut fungi, oft overlooked contributors to herbivore digestion. Can J Microbiol. 1993;39:1003–13.

    Article  CAS  Google Scholar 

  19. Jiang TT, Shao T-Y, Ang WXG, Kinder JM, Turner LH, Pham G, et al. Commensal fungi recapitulate the protective benefits of intestinal bacteria. Cell Host Microbe. 2017;22:809-816.e4.

    Article  CAS  Google Scholar 

  20. van Tilburg BE, Pettersen VK, Gutierrez MW, Laforest-Lapointe I, Jendzjowsky NG, Cavin J-B, et al. Intestinal fungi are causally implicated in microbiome assembly and immune development in mice. Nat Commun. 2020;11:2577.

    Article  Google Scholar 

  21. Yeung F, Chen Y-H, Lin J-D, Leung JM, McCauley C, Devlin JC, et al. Altered immunity of laboratory mice in the natural environment is associated with fungal colonization. Cell Host Microbe. 2020;27:809–822.e6.

  22. Getzke F, Thiergart T, Hacquard S. Contribution of bacterial-fungal balance to plant and animal health. Curr Opin Microbiol. 2019;49:66–72.

    Article  CAS  Google Scholar 

  23. Wheeler ML, Limon JJ, Bar AS, Leal CA, Gargus M, Tang J, et al. Immunological consequences of intestinal fungal dysbiosis. Cell Host Microbe. 2016;19:865–73.

    Article  CAS  Google Scholar 

  24. Barelli C, Albanese D, Stumpf RM, Asangba A, Donati C, Rovero F, et al. The gut microbiota communities of wild arboreal and ground-feeding tropical primates are affected differently by habitat disturbance. mSystems. 2020;5:e00061–20.

  25. Montoya-Ciriaco N, Gómez-Acata S, Muñoz-Arenas LC, Dendooven L, Estrada-Torres A, Díaz de la Vega-Pérez AH, et al. Dietary effects on gut microbiota of the mesquite lizard Sceloporus grammicus (Wiegmann, 1828) across different altitudes. Microbiome. 2020;8:6.

    Article  CAS  Google Scholar 

  26. Clayton JB, Vangay P, Huang H, Ward T, Hillmann BM, Al-Ghalith GA, et al. Captivity humanizes the primate microbiome. Proc Natl Acad Sci USA. 2016;113:10376–81.

    Article  CAS  Google Scholar 

  27. Bornbusch SL, Greene LK, Rahobilalaina S, Calkins S, Rothman RS, Clarke TA, et al. Gut microbiota of ring-tailed lemurs (Lemur catta) vary across natural and captive populations and correlate with environmental microbiota. Anim Microbiome. 2022;4:29.

    Article  CAS  Google Scholar 

  28. San Juan PA, Castro I, Dhami MK. Captivity reduces diversity and shifts composition of the Brown Kiwi microbiome. Anim Microbiome. 2021;3:48.

    Article  CAS  Google Scholar 

  29. Viney M. The gut microbiota of wild rodents: challenges and opportunities. Lab Anim. 2019;53:252–8.

    Article  CAS  Google Scholar 

  30. Kreisinger J, Čížková D, Vohánka J, Piálek J. Gastrointestinal microbiota of wild and inbred individuals of two house mouse subspecies assessed using high-throughput parallel pyrosequencing. Mol Ecol. 2014;23:5048–60.

    Article  CAS  Google Scholar 

  31. Underhill DM, Iliev ID. The mycobiota: interactions between commensal fungi and the host immune system. Nat Rev Immunol. 2014;14:405–16.

    Article  CAS  Google Scholar 

  32. Radwan J, Babik W, Kaufman J, Lenz TL, Winternitz J. Advances in the evolutionary understanding of MHC polymorphism. Trends Genet. 2020;36:298–311.

    Article  CAS  Google Scholar 

  33. Klein J. Natural history of the major histocompatibility complex. New York: Wiley; 1986.

    Google Scholar 

  34. Spurgin LG, Richardson DS. How pathogens drive genetic diversity: MHC, mechanisms and misunderstandings. Proc R Soc B. 2010;277:979–88.

    Article  CAS  Google Scholar 

  35. Blum JS, Wearsch PA, Cresswell P. Pathways of antigen processing. Annu Rev Immunol. 2013;31:443–73.

    Article  CAS  Google Scholar 

  36. Hepworth MR, Monticelli LA, Fung TC, Ziegler CGK, Grunberg S, Sinha R, et al. Innate lymphoid cells regulate CD4+ T-cell responses to intestinal commensal bacteria. Nature. 2013;498:113–7.

    Article  CAS  Google Scholar 

  37. Roland MM, Mohammed AD, Kubinak JL. How MHCII signaling promotes benign host-microbiota interactions. PLoS Pathog. 2020;16:e1008558.

    Article  CAS  Google Scholar 

  38. Neefjes J, Jongsma MLM, Paul P, Bakke O. Towards a systems understanding of MHC class I and MHC class II antigen presentation. Nat Rev Immunol. 2011;11:823–36.

    Article  CAS  Google Scholar 

  39. Biedrzycka A, O’Connor E, Sebastian A, Migalska M, Radwan J, Zając T, et al. Extreme MHC class I diversity in the sedge warbler (Acrocephalus schoenobaenus); selection patterns and allelic divergence suggest that different genes have different functions. BMC Evol Biol. 2017;17:159.

    Article  Google Scholar 

  40. Bolnick DI, Snowberg LK, Caporaso JG, Lauber C, Knight R, Stutz WE. Major histocompatibility complex class IIb polymorphism influences gut microbiota composition and diversity. Mol Ecol. 2014;23:4831–45.

    Article  CAS  Google Scholar 

  41. Montero BK, Wasimuddin, Schwensow N, Gillingham MAF, Ratovonamana YR, Rakotondranary SJ, et al. Evidence of MHC class I and II influencing viral and helminth infection via the microbiome in a non-human primate. PLoS Pathog. 2021;17:e1009675.

    Article  CAS  Google Scholar 

  42. Savage AE, Zamudio KR. MHC genotypes associate with resistance to a frog-killing fungus. Proc Natl Acad Sci USA. 2011;108:16705–10.

    Article  CAS  Google Scholar 

  43. Savage AE, Zamudio KR. Adaptive tolerance to a pathogenic fungus drives major histocompatibility complex evolution in natural amphibian populations. Proc Biol Sci. 2016;283:20153115.

    Google Scholar 

  44. Belasen AM, Bletz MC, da Silva Leite D, Toledo LF, James TY. Long-term habitat fragmentation is associated with reduced MHC IIB diversity and increased infections in amphibian hosts. Front Ecol Evol. 2019;6:236.

    Article  Google Scholar 

  45. Belasen AM, Riolo MA, Bletz MC, Lyra ML, Toledo LF, James TY. Geography, host genetics, and cross-domain microbial networks structure the skin microbiota of fragmented Brazilian Atlantic forest frog populations. Ecol Evol. 2021;11:9293–307.

    Article  Google Scholar 

  46. Derakhshani H, Plaizier JC, De Buck J, Barkema HW, Khafipour E. Association of bovine major histocompatibility complex (BoLA) gene polymorphism with colostrum and milk microbiota of dairy cows during the first week of lactation. Microbiome. 2018;6:203.

    Article  Google Scholar 

  47. Komdeur J, Piersma T, Kraaijeveld K, Kraaijeveld-Smit F, Richardson DS. Why Seychelles warblers fail to recolonize nearby islands: unwilling or unable to fly there?: reduced island colonization by Seychelles Warbler. Ibis. 2004;146:298–302.

    Article  Google Scholar 

  48. Spurgin LG, Wright DJ, Velde M, Collar NJ, Komdeur J, Burke T, et al. Museum DNA reveals the demographic history of the endangered Seychelles warbler. Evol Appl. 2014;7:1134–43.

    Article  Google Scholar 

  49. Richardson DS, Westerdahl H. MHC diversity in two Acrocephalus species: the outbred great reed warbler and the inbred Seychelles warbler. Mol Ecol. 2003;12:3523–9.

    Article  CAS  Google Scholar 

  50. Brouwer L, Barr I, van de Pol M, Burke T, Komdeur J, Richardson DS. MHC-dependent survival in a wild population: evidence for hidden genetic benefits gained through extra-pair fertilizations. Mol Ecol. 2010;19:3444–55.

    Article  Google Scholar 

  51. Davies CS, Taylor MI, Hammers M, Burke T, Komdeur J, Dugdale HL, et al. Contemporary evolution of the innate immune receptor gene TLR3 in an isolated vertebrate population. Mol Ecol. 2021;30:2528–42.

  52. Gilroy DL, van Oosterhout C, Komdeur J, Richardson DS. Toll-like receptor variation in the bottlenecked population of the endangered Seychelles warbler. Anim Conserv. 2017;20:235–50.

    Article  Google Scholar 

  53. Barton GM. Viral recognition by Toll-like receptors. Semin Immunol. 2007;19:33–40.

    Article  CAS  Google Scholar 

  54. Komdeur J. Importance of habitat saturation and territory quality for evolution of cooperative breeding in the Seychelles warbler. Nature. 1992;358:493–5.

    Article  Google Scholar 

  55. Hammers M, Kingma SA, Spurgin LG, Bebbington K, Dugdale HL, Burke T, et al. Breeders that receive help age more slowly in a cooperatively breeding bird. Nat Commun. 2019;10:1301.

    Article  Google Scholar 

  56. Komdeur J, Pels MD. Rescue of the Seychelles warbler on Cousin Island, Seychelles: the role of habitat restoration. Biol Conserv. 2005;124:15–26.

    Article  Google Scholar 

  57. Richardson DS, Jury FL, Blaakmeer K, Komdeur J, Burke T. Parentage assignment and extra-group paternity in a cooperative breeder: the Seychelles warbler (Acrocephalus sechellensis). Mol Ecol. 2001;10:2263–73.

    Article  CAS  Google Scholar 

  58. Komdeur J, Daan S. Breeding in the monsoon: semi-annual reproduction in the Seychelles warbler (Acrocephalus sechellensis). J Ornithol. 2005;146:305–13.

    Article  Google Scholar 

  59. Sparks AM, Spurgin LG, Velde M, Fairfield EA, Komdeur J, Burke T, et al. Telomere heritability and parental age at conception effects in a wild avian population. Mol Ecol. 2021;31:6324–38.

  60. Brouwer L, Richardson DS, Eikenaar C, Komdeur J. The role of group size and environmental factors on survival in a cooperatively breeding tropical passerine. J Anim Ecol. 2006;75:1321–9.

    Article  Google Scholar 

  61. Griffiths R, Double MC, Orr K, Dawson RJG. A DNA test to sex most birds. Mol Ecol. 1998;7:1071–5.

    Article  CAS  Google Scholar 

  62. Coulon A. genhet: an easy-to-use R function to estimate individual heterozygosity. Mol Ecol Resour. 2010;10:167–9.

    Article  CAS  Google Scholar 

  63. R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2020.

  64. Palmer JM, Jusino MA, Banik MT, Lindner DL. Non-biological synthetic spike-in controls and the AMPtk software pipeline improve mycobiome data. PeerJ. 2018;6:e4925.

    Article  Google Scholar 

  65. D’Amore R, Ijaz UZ, Schirmer M, Kenny JG, Gregory R, Darby AC, et al. A comprehensive benchmarking study of protocols and sequencing platforms for 16S rRNA community profiling. BMC Genomics. 2016;17:55.

    Article  Google Scholar 

  66. Ihrmark K, Bödeker ITM, Cruz-Martinez K, Friberg H, Kubartova A, Schenck J, et al. New primers to amplify the fungal ITS2 region - evaluation by 454-sequencing of artificial and natural communities. FEMS Microbiol Ecol. 2012;82:666–77.

    Article  CAS  Google Scholar 

  67. White TJ, Innis MA, Gelfand DH, Sninsky JJ, editors. PCR protocols: a guide to methods and applications. San Diego: Academic Press; 1990.

    Google Scholar 

  68. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet j. 2011;17:10.

    Article  Google Scholar 

  69. Joshi NA, Fass JN. Sickle: a sliding-window, adaptive, quality-based trimming tool for FastQ files. Available at https://github.com/najoshi/sickle. 2011.

  70. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37:852–7.

    Article  CAS  Google Scholar 

  71. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.

    Article  CAS  Google Scholar 

  72. Nilsson RH, Anslan S, Bahram M, Wurzbacher C, Baldrian P, Tedersoo L. Mycobiome diversity: high-throughput sequencing and identification of fungi. Nat Rev Microbiol. 2019;17:95–109.

    Article  CAS  Google Scholar 

  73. Nilsson RH, Larsson K-H, Taylor AFS, Bengtsson-Palme J, Jeppesen TS, Schigel D, et al. The UNITE database for molecular identification of fungi: handling dark taxa and parallel taxonomic classifications. Nucleic Acids Res. 2019;47:D259–64.

    Article  CAS  Google Scholar 

  74. McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE. 2013;8:e61217.

    Article  CAS  Google Scholar 

  75. Davis NM, Proctor DM, Holmes SP, Relman DA, Callahan BJ. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome. 2018;6:226.

    Article  Google Scholar 

  76. Fredricks DN, Smith C, Meier A. Comparison of six DNA extraction methods for recovery of fungal DNA as assessed by quantitative PCR. J Clin Microbiol. 2005;43:5122–8.

    Article  CAS  Google Scholar 

  77. Hsieh TC, Ma KH, Chao A. iNEXT: an R package for rarefaction and extrapolation of species diversity (Hill numbers). Methods Ecol Evol. 2016;7:1451–6.

    Article  Google Scholar 

  78. Burnham KP, Anderson DR, Huyvaert KP. AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons. Behav Ecol Sociobiol. 2011;65:23–35.

    Article  Google Scholar 

  79. Grueber CE, Nakagawa S, Laws RJ, Jamieson IG. Multimodel inference in ecology and evolution: challenges and solutions. J Evol Biol. 2011;24:699–711.

    Article  CAS  Google Scholar 

  80. Gelman A, Yu-Sung S. arm: data analysis using regression and multilevel/hierarchical models. R package version 1.11–2. https://CRAN.R-project.org/package=arm. 2020.

  81. Fox J, Weisberg S. An R companion to applied regression, third edition. Thousand Oaks CA: Sage; 2019. https://socialsciences.mcmaster.ca/jfox/Books/Companion/.

  82. Hartig F. DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. R package version 031. 2019. http://florianhartig.github.io/DHARMa/.

  83. Barton K. MuMIn: multi-model inference. R package version 1.44.3/r480. https://R-Forge.R-project.org/projects/mumin/. 2021.

  84. Venables WN, Ripley BD. Modern applied statistics with S. 4th ed. New York: Springer; 2002.

    Book  Google Scholar 

  85. Nowak MA, Tarczy-Hornoch K, Austyn JM. The optimal number of major histocompatibility complex molecules in an individual. Proc Natl Acad Sci USA. 1992;89:10896–9.

    Article  CAS  Google Scholar 

  86. Lahti L, Shetty S. microbiome. R package version 1.14.0. 2012. http://microbiome.github.io.

  87. Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ. Microbiome datasets are compositional: and this is not optional. Front Microbiol. 2017;8:2224.

    Article  Google Scholar 

  88. Okansen J, Guillaume Blanchet F, Friendly M, Kindt R, Legendre P, McGlinn D, et al. vegan: community ecology package. R package version 2.5–7. https://CRAN.R-project.org/package=vegan. 2020.

  89. Anderson MJ. A new method for non-parametric multivariate analysis of variance: non-parametric MANOVA for ecology. Austral Ecol. 2001;26:32–46.

    Google Scholar 

  90. Lin H, Peddada SD. Analysis of compositions of microbiomes with bias correction. Nat Commun. 2020;11:3514.

    Article  CAS  Google Scholar 

  91. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995;57:289–300.

    Google Scholar 

  92. Harrison XA, McDevitt AD, Dunn JC, Griffiths SM, Benvenuto C, Birtles R, et al. Fungal microbiomes are determined by host phylogeny and exhibit widespread associations with the bacterial microbiome. Proc R Soc B. 2021;288:20210552.

    Article  CAS  Google Scholar 

  93. Nash AK, Auchtung TA, Wong MC, Smith DP, Gesell JR, Ross MC, et al. The gut mycobiome of the Human Microbiome Project healthy cohort. Microbiome. 2017;5:153.

    Article  Google Scholar 

  94. Sun B, Xia Y, Garber PA, Amato KR, Gomez A, Xu X, et al. Captivity is associated with gut mycobiome composition in Tibetan macaques (Macaca thibetana). Front Microbiol. 2021;12:665853.

    Article  Google Scholar 

  95. Raimondi S, Amaretti A, Gozzoli C, Simone M, Righini L, Candeliere F, et al. Longitudinal survey of fungi in the human gut: ITS profiling, phenotyping, and colonization. Front Microbiol. 2019;10:1575.

    Article  Google Scholar 

  96. Hallen-Adams HE, Kachman SD, Kim J, Legge RM, Martínez I. Fungi inhabiting the healthy human gastrointestinal tract: a diverse and dynamic community. Fungal Ecol. 2015;15:9–17.

    Article  Google Scholar 

  97. Hoffmann C, Dollive S, Grunberg S, Chen J, Li H, Wu GD, et al. Archaea and fungi of the human gut microbiome: correlations with diet and bacterial residents. PLoS ONE. 2013;8:e66019.

    Article  CAS  Google Scholar 

  98. Li J, Li L, Jiang H, Yuan L, Zhang L, Ma J, et al. Fecal bacteriome and mycobiome in bats with diverse diets in South China. Curr Microbiol. 2018;75:1352–61.

    Article  CAS  Google Scholar 

  99. Bensch K, Braun U, Groenewald JZ, Crous PW. The genus Cladosporium. Stud Mycol. 2012;72:1–401.

    Article  CAS  Google Scholar 

  100. Chen L, Zhang L, Li Z-H, Hui F-L. Sympodiomycopsis yantaiensis sp. nov., a basidiomycetous yeast isolated from insect frass. Int J Syst Evol Microbiol. 2013;63:3501–5.

    Article  CAS  Google Scholar 

  101. Wei Y-H, Liou G-Y, Liu H-Y, Lee F-L. Sympodiomycopsis kandeliae sp. nov., a basidiomycetous anamorphic fungus from mangroves, and reclassification of Sympodiomycopsis lanaiensis as Jaminaea lanaiensis comb. nov. Int J Syst Evol Microbiol. 2011;61:469–73.

    Article  CAS  Google Scholar 

  102. Song SJ, Sanders JG, Delsuc F, Metcalf J, Amato K, Taylor MW, et al. Comparative analyses of vertebrate gut microbiomes reveal convergence between birds and bats. mBio. 2020;11:e02901-19. https://doi.org/10.1128/mBio.02901-19.atom.

    Article  CAS  Google Scholar 

  103. Youngblut ND, Reischer GH, Walters W, Schuster N, Walzer C, Stalder G, et al. Host diet and evolutionary history explain different aspects of gut microbiome diversity among vertebrate clades. Nat Commun. 2019;10:2200.

    Article  Google Scholar 

  104. Bodawatta KH, Koane B, Maiah G, Sam K, Poulsen M, Jønsson KA. Species-specific but not phylosymbiotic gut microbiomes of New Guinean passerine birds are shaped by diet and flight-associated gut modifications. Proc R Soc B. 2021;288:rspb.2021.0446, 20210446.

  105. Lavrinienko A, Scholier T, Bates ST, Miller AN, Watts PC. Defining gut mycobiota for wild animals: a need for caution in assigning authentic resident fungal taxa. Anim Microbiome. 2021;3:75.

    Article  CAS  Google Scholar 

  106. Fiers WD, Gao IH, Iliev ID. Gut mycobiota under scrutiny: fungal symbionts or environmental transients? Curr Opin Microbiol. 2019;50:79–86.

    Article  CAS  Google Scholar 

  107. Zhang C, Derrien M, Levenez F, Brazeilles R, Ballal SA, Kim J, et al. Ecological robustness of the gut microbiota in response to ingestion of transient food-borne microbes. ISME J. 2016;10:2235–45.

    Article  Google Scholar 

  108. Santus W, Devlin JR, Behnsen J. Crossing kingdoms: how the mycobiota and fungal-bacterial interactions impact host health and disease. Infect Immun. 2021;89:e00648-e720.

    Article  CAS  Google Scholar 

  109. Tedersoo L, Lindahl B. Fungal identification biases in microbiome projects: fungal identification biases in microbiome projects. Environ Microbiol Rep. 2016;8:774–9.

    Article  Google Scholar 

  110. Tedersoo L, Anslan S, Bahram M, Põlme S, Riit T, Liiv I, et al. Shotgun metagenomes and multiple primer pair-barcode combinations of amplicons reveal biases in metabarcoding analyses of fungi. MycoKeys. 2015;10:1–43.

    Article  Google Scholar 

  111. Li S, Deng Y, Wang Z, Zhang Z, Kong X, Zhou W, et al. Exploring the accuracy of amplicon-based internal transcribed spacer markers for a fungal community. Mol Ecol Resour. 2020;20:170–84.

    Article  CAS  Google Scholar 

  112. Vargas-Gastélum L, Romero-Olivares AL, Escalante AE, Rocha-Olivares A, Brizuela C, Riquelme M. Impact of seasonal changes on fungal diversity of a semi-arid ecosystem revealed by 454 pyrosequencing. FEMS Microbiol Ecol. 2015;91:fiv044.

  113. Murillo T, Schneider D, Fichtel C, Daniel R. Dietary shifts and social interactions drive temporal fluctuations of the gut microbiome from wild redfronted lemurs. ISME Commun. 2022;2:3.

    Article  Google Scholar 

  114. Komdeur J. Seasonal timing of reproduction in a tropical bird, the Seychelles warbler: a field experiment using translocation. J Biol Rhythms. 1996;11:333–46.

    Article  CAS  Google Scholar 

  115. Escallón C, Belden LK, Moore IT. The cloacal microbiome changes with the breeding season in a wild bird. Int Org Biol. 2019;1:oby009.

    Google Scholar 

  116. Escallón C, Becker MH, Walke JB, Jensen RV, Cormier G, Belden LK, et al. Testosterone levels are positively correlated with cloacal bacterial diversity and the relative abundance of Chlamydiae in breeding male rufous-collared sparrows. Funct Ecol. 2017;31:192–203.

    Article  Google Scholar 

  117. Chi W-C, Chen W, He C-C, Guo S-Y, Cha H-J, Tsang LM, et al. A highly diverse fungal community associated with leaves of the mangrove plant Acanthus ilicifolius var. xiamenensis revealed by isolation and metabarcoding analyses. PeerJ. 2019;7:7293.

    Article  Google Scholar 

  118. Khan SA, Hamayun M, Yoon H, Kim H-Y, Suh S-J, Hwang S-K, et al. Plant growth promotion and Penicillium citrinum. BMC Microbiol. 2008;8:231.

    Article  Google Scholar 

  119. Visagie CM, Houbraken J, Frisvad JC, Hong S-B, Klaassen CHW, Perrone G, et al. Identification and nomenclature of the genus Penicillium. Stud Mycol. 2014;78:343–71.

    Article  CAS  Google Scholar 

  120. Costello EK, Stagaman K, Dethlefsen L, Bohannan BJM, Relman DA. The application of ecological theory toward an understanding of the human microbiome. Science. 2012;336:1255–62.

    Article  CAS  Google Scholar 

  121. Uren Webster TM, Rodriguez-Barreto D, Castaldo G, Gough P, Consuegra S, Garcia de Leaniz C. Environmental plasticity and colonisation history in the Atlantic salmon microbiome: a translocation experiment. Mol Ecol. 2020;29:886–98.

    Article  Google Scholar 

  122. Brown TJ, Spurgin LG, Dugdale HL, Komdeur J, Burke T, Richardson DS. Causes and consequences of telomere lengthening in a wild vertebrate population. Mol Ecol. 2021;31:5933–45.

  123. van de Crommenacker J, Komdeur J, Burke T, Richardson DS. Spatio-temporal variation in territory quality and oxidative status: a natural experiment in the Seychelles warbler (Acrocephalus sechellensis): territory quality-related oxidative costs in a wild passerine. J Anim Ecol. 2011;80:668–80.

    Article  Google Scholar 

  124. Noguera JC, Aira M, Pérez-Losada M, Domínguez J, Velando A. Glucocorticoids modulate gastrointestinal microbiome in a wild bird. R Soc open sci. 2018;5:171743.

    Article  Google Scholar 

  125. Knutie SA. Food supplementation affects gut microbiota and immunological resistance to parasites in a wild bird species. J Appl Ecol. 2020;57:536–47.

    Article  CAS  Google Scholar 

  126. Madden AA, Oliverio AM, Kearns PJ, Henley JB, Fierer N, Starks PTB, et al. Chronic stress and captivity alter the cloacal microbiome of a wild songbird. J Exp Biol. 2022;225:jeb243176.

    Article  Google Scholar 

  127. Richardson DS, Komdeur J, Burke T. Inbreeding in the Seychelles warbler: environment-dependent maternal effects. Evolution. 2004;58:2037–48.

    Google Scholar 

  128. Brouwer L, Komdeur J, Richardson DS. Heterozygosity-fitness correlations in a bottlenecked island species: a case study on the Seychelles warbler. Mol Ecol. 2007;16:3134–44.

    Article  CAS  Google Scholar 

  129. Bebbington K, Spurgin LG, Fairfield EA, Dugdale HL, Komdeur J, Burke T, et al. Telomere length reveals cumulative individual and transgenerational inbreeding effects in a passerine bird. Mol Ecol. 2016;25:2949–60.

    Article  CAS  Google Scholar 

  130. Levitz SM. Interactions of Toll-like receptors with fungi. Microbes Infect. 2004;6:1351–5.

    Article  CAS  Google Scholar 

  131. Khan AA, Yurkovetskiy L, O’Grady K, Pickard JM, de Pooter R, Antonopoulos DA, et al. Polymorphic immune mechanisms regulate commensal repertoire. Cell Rep. 2019;29:541-550.e4.

    Article  CAS  Google Scholar 

  132. Gaigher A, Burri R, San-Jose LM, Roulin A, Fumagalli L. Lack of statistical power as a major limitation in understanding MHC-mediated immunocompetence in wild vertebrate populations. Mol Ecol. 2019;28:5115–32.

    Article  Google Scholar 

  133. Slippers B, Crous PW, Jami F, Groenewald JZ, Wingfield MJ. Diversity in the Botryosphaeriales: looking back, looking forward. Fungal Biol. 2017;121:307–21.

    Article  Google Scholar 

  134. da Silva RT, Guimarães DA, Camargo ZP, Rodrigues AM, Maceira JP, Bernardes-Engemann AR, et al. Cutaneous murine model of infection caused by Neoscytalidium dimidiatum: a preliminary study of an emerging human pathogen. Med Mycol. 2016;54:890–8.

    Article  Google Scholar 

  135. Hamad I, Raoult D, Bittar F. Repertory of eukaryotes (eukaryome) in the human gastrointestinal tract: taxonomy and detection methods. Parasite Immunol. 2016;38:12–36.

    Article  CAS  Google Scholar 

  136. Parfrey LW, Walters WA, Knight R. Microbial eukaryotes in the human microbiome: ecology, evolution, and future directions. Front Microbio. 2011;2:153.

  137. Shuai M, Fu Y, Zhong H-L, Gou W, Jiang Z, Liang Y, et al. Mapping the human gut mycobiome in middle-aged and elderly adults: multiomics insights and implications for host metabolic health. Gut. 2022;71:1812–20.

  138. Zhang W, Collins A, Gibson J, Tapper WJ, Hunt S, Deloukas P, et al. Impact of population structure, effective bottleneck time, and allele frequency on linkage disequilibrium maps. Proc Natl Acad Sci USA. 2004;101:18075–80.

    Article  CAS  Google Scholar 

  139. Lenz TL. Computational prediction of MHC ii-antigen binding supports divergent allele advantage and explains trans-species polymorphism. Evolution. 2011;65:2380–90.

    Article  Google Scholar 

  140. Tipton L, Müller CL, Kurtz ZD, Huang L, Kleerup E, Morris A, et al. Fungi stabilize connectivity in the lung and skin microbial ecosystems. Microbiome. 2018;6:12.

    Article  Google Scholar 

Download references

Acknowledgements

We would like to thank the Seychelles Bureau of Standards and the Department of Environment for providing permission to conduct fieldwork and Nature Seychelles for facilitating the fieldwork on Cousin Island. This study would not have been possible without the contribution of exceptional fieldworkers and technicians associated with the Seychelles Warbler Project. Microbiome sequencing data was generated by the Centre for Genomic Research, University of Liverpool. We thank Gavin Horsburgh and Kathryn Maher for their assistance with the MHC sequencing work and Marco van der Velde for the microsatellite genotyping. The research presented in this paper was carried out on the high-performance computing clusters supported by the Research and Specialist Computing Support service at the University of East Anglia and IT Services at the University of Sheffield.

Funding

Microbiome sequencing was supported by a Natural Environment Research Council (NERC) NBAF Pilot Scheme Grant (NBAF1092) awarded to DSR and a NERC grant (NE/S010939/1) awarded to DSR and HLD. The MHC laboratory work was supported by the NERC Biomolecular Analysis Facility at the University of Sheffield (NBAF1150). CSD was funded by a NERC PhD studentship (NERC EnvEast Doctoral Training Programme grant NE/L002582/1).

Author information

Authors and Affiliations

Authors

Contributions

The study was conceived by SFW and DSR. SFW, CSD, and DSR performed the fieldwork. SFW, CSD, and MEM conducted the microbiome laboratory work. CSD conducted the MHC laboratory work and genotyping whilst affiliated with NBAF Sheffield. SFW performed the bioinformatics and statistical analyses and drafted the manuscript with input from DSR. DSR, HLD, and JK managed the Seychelles Warbler Project. The authors read and approved the final manuscript.

Corresponding authors

Correspondence to Sarah F. Worsley or David S. Richardson.

Ethics declarations

Ethics approval and consent to participate

Fieldwork was carried out in accordance with local ethical regulations and agreements. The Seychelles Department of Environment and the Seychelles Bureau of Standards approved the fieldwork.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

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 licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence 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 licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Worsley, S.F., Davies, C.S., Mannarelli, ME. et al. Assessing the causes and consequences of gut mycobiome variation in a wild population of the Seychelles warbler. Microbiome 10, 242 (2022). https://doi.org/10.1186/s40168-022-01432-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40168-022-01432-7

Keywords

  • Gut microbiome
  • Fungi
  • Major histocompatibility complex
  • Genetic variation
  • Fitness
  • Acrocephalus sechellensis