Population bottlenecks constrain host microbiome diversity and genetic variation impeding fitness.

It is becoming increasingly clear that microbial symbionts influence key aspects of their host's fitness, and vice versa. This may fundamentally change our thinking about how microbes and hosts interact in influencing fitness and adaptation to changing environments. Here we explore how reductions in population size commonly experienced by threatened species influence microbiome diversity. Consequences of such reductions are normally interpreted in terms of a loss of genetic variation, increased inbreeding and associated inbreeding depression. However, fitness effects of population bottlenecks might also be mediated through microbiome diversity, such as through loss of functionally important microbes. Here we utilise 50 Drosophila melanogaster lines with different histories of population bottlenecks to explore these questions. The lines were phenotyped for egg-to-adult viability and their genomes sequenced to estimate genetic variation. The bacterial 16S rRNA gene was amplified in these lines to investigate microbial diversity. We found that 1) host population bottlenecks constrained microbiome richness and diversity, 2) core microbiomes of hosts with low genetic variation were constituted from subsets of microbiomes found in flies with higher genetic variation, 3) both microbiome diversity and host genetic variation contributed to host population fitness, 4) connectivity and robustness of bacterial networks was low in the inbred lines regardless of host genetic variation, 5) reduced microbial diversity was associated with weaker evolutionary responses of hosts in stressful environments, and 6) these effects were unrelated to Wolbachia density. These findings suggest that population bottlenecks reduce hologenomic variation (combined host and microbial genetic variation). Thus, while the current biodiversity crisis focuses on population sizes and genetic variation of eukaryotes, an additional focal point should be the microbial diversity carried by the eukaryotes, which in turn may influence host fitness and adaptability with consequences for the persistence of populations.

Introduction

It is becoming increasingly clear that most eukaryotes live in intimate and complex relationships with microbial communities both in their external environment and on or within their body including their gut [1]. Numerous studies have documented how the presence and abundance of certain microbes can influence key aspects of host fitness, such as lifespan, fecundity, immune responses, metabolic health, behaviour, and thermal stress tolerance traits [2-8]. Conversely, the host can also control microbial composition to some extent, such as through changing nutrient availability by diet choice or host metabolism [9,10], triggering immune factors [8,11], or controlling the gut mechanically such as through peristalsis [12]. The genetic background of the host can also interact with the microbiome [13,14]. For instance, host genes can affect the abundance of certain bacteria, allowing the microbial composition of hosts to be treated as a quantitative trait in genetic analysis [15]. These interactions between the host and its microbiota can in turn have a substantial impact on host fitness [16-18].

Because interactions between microbes and their hosts shape so many aspects of life, the impact on core biological processes including evolutionary adaptation may need to be reconsidered [19-22]. The ‘holobiont’ concept reflects the idea that eukaryote individuals do not act as autonomous units, but rather as networks consisting of the host and all its associated microbiota, and that their collective genomes–the ‘hologenome’–forms a cohesive unit of selection [23-25]. Some experimental support for this idea is emerging. For instance, host selection for thermal tolerance resulted in an altered microbial composition and modulated the microbes’ response to temperature [26]. While resident microbes respond to stressful environmental conditions, they can also subsequently aid the response of the host to such conditions [7]. More generally, the microbial community may play a so-far underappreciated role in the broader context of population persistence which is important for research areas like conservation biology [18,19,27,28].

These conjectures raise the issue of how microbes respond to changes in the host’s population size. Populations that undergo repeated or persistent reductions in size frequently suffer from inbreeding depression and genetic drift resulting in lowered fitness and reduced genetic variations, ultimately impeding evolutionary capacity and increasing the risk of extinction [29-34]. These patterns have been investigated from the perspective of the nuclear genetic background rather than the hologenome; for instance, inbreeding depression is normally assumed to reflect the expression of deleterious alleles in their homozygotic form. However, host population size decreases may also lead to microbial perturbations, with potential feed-back effects on host fitness which could contribute to additional lowering of population size [27,28,35]. These interacting effects of microbial diversity, genetic diversity and fitness have not yet been widely tested.

Here we investigate microbial composition and abundance in 50 populations of Drosophila melanogaster that have experienced a variable number of generations of population bottlenecks. These populations were subsequently phenotyped for fitness components and sequenced to obtain a genome-wide molecular measure of genetic variation [33]. This set of lines provides a unique resource to investigate associations between host genetic variation, host fitness and microbial diversity (). We consider two main questions. 1) How is microbiome richness and diversity associated with host genetic variation imposed by population bottlenecks? 2) Are both host microbiome diversity and genetic variation associated with host fitness, and if so, do they contribute independently? Answers to these questions increase our understanding of the importance of microbiomes in conservation genetics and evolutionary biology. We hypothesise that reduced host genetic variation provides a less diverse environment for microbes leading to reduced microbial diversity through fewer host genotype-microbial taxon-specific interactions (). Evidence suggests that increased microbiome diversity is adaptive for the host and that increased hologenome genetic variation facilitates faster evolutionary responses to selection [22,25,36]. Our work further highlights the importance of considering both host and microbiome diversity for predicting evolutionary consequences of population bottlenecks in small and threatened populations.

Hypothetical associations between host genetic variation, microbiome diversity and host fitness.

Conceptual illustration of hypothetical associations between host genetic variation, microbiome diversity and host fitness in the experimental set of lines. Host genomic variation was manipulated through exposing populations to a variable number of bottlenecks which is assumed to have a causal effect on fitness and microbial diversity. The dark blue arrows represent potential unidirectional effects, and light blue arrows represent potential bidirectional effects. In one line of causality, low genetic variation in the host constrains the diversity (species richness and abundance) of the host-associated microbiome, which in turn affects host fitness directly, i.e. through supporting fewer host genotype-microbial taxon-specific interactions that are functionally important, and/or indirectly through changing the composition and relative abundance of beneficial or pathogenic microbes. In another line of causality, low host genetic variation directly affects host fitness which in turn leads to a less diverse microbiome through physiological or metabolic changes in the host “environment”, or indirectly such as through a higher abundance of detrimental pathogenic microbes resulting from a weakened immune response. Sources: Wikimedia Commons (Drosophila: Togopic, DNA: Kadumago), and authors’ own illustrations (microbiota).

Results

We selected 50 lines for microbiome analysis, which represent a subset of 109 inbred lines that have been through marked reductions in population size for a varying number of generations as well as 10 outbred control lines kept at a large population size (for details see Ørsted et al. 2019 [33]). These original total 119 lines had been sequenced using Genotyping-By-Sequencing (GBS) to allow nucleotide diversity (π) to be estimated. In addition, the fitness component egg-to-adult viability (hereafter just viability) was determined in each line (see Ørsted et al. 2019 [33] for details on genomic variation and phenotypes assessed). We found a positive correlation between viability and π for these 119 lines (rs = 0.45; p < 0.001; Fig A in ) indicating inbreeding depression for viability which is used here as a proxy for Darwinian fitness. For the microbiome characterization, two groups each consisting of 25 lines were selected based on measures of standardized viability and π (sum of Z-scores, see Materials and Methods). One group of 25 lines had low genetic variation and low viability, and another group of 25 lines had high genetic variation and high viability. The ‘high’ group included nine of the outbred controls as the outbred lines generally had very high viability and nucleotide variability (). The inbreeding procedure resulted in varying degrees of nucleotide diversity in the resulting lines, and as such the genetic variation of the outbred lines and the inbred lines in the ‘high’ group did not significantly differ (Fig B in ). However, to distinguish between effects of genetic variation within inbred lines from the effects of bottlenecks itself (inbreeding/outbreeding), we separated lines into three categories: ‘low genetic variation’ and ‘high genetic variation’ inbred lines and ‘outbred’ (OB) control lines (Fig B in ). Six lines evenly distributed among the genetic variation groups were removed due to a high relative abundance of the microbial endosymbiont Wolbachia, which might affect host fitness and/or abundance of other microbial taxa, resulting in 44 lines being analysed (see Materials and Methods for details).

Loss of host genomic variation decreases microbiome richness and diversity

Loss of host genetic variation was associated with decreased microbiome richness and diversity. The microbiomes from the outbred flies and the high genetic variation lines had a higher level of community diversity than the microbiomes of flies from the low genetic variation group (). Generally, we observed an increase in microbiome diversity with increasing fly host genetic variation and viability, regardless of which measure we used to quantify the microbiome diversity. Richness indices, namely observed amplified sequence variant (ASV) richness and chao1 estimated ASV richness (), showed a stronger trend than the diversity indices accounting for relative abundances (Shannon-Wiener and Simpson’s; ). Even at the class level, the microbiomes from the low genetic variation group were less taxonomically heterogeneous than those from the high genetic variation and outbred groups (Fig C in ). Most bacteria in the D. melanogaster microbiomes belonged to the Alphaprotoebacteria, Bacilli, and to a lesser extent Bacteroidia and Actinobacteria. The increasing level of taxonomic diversity with increasing host genetic variation was further evident at the genus level, with the low genetic variation group harbouring the fewest number of genera, while the outbred flies had the highest number (Fig D in ). Non-metric multidimensional scaling (NMDS) analysis based on the weighted and unweighted UniFrac distances () showed that the microbiomes differed depending on the fly host’s level of genetic variation. Generally, there was better segregation of the fly groups when only the microbiome membership (unweighted UniFrac) was considered (, with a clear trend across the viability and nucleotide diversity gradient across the first dimension of ordination space. Meanwhile the relative abundance patterns of the ASVs across the lines created noise, notably among the high and outbred groups (. However, the community membership was more similar among the outbred lines and more dispersed within the inbred lines (), while the abundance patterns showed greater similarity among the low genetic variation flies than the higher variation and outbred flies, due to the high relative abundance of a few ASVs in common in the low variation flies (). PERMANOVA analysis showed that the fly groupings explained 22.8% and 43.0% of the microbial community membership and structure (i.e. membership + relative abundance using weighted UniFrac), respectively (p < 0.001). The host traits, nucleotide diversity and viability were significantly associated with the microbiome community patterns (significant envfit arrows, ). Viability explained a greater proportion of the microbiome composition than did nucleotide diversity, where nucleotide diversity and viability explained 15.9% and 17.8% of the community membership, while explaining 27.7% and 35.3% of the community structure.

Host population bottlenecks decrease microbiome diversity.

Boxplots of microbiome diversity in the three host genetic variation groups; low genetic variation (Low; red), high genetic variation (High; green), and outbred (blue) measured as either richness indices (A. observed ASV richness (Alpha richness), and B. estimated ASV richness (chao1 estimation)) or measured using indices accounting for individual ASV abundances as well (C. Shannon-Wiener Index, and D. Simpson’s Index). In all panels, the p values of a Kruskal-Wallis test show a significant effect of group, while asterisks denote the results of pairwise Wilcox’s t-tests between groups; *** p < 0.001; ** p < 0.01; * p < 0.05; and ns: p > 0.05.

Microbiome composition differ depending on the host’s level of genetic variation.

Non-metric multidimensional scaling (NMDS) plots based on the UniFrac distances (weighted; A-B, and unweighted; C-D) between the 44 D. melanogaster lines (low genetic variation (red), high genetic variation (green), and outbred (blue). UniFrac distances account for the relative relatedness of community members, where weighted unifrac incorporates the abundance of observed organisms, while unweighted unifrac only considers presence or absence. Isolines of associated covariates are shown for nucleotide diversity (A and C) and egg-to-adult viability (B and D). Envfit values of significant drivers (p < 0.05) are shown as arrows (viability and nucleotide diversity, respectively).

Both host genetic variation and microbiome diversity contribute to host fitness

To identify drivers of host fitness, we fitted generalized linear mixed effects models (GLMMs) with viability as a function of host nucleotide diversity (π) and microbiome diversity and their interaction (one model for each of the four different measures of microbiome diversity; see Materials and Methods). For all microbiome diversity metrics, both π and microbiome diversity contributed to host fitness, with π contributing to a greater extent (~2x larger scaled effect sizes; ). For Alpha richness, we observed a significant positive interaction between host genetic variation and microbial diversity on host fitness, meaning that they did not act independently. In fact, the positive interaction coefficients suggested synergism; i.e. in lines with high genetic variation, the effect of increasing microbiome richness was greater than in lines with low genetic variation and vice versa (Fig E in ), while for chao1 ASV richness, Shannon-Wiener and Simpson’s indices, host genetic variation and microbial diversity contributed independently to host fitness (Table A in ). In all cases, the full models had higher explanatory power than each individual model (π or microbiome diversity alone), based on χ2 tests ( ). Using data from Ørsted et al. (2019) [33] for the set of lines used in the present study, we also associated microbial diversity with evolutionary responses in two traits (dry body mass and productivity measured as eggs per female per day) here defined as the slope of an ordinary linear regression across 10 generations of rearing on a stressful medium. Interestingly, we found an effect of microbial diversity on evolutionary responses for both traits (Fig F in ). For details on assessment of these phenotypes, see Ørsted et al. (2019) [33].

Host genetic variation is associated with differential bacterial relative abundance

The subset of the bacterial taxa that were differentially more abundant in the outbred flies were more diverse than in the inbred lines (i.e. the low and high genetic variation groups), with as many as 54 ASVs belonging to five bacterial classes and 12 orders (). In comparison, only one Acetobacter ASV was significantly more abundant in the bottlenecked flies, i.e. the low and high genetic variation groups. More generally, at both the ASV and genus levels, Acetobacter relative abundances were particularly high in the low genetic variation lines compared to the lines with higher genetic variation (), with Acetobacter ASVs making up the majority of the communities in the low genetic variation lines (median 94.0% for ASVs 1+2 and 96.9% for the genus). In contrast, five ASVs belonging to Lactobacillus and Enterococcus were significantly more abundant in the high genetic variation and outbred flies, while at the genus level, Enterococcus but not Lactobacillus showed this same trend ( and ). In addition, certain ASVs displayed incremental increases in relative abundances with increasing level of genetic variation and outbreeding, despite DESeq results being non-significant (). Since relative abundances are used throughout the present study, these differential abundances do not correspond to simple changes in diversity, but actual changes in certain ASVs and genera shifting between genetic variation groups.

Host genetic variation is associated with differential abundance of microbial taxa.

Relative abundance of differentially enriched bacterial ASVs and genera A-C. Acetobacter ASVs and genus, D-E. Lactobacillus ASVs and genus, and F-G. Corynebacteriales ASVs in the three host genetic variation groups; low genetic variation (Low), high genetic variation (High), and outbred. The abundance axes are not scaled the same for all of the ASVs and genera because Acetobacter constitute the majority of the microbiome of the low genetic variation group. Letters denote significant differences in relative abundance of ASVs or genera between fly groups (DESeq2 Wald test adjusted p < 0.05). The percent reads are relative to the number to which all lines were initially rarefied (14,488 reads/sample).

Core microbiomes of the low genetic variation lines are subsets of the high genetic variation lines

The number of ASVs that were present in at least 80% of lines was defined as the core microbiome. The number of ASVs belonging to this core decreased with decreasing host genetic diversity in the fly groups, with 48, 22, and 11 ASVs in outbred, high, and low genetic variation flies respectively (). Eleven of these ASVs persisted as part of the core microbiome across all of the fly groups, and belonged to the Acetobacter, Enterococcus, Lactobacillus, Leuconostoc, and Bacilli. Meanwhile, nine ASVs belonging to the Corynebacterium, Empedobacter, Nocardiaceae, Enterococcus, Mesorhizobium, and Lactobacillus, were only found in the core microbiomes of the high and outbred flies. The core microbiome of the lower genetic variation lines were subsets of those present in the higher genetic variation lines and outbreds, with 11/11 low genetic variation line ASVs represented in the high genetic variation lines, and 20/22 high genetic variation line ASVs represented in the outbred lines ().

Core microbiomes of the low genetic variation lines are subsets of the high genetic variation lines.

Venn-diagram showing that the core microbiomes of the low genetic variation lines (Low; red) are subsets of those present in the higher genetic variation lines (High; green) and outbred lines (blue), with 11/11 low genetic variation group ASVs represented in the high genetic variation group, and 20/22 high genetic variation group ASVs represented in the outbred group. The taxonomy of individual ASVs can be seen in .

Bacterial co-abundance networks decreased in complexity with host inbreeding

In parallel to the DESeq analysis, which highlighted ASVs and genera, which were categorically enriched in the different fly lines, co-abundance network analysis across the microbiome dataset was performed to identify bacterial groups that were co-varying in relative abundance with the host fitness traits of the lines and co-varying with each other (). Fly viability and nucleotide diversity co-varied with a large cluster of bacterial ASVs within the same modular cluster (). Generally, viability co-varied with a larger group of bacterial ASVs than nucleotide diversity, giving 40 and 13 associations, respectively, while 12 of these ASVs co-varied with both viability and nucleotide diversity. Interestingly, most of the ASVs that co-varied with viability were either differentially enriched in the outbred flies or in both the outbred and high genetic variation flies (). The ASVs that were differentially enriched in both outbred and high genetic variation flies belonged to 24 different taxonomic groups ( ). In contrast, the two Acetobacter ASVs that were significantly enriched in the low genetic variation flies belonged to separate modular clusters as highly prevalent negatively correlated ASVs (Fig G in ). The smaller modules were often exclusively made up of related ASVs, notably modules with Lactobacillus and Enterococcus species. The networks of each individual fly group () revealed that inbreeding the flies resulted in a lower number of co-varying ASVs, and decreased the degree of connectivity (i.e. the number of significant correlations) and overall complexity of the network, with the latter consisting of only a few ASVs in a modular structure. Having a higher host genetic variation did not result in a higher degree of covariation among ASVs among inbred lines. Host genetic variation was both positively and negatively correlated with the abundance of bacteria in the outbred network, while host genetic variation did not significantly co-vary with bacterial abundance in the high and low genetic variation networks (at correlation adjusted p < 0.05). Moreover, the ASVs belonging to the Acetobacter, that were significantly enriched in the low genetic variation flies, displayed covariation only with Lactobacillus among the inbred flies, (Fig 6E), in contrast to those ASVs in the outbred flies.

Fig 6

Bacterial co-abundance networks decreases in complexity with lower host genetic variation.

Co-abundance networks of the fly microbiome. ASVs present in at least 30 reads in total and in at least three lines; correlations > 0.5 (or <-0.5) and fdr-corrected p-values < 0.05 are shown. The nodes are individual ASVs and the host fitness traits, egg-to-adult viability (Viability) and nucleotide diversity (Nucl. diversity), while the edges represent positive and negative correlations, and correlations linking host fitness traits and bacterial ASVs (which were positive correlations). The network containing lines from all of the fly groups (A, B), the outbred (C), high genetic variation (D), and low genetic variation (E) groups are shown. In A, C, D, E, node sizes display the degree of connectivity (relative within each figure) and colors mark the modularity structure (i.e. in which community, or cluster, the ASVs belong to based on the Leiden algorithm). The orange and green edges connecting Viability and Nucleotide diversity were positive except for two instances in C where the negative correlations are displayed with dotted lines. In B, node sizes display the number of lines within which the ASVs are present, and colors mark the major groups in which the ASVs are differentially enriched. The taxonomic assignment of each node can be seen in Fig G in .

Discussion

In this study, we investigated the effects of host genetic variation on the microbiome diversity in D. melanogaster lines with experimentally manipulated levels of genetic variation. Consistent with a growing body of literature [2-8], we show that the microbial community of the host is strongly linked to the fitness of the host, in our case the viability of offspring. In addition, we link reductions in host genetic variation with microbiome diversity and clearly show that microbiome diversity is lower in lines with reduced genetic variation, regardless of whether we use richness or relative abundance metrics (). Only a few previous studies have investigated the association between components of the microbiota and host genetics in invertebrates, and those that did have mainly focused on genetic variability in the host’s ability to control the presence and/or abundance of specific endosymbionts, e.g. in [37-40]. Conversely, a few studies have investigated how endosymbionts can affect genetic variation in the host by distorting sex-ratios, such as through endosymbiont induction of cytoplasmic incompatibility or parthenogenesis [41], potentially decreasing the effective population size and increasing genetic drift [42,43] or by influencing population dynamics and dispersal possibly affecting host genetic variation [44,45].

Despite recent calls for an integration of microbiome research in evolutionary and conservation biology [18,21,22,27,28,36,46,47], little progress has been made on experimentally testing fundamental associations between population size, host genetic variation and microbial diversity. Here we provide a novel demonstration that restricting host population size has a profound impact on both the richness and the diversity of host-associated microbiomes, and that effects of host genetic variation and microbial richness interact synergistically in affecting host fitness. We also provide evidence that increased microbial diversity in the lines is associated with a stronger evolutionary response to stressful environments compared to responses seen in lines with less microbial diversity (Fig F in ). Thus, we show that populations with high host genomic variation harbor the most diverse microbial community, and vice versa for populations with low genetic variation, and that these lines are the least resilient to evolutionary pressure. According to the hologenome hypothesis, the holobiont constitutes a single unit of selection, thus our results suggest that small and fragmented populations have a reduced potential for responding to selection pressures not only due to reduced genetic variation and high rates of inbreeding but also due to reduced microbiome variation, effectively reducing hologenomic variation. Maintaining a diverse microbiome community may therefore be crucial for a host’s responsiveness to environmental change [22-24,42].

Our study sheds light on the ‘inheritance’ of microbiomes and on hosts as ‘environments’ for the microbiota. While many animals, including a range of insects, have transovarial vertical transmission of bacterial symbionts, i.e. via the egg to the embryo [1,48], the microbiome of captive D. melanogaster, especially the gut bacteria, are mainly horizontally acquired [49], and consist mostly of microbes from the diet and from microbes expelled into the immediate environment by conspecifics or predecessors [50]. Interestingly, our results show that the ‘high genetic variation’ bottlenecked populations had a markedly reduced diversity of microbes, low numbers of co-varying bacteria, and lower fitness compared to the outbred control group, despite there being no difference in genetic variation between the two groups. This could indicate that a reduction in host population size during a population bottleneck constrains the amount of microbial diversity available for random ‘sampling’ for the next generation, similar to the effects of genetic drift on host genetic variation (see ref [33] for an elaborated discussion of sampling of alleles during experimental bottlenecks). Indeed, the larger between-line variation in the microbiome community membership, i.e., unweighted UniFrac-based ordination, among the inbred flies as compared to the outbred group strongly suggests a community shift driven by repeated population bottlenecks (). Similar to drift, these effects could be exacerbated in small populations, but unlike genetic drift, which is a stochastic process, the effects of bottlenecks yield seemingly directional and predictable effects on microbiome diversity.

The trend of decreasing microbiome diversity with increasing host inbreeding has also been found in inbred populations of Diannan small-ear pigs and tortoises [51,52]. This supports that a decreased ability of inbred hosts to harbor diverse and stable microbiomes could be a general pattern, although these studies did not link microbial diversity and genetic variation to fitness effects or evolutionary responses. Previous experiments on Drosophila have highlighted the impact of fly and larval density on microbial communities and how this can impact fitness, producing a type of Allee effect [53]. Simultaneously, a reduction in fitness due to the genetic effects of population bottlenecks could in turn result in flies that cannot harbor some components of the bacterial community, which could be especially critical if microbes with specific metabolic functions are lost [18]. We are unable to tease apart these effects in the current study.

Despite this, we show that while host genetic variation has the strongest association with fitness, the microbiome diversity also contributes synergistically to host effects. We envision that increased homozygosity in hosts in a population provides a less optimal habitat for maintaining stable and diverse microbial communities. In this respect, a number of studies have previously reported that hosts which are unhealthy and/or ill, harbor different and less diverse microbiomes compared to their healthier counterparts [16,54-57]. Moreover, host genotypes have been found to profoundly affect the microbiome composition, as well as their stability and beneficial effects on the host [58-60], while other studies have shown that microbiomes affect host fitness, development, and even allele frequency and host evolution [25,61-63]. Given that flies from inbred lines with high nucleotide diversity in our study harbored less diverse microbiomes and had lower fitness than flies from outbred lines which had similar nucleotide diversity, we suggest that the microbiome is important in reducing fly fitness at high levels of genetic variation typical in natural populations. On the other hand, lower nucleotide diversity correlated with even lower microbial diversity and host fitness, suggesting that, when inbreeding results in low host genetic variation, microbial diversity loss is primarily caused by reduced host fitness. Perhaps the less optimal host ‘habitat’ for microbes in lines with low genetic diversity results in less favorable growth conditions. It has been known for a long time that productivity of Drosophila sp. can increase when there is high genetic heterogeneity in cultures, presumably because multiple genotypes can better utilize different aspects of the environment [64,65]. The genetically diverse host populations could provide more variable environments for microbes, allowing a greater diversity of microbes to persist, or support greater microbiome diversity through a higher number of functionally important host genotype-microbial taxon-specific interactions. Tests of these ideas in the context of small populations with restricted genetic variation require further studies that might examine the impact of e.g. antibiotics on host and microbe associations with fitness. Other possibilities for future research in this area involve the transfer of beneficial microbiomes to see if these can “rescue” low fitness populations or individuals suffering from disease. Such experiments have been suggested in other contexts to improve robustness and productivity of livestock and agricultural cultivars and for treating human diseases [4,66]. It remains unexplored whether such fitness improvements by microbiome transfer can alleviate the negative impacts commonly associated with low genetic variation like inbreeding depression.

The outbred lines investigated in our study have a much more complex network of co-varying, and interacting microbiomes compared to the low and high genetic variation groups (). This suggests that outbred lines harbour healthier and more robust microbial ‘ecosystems’ than bottlenecked lines, where functional redundancy within the diverse microbial community of outbred flies promotes community stability and subsequent host health. This means that removal of individual microbial species is expected to have less impact compared to in low and high genetic variation lines, where the removal of one key species may cause more dramatic effects. The importance of functional redundancy in host-associated microbiome stability and host health has been previously demonstrated in various hosts as well as in many other ecosystems [67,68]. Notably, it is now widely acknowledged that taxonomically diverse microbiomes harbor robust and stable functional redundancy, where disturbance in the environment is countered, up to a point, by the resilience built upon functional redundancy and high taxonomic diversity in a community [67,68]. The global co-variation network in our flies suggests that microbial diversity shifts gradually across a continuum of host genetic variation. However, co-variation results from individual fly lines with different levels of genetic diversity more clearly show that the tipping point between community resilience and functional redundancy and communities suffering from more stochastic microbe-microbe associations, and hence reduced functional redundancy, relates more to bottleneck treatment effects than genetic variation within bottleneck treatment. This interpretation is affected by whether lines with high relative abundance of Wolbachia was included or not, as Wolbachia can affect both host fitness and the presence/abundances of other microbial taxa [40,69-71] (see also Fig H in ).

This bears some resemblance to ecological patterns observed in microbiomes of hosts living in disturbed habitats [72-74]. For instance, howler monkeys inhabiting degraded, more homogenous forests has a less diverse gut microbiome compared to conspecifics in non-fragmented forests, and as a result has lost the microbial metabolic pathway to detoxify plant compounds in the leaves of their diet [73]. In the same way in which diversity confers resilience in macro-ecological systems [75,76], where species-rich communities are less susceptible to invasion as more niches are occupied and limiting resources are used more efficiently, such processes could be important for the robustness of microbial ecosystems within hosts [77,78]. A high diversity of commensal microbes could mean functional redundancy [18,36], thus an increased resilience towards shifts in functional diversity, and better protection against pathogens [20,79].

In conclusion, we demonstrate that restricting host population sizes and thereby genetic variation is associated with a reduction in diversity of host-associated microbiomes. We observe effects of both host genetic variation and microbial diversity in explaining host fitness, however the patterns of causality remain unclear, i.e. we suggest that it is not an all-or-nothing effect, but rather a continuum across nucleotide diversity where fitness and microbiome depletion become relatively more important drivers. It is similarly difficult to establish whether fitness differences are due to beneficial or potential pathogenic bacteria, partly because pathogenicity depends on many factors including composition of the whole microbial community present [40], and age of the host [80]. Despite this, we show clear effects of host population bottlenecks on the diversity of the microbiome, similar to effects normally observed in sampling of genetic alleles during genetic drift, and we show that both the reduction in microbiome diversity and genetic variation is tightly linked to host fitness and evolutionary potential. Therefore, microbial diversity in e.g. whole insect, fecal or environmental samples could be a useful proxy for population fitness and potentially adaptability. Lastly, our results open up multiple avenues for further studies, such as transplantation of microbiomes as a means of ‘rescuing’ populations that suffer from inbreeding depression, potentially relevant to species of conservation concern.

Materials and methods

Fly population bottlenecks

The D. melanogaster flies used in the study originated from 232 wild females caught at Oakridge winery in the Yarra Valley, Victoria, Australia. These females contributed equally to a mass bred population kept at a population size of approximately 1000 individuals. The flies were maintained on a 12:12 h light:dark cycle at 19°C on a standard laboratory food composed of yeast, sucrose, oatmeal, and agar mixed with tap water. Nipagen and acetic acid were added to prevent fungal growth (Table B in ). From the mass bred population, we created 40 lines of each of three different expected levels of inbreeding for a total of 120 lines, from which we could obtain nucleotide diversity measures for 109 lines. This is described in detail in Ørsted et al. (2019) [33]. In summary, lines from the three levels of inbreeding experienced 2, 3, and 5 consecutive generations of bottlenecks each consisting of two males and two females (census size, N = 4) resulting in expected inbreeding coefficients of F = 0.125, 0.219 and 0.381, respectively. After reaching the desired inbreeding level, we flushed the population sizes to 200 individuals. Simultaneously, we established 10 control lines that were assumed outbred and maintained at minimum 500 individuals.

Fly lines selected for microbiome analysis

We obtained a measure of the realised genetic variation in all inbred and outbred lines using genotyping-by-sequencing (GBS) to calculate nucleotide diversity (π) from genomic SNPs (described in detail in Ørsted et al. (2019) [33]; average number of SNPs ± sd = 26,877 ± 4,061). For the microbiome analysis, we aimed at comparing two groups of fly lines: one group with low genetic variation and overall low performance and one with high genetic variation and high overall performance. For a quantitative selection of these groups, we calculated a composite measure of performance as the sum of standardised values of viability and nucleotide diversity (Z-score). We selected the 25 lowest ranking lines and the 25 highest ranking lines based on a summed Z-score (50 lines total; Fig A in ). Because we selected lines regardless of their inbred/outbred status, nine of the ten outbred control lines were included in the ‘high genetic variation’ lines due to their high Z-scores sum. In all analysis, we differentiate between three groups of flies: ‘low’ and ‘high’ genetic variation lines within the inbred lines and ‘outbred’ (OB) lines that have not been through population bottlenecks. Following establishing the lines, population sizes were flushed to ca. 200 individuals per line, and in the F1 generation egg-to-adult viability was assessed by distributing 15 eggs into each of five vials per line and calculating the proportion of eggs developing successfully to the adult stage [33]. Flies used for microbiome analysis consisted of male flies emerging from the egg-to-adult viability assessment (merged across the 5 vials) that were snap frozen in liquid nitrogen at 2–3 days of age and subsequently stored at -80°C.

16S rRNA gene sample preparation and sequencing

The whole fly genomic DNA was extracted following the same protocol as has been described in Ørsted et al. (2019) [33]. In brief, 15 randomly collected male flies from each experimental line was homogenized by bead beating at 2x6 s cycles and 6500 rpm using a Precellys mechanical homogenizer (Bertin Techologies, Montigny le Bretonneux, France), and the DNA was purified with the DNeasy Blood & Tissue kit (QIAGEN, Hilden, Germany) with modifications specific for extracting insect tissues. The V1-V3 hypervariable regions of the bacterial 16S rRNA gene was amplified and multiplexed for Illumina sequencing according to Albertsen et al. (2015) [81], and the library pool was sequenced on a MiSeq sequencer using the MiSeq Reagent kit v3, 2x300bp paired-end configuration, and 20% PhiX control spike-in. The raw demultiplexed sequenced reads have been deposited under SRA Bioproject (accession number PRJNA813140) at NCBI.

Microbiome analysis

The demultiplexed paired-end reads were quality filtered, assembled to make consensus amplicon reads, clustered into ASVs, and taxonomically assigned using a custom workflow AmpProc version 5.1.0.beta2.11.1 (https://github.com/eyashiro/AmpProc), which relies on USEARCH version 11.0.66_i86linux64 [82] sequenced reads processing and FastTree version 2.1.10 [83] for tree building. When present in insects, the endosymbiont Wolbachia can affect host fitness and/or the presence/abundance of other microbial taxa, complicating interpretation of results [40,69-71]. Therefore, we removed six lines with a relative Wolbachia abundance >85%. These six lines were evenly distributed among the genetic variation groups (Fig I in ). This threshold was chosen to maintain a relatively high number of reads per line/sample, i.e., at least 14,000 reads per sample, after Wolbachia reads were removed. Prior to further analysis, we also tested whether removing Wolbachia from the frequency table would skew our results and interpretations. There was no difference in relative abundance of Wolbachia between genetic variation groups (Kruskal-Wallis Rank Sum test; χ2 = 2.54, df = 2, p = 0.28). Further, we found no significant correlations between Wolbachia abundance and nucleotide diversity (rs = 0.24, t48 = 1.68, p = 0.10) or fitness (rs = 0.23, t48 = 1.67, p = 0.10), and there was no significant difference between the six removed lines and the 44 remaining lines in mean nucleotide diversity (Wilcoxon rank sum test; W = 110, p = 0.53) or mean fitness (W = 116, p = 0.64). Thus, our main results namely that lines with low genetic variation and fitness harbor a less diverse microbiome were not affected by removing lines high in Wolbachia abundance. However, the prevalence of Wolbachia reads above 10% per sample distorted the relative abundance of other bacterial taxonomic groups to a small extent (Fig I in ), thereby leading to slightly different interpretations of the co-abundance networks (see Fig H in for details).

Next, QIIME version 1.9.1 [84] was used to rarefy all of the samples to the smallest acceptable sample size i.e. 14,488 reads per sample (single_rarefaction.py), and to generate the observed and estimated chao1 richness and Shannon-Wiener and Simpsons diversity values for each sample (alpha_diversity.py), weighted and unweighted UniFrac matrices (beta_diversity.py), and core microbiome groups (compute_core_microbiome.py). R version 3.5.0 was used for downstream analysis and statistics. To assess the composition of core microbiomes, we first normalized the dataset to eight lines per fly group (which is the number of outbred lines in the study after filtering samples with low read counts and removing Wolbachia) by randomly keeping three, two, and three lines from the high, medium, and low inbreeding levels, respectively, from the high and low genetic variation groups [33], and all of the outbred lines. Next, we considered an ASV as belonging to the core microbiome, the percentage to the nearest tenth if we were to accept a prevalence of an ASV in 80% of the lines in each fly group.

Non-metric multidimensional plots based on the UniFrac matrices were generated with the R-package ‘vegan’ (version 2.5–2 [85]). The effect of two host traits (nucleotide diversity, egg-to-adult viability; [33]) were calculated with the envfit function and the significant envfit values (p < 0.05) of the fly host traits were displayed as arrows and fly host trait gradients across the ordination space as isolines using the ordisurf function. PERMANOVA analysis using vegan’s adonis function was undertaken for the UniFrac matrices as a function of fly lines, nucleotide diversity and viability. Differentially enriched ASVs were identified with DESeq2 (version 1.22.2 [86]) using the Wald test, contrasting across three levels (low and high genetic variation and outbred fly groups), and accepting the enriched ASVs with adjusted p < 0.05. The same approach was used on the genus-level frequency table. To test for differences in measures of microbiome diversity between the genetic variation groups, we performed Kruskal-Wallis tests and pairwise Wilcox’s t-tests between groups.

Generalised linear mixed models

To assess drivers of fitness measured as egg-to-adult viability, we fitted generalised linear mixed effect models (GLMMs) in the R-package ‘lme4’ [87] assuming a binomial distribution with logit link function. We fitted viability as a function of nucleotide diversity and microbiome diversity (four different measures as described above) and their interaction as fixed effects; replicate vials were included as a random effect, as flies from the same vial are not independent. The full models were compared with individual univariate models of either nucleotide diversity or microbiome diversity by χ2 difference tests. We detected no over-dispersion in any of the models (residuals to df ratio > 0.31; χ2 > 76.02; p > 0.99). Conditional coefficients of determination of the GLMMs interpreted as the variance explained by the entire model, including both fixed and random effects, were calculated as , where are the variances of the fixed effect components, the random effects and the residual variance, respectively (see ‘delta-method’ in Nakagawa et al. (2017) [88]).

Co-occurrence network analysis

For co-occurrence network analysis, the ASV table of all of the fly lines was used for the global network, and the subset of ASV tables generated for the core microbiomes analysis for outbred, high genetic variation, and low genetic variation groups were used as the starting frequency tables. In each of these ASV frequency tables, the fly nucleotide diversity and egg-to-adult viability measurements were included. The ASVs present in fewer than three lines and ASVs with less than 30 reads in the respective frequency tables were removed. The ‘Hmisc’ R-package [89] was used to generate Spearman correlation matrices, and the upper triangle was used to generate the network edge list. Only edges with FDR-corrected p < 0.05 and correlation coefficients > 0.5 (or < -0.5) were kept. In Gephi version 0.9.2, the network was generated with the Fruchterman Reingold algorithm and minor crossing of edges were corrected manually to improve the visualization of the clusters. From among Gephi’s built-in features, the Leiden algorithm was used to calculate the modularity structure using edge weights and default settings, and the degree of connectivity was calculated as the number of nodes with which each node is connected by an edge.

Supporting figures and tables.

(DOCX)

Click here for additional data file.

ASV table, nucleotide diversity, fitness and microbiome diversity of fly lines.

S1 Table found online as a Microsoft Excel spreadsheet. Nucleotide diversity, fitness and microbiome diversity (observed alpha richness, estimated chao1 richness, Shannon-Wiener and Simpson’s indices), and ASV table rarefied to 14,488 reads per sample of the 44 fly lines used in the main results, as well as ASV table rarefied to 29,000 reads per sample where Wolbachia was not prefiltered, and which includes the six lines that were removed from the original 50 lines due to high relative Wolbachia abundance (in a separate tab). Headers of the metadata sheets are described in the last tab.

(XLSX)

Click here for additional data file.

DESeq2 table showing differential abundance of microbial taxa.

DESeq2 table is found online as a Microsoft Excel spreadsheet. The log fold change calculations are done on the ASV (L7 tabs in the spreadsheet) and genus level frequency tables (L6 tabs in the spreadsheet). The DESeq contrasts were done among outbred, high and low genetic variation fly groups, and the pairwise contrasts are displayed in separate sheets. For the most significant pairwise comparisons (p < 0.1), the fly group in which the ASVs or genera are enriched is indicated.

(XLSX)

Click here for additional data file.

Table 1

Host genetic variation and microbiome richness interact synergistically on host fitness.

Microbiome diversity measure	Fixed effects	Estimate	Std. Error	z value	p
Observed ASV Richness (ObsAlpha)	Intercept	0.106	0.081	1.311	0.1899
R2full model = 0.355	NuclDiv	0.784	0.088	8.883	6.49E-19***
	ObsAlpha	0.499	0.087	5.740	9.49E-09*
	NuclDiv*ObsAlpha	0.194	0.089	2.182	0.0291*
	Random effects	Std.Dev
	Vial	0.683
	Comparison with individual models	χ2 df	χ 2	p
	NuclDiv; ObsAlpha; Full model	2	70.322	5.37E-16	***

Results of the general linear mixed model (GLMMs) of egg-to-adult viability as a function of nucleotide diversity (NuclDiv) and microbiome diversity (Alpha richness; ObsAlpha) and their interaction as fixed effects. Both dependent and independent variables are scaled (Z-standardization) to allow direct comparison of effect sizes. Replicate vial IDs were included as a random effect, as flies from the same vial are not considered independent. Conditional coefficients of determination of the GLMMs interpreted as the variance explained by the entire model, including both fixed and random effects, are shown. Asterisks denote the significance of individual variables or interactions; *** p < 0.001; ** p < 0.01; and * p < 0.05. The full model including both dependent variables and their interaction is compared with individual models with either nucleotide diversity or alpha diversity with a χ2 test.

Table 2

Core microbiomes of the three host genetic variations groups.

ASVs	Class	Order	Family	Genus	Species	Low	High	OB
ASV1	Alphaproteobacteria	Acetobacterales	Acetobacteraceae	Acetobacter	sp.	+	+	+
ASV2	Alphaproteobacteria	Acetobacterales	Acetobacteraceae	Acetobacter	sp.	+	+	+
ASV4	Alphaproteobacteria	Acetobacterales	Acetobacteraceae	Acetobacter	sp.	+	+	+
ASV5	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus	sp.	+	+	+
ASV7	Bacilli	Lactobacillales	Lactobacillaceae	Lactobacillus	plantarum	+	+	+
ASV8	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus	sp.	-	-	+
ASV9	Bacilli	Lactobacillales	Leuconostocaceae	Leuconostoc	sp.	+	+	+
ASV10	Actinobacteria	Corynebacteriales	Corynebacteriaceae	Corynebacterium	sp.	-	+	+
ASV11	Bacteroidia	Flavobacteriales	Weeksellaceae	Empedobacter	sp.	-	+	+
ASV12	Bacilli	Lactobacillales	Lactobacillaceae	Lactobacillus	sp.	+	+	+
ASV13	Actinobacteria	Corynebacteriales	Nocardiaceae	g__	sp.	-	+	+
ASV14	Actinobacteria	Corynebacteriales	f__	g__	sp.	-	+	+
ASV15	Bacilli	Lactobacillales	Lactobacillaceae	Lactobacillus	sp.	+	+	+
ASV16	Bacilli	Lactobacillales	Lactobacillaceae	Lactobacillus	plantarum	+	+	+
ASV17	Actinobacteria	Corynebacteriales	Nocardiaceae	g__	sp.	-	+	+
ASV18	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus	sp.	-	+	+
ASV19	Actinobacteria	Corynebacteriales	Corynebacteriaceae	Corynebacterium	sp.	-	-	+
ASV21	Alphaproteobacteria	Caulobacterales	Caulobacteraceae	Brevundimonas	sp.	-	-	+
ASV22	Bacteroidia	Sphingobacteriales	Sphingobacteriaceae	Sphingobacterium	mizutaii	-	-	+
ASV23	Gammaproteobacteria	Burkholderiales	Alcaligenaceae	Alcaligenes	sp.	-	-	+
ASV24	Bacteroidia	Sphingobacteriales	Sphingobacteriaceae	Sphingobacterium	mizutaii	-	-	+
ASV25	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Mesorhizobium	sp.	-	+	+
ASV27	Bacilli	Lactobacillales	Lactobacillaceae	Lactobacillus	brevis	+	+	+
ASV28	Gammaproteobacteria	Burkholderiales	Alcaligenaceae	Alcaligenes	sp.	-	-	+
ASV29	Actinobacteria	Corynebacteriales	Nocardiaceae	Rhodococcus	sp.	-	-	+
ASV30	Actinobacteria	Micrococcales	Microbacteriaceae	Leucobacter	sp.	-	-	+
ASV32	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	g__	sp.	-	-	+
ASV36	Actinobacteria	Micrococcales	Micrococcaceae	Glutamicibacter	sp.	-	-	+
ASV40	Actinobacteria	Micrococcales	Microbacteriaceae	Microbacterium	sp.	-	-	+
ASV45	Alphaproteobacteria	Caulobacterales	Caulobacteraceae	Brevundimonas	sp.	-	-	+
ASV46	Actinobacteria	Micrococcales	Brevibacteriaceae	Brevibacterium	Brevibacterium	-	-	+
ASV48	Bacilli	Bacillales	Planococcaceae	Lysinibacillus	sp.	-	-	+
ASV53	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus	sp.	-	+	+
ASV56	Bacilli	Lactobacillales	Lactobacillaceae	Lactobacillus	brevis	-	+	+
ASV68	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus	sp.	-	-	+
ASV72	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Mesorhizobium	sp.	-	-	+
ASV89	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus	sp.	-	-	+
ASV91	Actinobacteria	Corynebacteriales	Corynebacteriaceae	Corynebacterium	flavescens	-	-	+
ASV93	Bacilli	Bacillales	Planococcaceae	Lysinibacillus	sp.	-	-	+
ASV94	Bacilli	Lactobacillales	Lactobacillaceae	Lactobacillus	sp.	-	+	-
ASV105	Bacilli	Bacillales	Planococcaceae	Lysinibacillus	sp.	-	-	+
ASV108	Bacilli	Bacillales	Planococcaceae	Lysinibacillus	sp.	-	-	+
ASV110	Bacilli	Bacillales	Planococcaceae	Lysinibacillus	sp.	-	-	+
ASV114	Bacilli	Bacillales	Planococcaceae	Lysinibacillus	sp.	-	-	+
ASV115	Bacilli	Bacillales	Planococcaceae	Lysinibacillus	sp.	-	-	+
ASV117	Bacilli	Bacillales	Planococcaceae	Lysinibacillus	sp.	-	-	+
ASV152	Bacilli	o__	f__	g__	sp.	+	+	+
ASV189	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus	sp.	-	-	+
ASV193	c__	o__	f__	g__	sp.	-	-	+
ASV383	Bacilli	Lactobacillales	Leuconostocaceae	Leuconostoc	sp.	-	+	-

Core microbiome ASVs in each of the three host genetic variations groups: low genetic variation (Low), high genetic variation (High), and outbred (OB). The ASVs belonging to the core microbiome for each fly group is represented by + (core) and–(not core). Here we define presence in the core microbiome if an ASV is present in at least 80% of lines of a particular group. The lowest taxonomy level is listed (o__: order, f__: family, g__:genus).