We used the B6.129X1-Fut2^tm1Sdo/J mouse model (Domino et al., 2001), kept in independently ventilated cages under specific pathogen free conditions and with ad libitum water and chow. The initial B6.129X1-Fut2^tm1Sdo/J (Fut2^-/-) mice were purchased from the Jackson Laboratory and interbred with wildtype C57BL/6J (Fut2^+/+) mice. Animals were intercrossed by mating a Fut2^-/- male or female (which we distinguish by “breeding direction”) with a respective wildtype mouse. The resulting heterozygous offspring were mated to obtain the experimental animals of different Fut2 genotypes and the offspring kept together according to gender. Feces were sampled first maximally 5 days after weaning (4 weeks after birth), and then every 7 days for 11 weeks and dissected after the final sampling (see Supplementary Table S1 for an overview of mice and respective breeding and sample information; TP-time points). Fecal samples were collected on ice and immediately transferred to -80°C for storage. To obtain the mucosa associated microbial communities after dissection (aged 11 weeks after weaning, “TP11”) the tissue was gently washed in 4 ml RNAlater^® (Ambion^®; Carlsbad, CA, United States) and stored separately from the luminal content in a fresh tube containing 1.5 ml RNAlater^®. Samples preserved in RNAlater^® were left over night at 4°C, were spun down, and depleted of supernatant before storage at -80°C. To avoid cross-contamination, instruments were rinsed and cleaned with 70% ethanol between anatomic sites. The approval for mouse husbandry was obtained from the local veterinary office. All animals were kept under a 12 h light-dark cycle and ad libitum water and standard chow diet. All experiments were performed following the German Regulations of Animal Welfare as approved by the “Ministerium für Landwirtschaft, Umwelt und ländliche Räume” (Kiel, Germany).

Tissue preserved in RNAlater^® was washed twice in PBS, fixed in 10% formalin for 24 h and embedded in paraffin. Five micrometer sections were deparaffinized and subjected to antigen retrieval in citrate buffer. Sections were stained with rhodamine-labeled Ulex europaeus agglutinin-1 (UEA-1; 1:100, Vector Laboratories, Burlingame, CA, United States) for 1h at room temperature, washed and mounted in Prolong Gold containing DAPI (Molecular Probes) and investigated at 100 × magnification.

DNA was extracted from feces and RNAlater^®-preserved mucosa samples (stored at -80°C) with the Qiagen Stool DNA Isolation Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. The 16S rRNA gene was amplified using forward (5′-CTATGCGCCTTGCCAGCCCGCTCAGTCAGAGTTTGATCCTGGCTCAG-3′) and reverse (5′-CGTATCGCCTCCCTCGCGCCATCAGXXXXXXXXXXCATGCTGCCTCCCGTAGGAGT-3′) primers flanking the V1 and V2 hypervariable regions (27F-338R) and was sequenced following the methods described in Rausch et al. (2011). The 454 Life Sciences primer B (forward) and A (reverse) adapter sequences are denoted in boldface, and the underlined sequences represent the broadly conserved bacterial primers 27F and 338R. A 2-base linker sequence (TC/CA; shown in italics) was added as recommended by Roche^® (454^®) (Margulies et al., 2005). A unique 10-base multiplex identifier (designated as XXXXXXXXXX) was added to the reverse primer to tag each PCR product. Raw sequence information including metadata can be accessed at the European Bioinformatics Institute ¹ under the accession number PRJEB13483.

Raw sequences were trimmed using mothur 1.31.2 (Schloss et al., 2009) requiring no ambiguous bases, a mean quality score of ≥ 35 and a minimum length of 200 nucleotides for the V1-V2 region. Chimeric sequences were detected by Usearch 4.25 (overlap of de novo and database-driven detection) (Edgar, 2010) and removed. Sequences were subsampled to 1500 sequences per sample, classified via RDP classifier (RDP 9, as provided by P. Schloss) with ≥ 60% bootstrap threshold (Cole et al., 2003; Wang et al., 2007) and aligned to the SILVA core database. Operational Taxonomic Units (OTUs) were binned in mothur using the average distance clustering. Phylogenetic tree construction on representative OTU sequences was carried out using FastTree 2.1 with a gamma corrected CAT substitution model (Price et al., 2010). With an average sequence count of 1494.45 ± 36.78 SD per sample, we achieved a high species coverage of 88.49 ± 6.15% over all samples, for fecal samples 87.70 ± 4.02% at TP1, 83.81 ± 7.32% at TP3, 86.13 ± 6.96% at TP5, and 84.26 ± 4.26% at TP11, respectively, where “TP#” stands for the number of weeks after weaning. Mucosa associated microbial communities were covered at the species level on average by 87.31 ± 3.81% for cecal tissue, 92.17 ± 3.13% for colonic tissue, 94.11 ± 2.45% for the jejunum, and 94.98 ± 2.39 for ileal tissue.

Species diversity indices (species richness, Shannon–Weaver index, as well as the phylogenetic diversities) were calculated in R 3.0.1 using the vegan and picante packages (Kembel et al., 2010; Oksanen et al., 2011; R Core Team, 2016). Phylogenetic measures were derived using species occurrences and their respective relative phylogenetic relatedness, resulting in unweighted NRI (Net Relatedness Index) and NTI (Nearest Taxon Index) (Webb, 2000). These measures represent phylogenetic effect sizes inferred by contrasting the observed relatedness patterns in a respective sample and a null model, by maintaining species occurrence frequency and sample species richness over 999 iterations (Gotelli, 2000). The phylogenetic measures of beta diversity, unweighted UniFrac, were calculated in mothur and provide insight into differentially present or abundant phylogenetic lineages between communities (Lozupone and Knight, 2005). Beta diversity metrics based on shared OTU presence (Jaccard distance), or shared abundance (Bray–Curtis distance) were calculated in the vegan package for R. For constrained ordination (Redundancy Analysis, RDA) the OTU tables were Hellinger-transformed (Rao, 1964; Legendre and Gallagher, 2001). We further employed distance-based Redundancy Analysis based on Jaccard-, UniFrac-, and Bray–Curtis distances (Legendre and Anderson, 1999). The significance of factors and axes was assessed via a permutative ANOVA approach (5000 permutations, stratified by time point or gastrointestinal (GIT) location in global analyses). All distances were square root-transformed to avoid negative eigenvalues. The analyses of community variability were performed using the betadisper function with a correction for subset size, as implemented in vegan (Anderson, 2006). Univariate analyses of repeated measurements (i.e., time course, whole GIT) was conducted using linear mixed models (LMM) with mouse ID as a random variable and a cage-dependent variance structure, and optimized by model selection using cAIC and its weights (Bartoń, 2016) (for example code see Supplementary Methods). For the analyses within a certain time point or GIT location, cage was identified as an important random variable for LMM analysis. All models were reduced under normal maximum likelihood (ML) and finally refitted using REML (Zuur et al., 2009; Pinheiro et al., 2011) (for example code see Supplementary Methods). Indicator species analysis was based on 10000 permutations using the generalized indicator value (IndVal.g) to assess the predictive value of a taxon for each respective host phenotype/category (occurrence in at least 10% of samples) (De Cáceres et al., 2010). P-values of the genus and OTU associations were deemed significant below a q-value of 0.05 (Storey et al., 2015).

Consensus genera and species-level OTU networks were generated using the SparCC algorithm (100 iterations, 10000 permutations, without singleton OTUs) implemented in mothur to avoid spurious correlations induced by compositionality (Friedman and Alm, 2012). The networks were constructed from the correlation matrices and weighted by the correlation coefficient in igraph 1.0.1 for R using only associations with P ≤ 0.005 (genera) or P ≤ 0.001 (species level OTUs) (Csardi and Nepusz, 2006). Calculation of relative centrality scores and network manipulation were also carried out in igraph, and weights were transformed to real values for the derivation of centralities (Freeman, 1979), modularity (Clauset et al., 2004; Newman and Girvan, 2004) and other network characteristics. Network robustness tests were performed by random attacks on networks by sequentially removing 25% of network nodes randomly over 1000 iterations, and mean values of the network characteristics for each fraction were used for further analysis. Targeted attacks were performed by sequentially removing the randomized 25% highest connected bacteria (highest node degree/highest number of connections) from the networks. We further simulated network attacks based on the association strength of bacteria to a host characteristic (Fut2 genotype, breeding direction), as measured by IndVal.g (generalized indicator value) (De Cáceres et al., 2010). The top 25% of associated taxa were used for sequential removal and permuted 1000 times to exclude any effect of removal order. We identified the most influential bacteria on network stability by the average real values of relative change in network characteristics after removal (diameter, number of network clusters, size of network clusters, transitivity, closeness). Random networks of similar sizes were repeatedly simulated (100 iterations) based on a similar degree distribution as the empirical network (Viger and Latapy, 2005) by randomly and evenly distributing associations between nodes (Erdős and Rényi, 1959) and by partially rewiring random networks (k = 0.6, 0.8; i.e., small world network with k as rewiring probability). The scale-free networks of similar size were constructed with a power law degree distribution (p = 1, 2, 4) (Goh et al., 2001), as were the exponential networks (p = 4, 6, 8) (Barabási and Albert, 1999). These networks were subjected to random and targeted attack as described above (100 iterations in each constructed network) and the resulting decay distribution was compared to original networks with the two-sided Kolmogorov–Smirnov test.

To explore the microbial communities at a basic level, we first investigated the dynamics of the major microbial taxonomic groups over time (feces) and among locations of the gastrointestinal tract (GIT; intestinal tissue sampled in mature mice, approx. 15 weeks of age/11 weeks after weaning). These analyses, based on LMMs considering cage effects, reveal interesting patterns among the major phyla Bacteroidetes, Firmicutes, and Proteobacteria with respect to secretor status and initial breeding direction (see Figure 1 and Table 1). Firmicutes and Bacteroidetes develop significant differences over time between secretor types depending on breeding direction, specifically toward the last time point. Proteobacteria appear to be mainly influenced by the initial breeding direction and decrease over time in mice with a Fut2^-/- grand dam, resulting in a significant difference toward the end of the time course (see Figure 1C and Table 1).

Phylum abundances within the mucosal communities show a weak decrease in Firmicutes and increase in Bacteroidetes in the ileal mucosa of Fut2^+/- mice. Proteobacteria abundances in the ileal mucosa on the other hand display an interaction between secretor type and breeding direction, with a higher abundance in secretors of the secretor-founded mouse line compared to non-secretors, and the opposite trend in the non-secretor-founded mouse line. Proteobacteria are the only bacterial group in the jejunum which show an effect of Fut2 genotype, particularly an increase in abundance with the number of functioning Fut2 alleles. In the proximal regions of the GIT (cecum, colon), the effects of the Fut2 genotype/secretor status on bacterial phyla decreases. Only Proteobacteria in the cecum are influenced by Fut2 genotype, with a lower abundance of Proteobacteria in non-secretor animals (Supplementary Table S2). The major phylum in the cecum, Bacteroidetes, is less abundant in the mucosal community of animals inoculated by Fut2^-/- grand dams (Supplementary Table S2).

Next, we further investigated the microbial communities for species and genera that are characteristic of Fut2 genotype, secretor status, or breeding direction via indicator species analysis (De Cáceres et al., 2010). The majority of significant associations is present with regard to breeding direction and include taxa such as Prevotella, Parabacteroides, Lactobacillus, and Lachnobacterium (Table 2 and Supplementary Tables S3, S4 and Figures S2, S3). However, several significant associations are also present for Fut2 genotype and/or secretor status, including Odoribacter and Propionibacterium (Fut2^+/+) and Prevotella, Lachnobacterium and Ruminococcus (non-secretors) (Table 2 and Supplementary Table S4).

Community complexity is reflected by the diversity of taxa within samples (i.e., alpha diversity) and is an important factor to understand ecological resilience and productivity. We measured the complexity of our community in complementary ways by describing the number of observed species (species richness), their distribution (Shannon entropy), and their relative phylogenetic relatedness (Net Relatedness Index- NRI, Nearest Taxon Index- NTI) (Webb, 2000). We identified an influence of maternal transmission on fecal microbial communities by all three alpha diversity measures considering potential cage effects, showing comparable starting diversities, strong deviations at subsequent time points and a final convergence of diversity at the end of the time course (Figure 2 and Table 3). Species richness decreases after weaning in animals whose grand dam was a non-secretor/Fut2^-/-, but equilibrates between breeding lines toward the end of the time course (TP11; Supplementary Figure S4 and Table S5). The orderedness of the species abundance distributions, as described by Shannon entropy, shows different trajectories over time according to breeding direction and Fut2 genotype. This results in significant differences in diversity among Fut2 genotypes at the last time point (Fut2^+/+ > Fut2^+/- > Fut2^-/-; Supplementary Figure S4 and Table S5) and evident differences of Fut2^+/+ mice between breeding directions.

With respect to the phylogenetic structure of the bacterial communities, we found stronger phylogenetic clustering in secretor- than in non-secretor mice. This phylogenetic clustering deteriorates over time until the communities become phylogenetically unstructured (secretors, NRI ∼ 0) to overdispersed in non-secretors (NRI < 0). Secretor status thus influences the trajectory of community assembly on a broad phylogenetic scale. The analysis of NRI shows that succession of bacterial communities results in phylogenetically less restricted communities (phylogenetic dispersion), specifically in non-secretors. The phylogenetic relatedness among closely related species (NTI), however, is increased in mice derived from a Fut2^+/+ grand dam, which indicates a phylogenetically restricted set of bacteria. These communities cluster quickly and remain clustered, while the bacterial communities derived from a Fut2^-/- grand dam diversify (decrease of NTI) and become more restricted (increase of NTI) in later time points. The microbial communities passed from either breeding line converge to the end of the time course, but take different trajectories during this development. However, analyzing the diversification of phylogenetic groups over time (LMM for NRI: time span- F_1,101 = 21.018, P < 0.0001, R²_marginal = 0.067; Spearman correlation: ρ = -0.3298, P = 8.80 × 10^-5) indicates colonization and the establishment of more distantly related groups (e.g., classes, orders). The corresponding increase of NTI over time (LMM for NTI: time points- F_1,101 = 23.153, P < 0.0001, R²_marginal = 0.250; Spearman correlation: ρ = 0.2967, P = 0.0005), however, indicates a decrease in phylogenetic distance among species and genera within a respective host. Mucosa-associated microbial communities at different locations along the GIT mostly do not differ according to breeding direction or genotype. The breeding lineage founded by Fut2^+/+ grand dams shows a significant decrease in species richness and entropy only in the ileum (Supplementary Table S6 and Figure S5).

In summary, community complexity significantly changes with respect to time point after weaning and host characteristics. This implies that the succession of microbial communities take different trajectories depending on the initial direction of the experimental cross, but is also influenced by Fut2 genotype within the respective breeding lineages.

To investigate differences between microbial communities we employed beta diversity measures based on differential species occurrence, abundance, and phylogenetic relatedness among samples (Jaccard, Bray–Curtis, and UniFrac, respectively). Interestingly, formal tests of community differentiation according to Fut2 genotype or secretor status revealed no significant differences between communities in the overall dataset [Redundancy Analysis (Rao, 1964), distance-based Redundancy Analysis (Legendre and Anderson, 1999)]. However, breeding direction consistently influences microbial communities in feces and the mucosa. Further, when we consider the interaction of Fut2 genotype and breeding direction we find mainly differences according to differential abundance in the fecal bacterial communities 1-, 3-, and 5-weeks after weaning (Figure 3 and Supplementary Tables S7, S8). The effect of this interaction decreases toward the end of the time course, while breeding direction remains significant (Supplementary Tables S7, S8). This differentiation, specifically between breeding directions is correlated with time (Bray–Curtis: direction- F_1,132 = 3.0753, P = 0.0002, time span- F_1,132 = 3.9970, P = 0.0002, direction × time span- F_1,132 = 1.4997, P = 0.0144, adj. R² = 0.0396; Jaccard: direction- F_1,132 = 2.1468, P = 0.0002, time span- F_1,132 = 1.8771, P = 0.0002, direction × time span- F_1,132 = 1.2623, P = 0.0046, adj. R² = 0.0167). However, the analysis of all time points together (permutations stratified by time points) reveals significant community differentiation with respect to Fut2 genotype within and between breeding directions (Table 4 and Supplementary Figure S6). Analyses of the mucosal microbial communities of single sites in the GIT reveals mainly an effect of breeding direction on the mucosal microbial communities (Supplementary Table S9), while in the cecal mucosa differences arise between secretor and non-secretor individuals (Supplementary Table S10). A combined analysis of all mucosal microbial communities reveals genotype effects in addition to interactions with breeding direction, further emphasizing the direct and indirect influence of ABH antigens in the mucosa (Table 4).

In addition to the differentiation between communities over time, we observe a decrease in overall community variability, and thus a homogenization of bacterial community composition. Differences in community variability become stronger between Fut2 genotypes and secretor status over time (TP11- secretor status: Jaccard: P = 0.0344, Bray–Curtis: P = 0.0342, UniFrac: P = 0.3002; Euclidean: P = 0.0435). These differences occur specifically between homozygous secretors (Fut2^+/+) and non-secretors (Fut2^-/-; Jaccard: P = 0.0189; Bray–Curtis: P = 0.0217; UniFrac: P = 0.1548; Euclidean: P = 0.0756; Supplementary Figure S7). This implies a less variable community composition in Fut2^-/- mice. Breeding direction, on the other hand has almost no observable effect on the variability of fecal communities. In summary, we observe strong community differentiation between breeding lines, with an additional separation of communities according to Fut2 genotype, the latter of which is limited to the mucosa.

To infer potential interactions among bacterial community members we used co-occurrence networks of bacterial consensus genera and species, based on the correlation of their abundances. To construct networks for each time point, we split the abundance table accordingly and calculated measures of determination on the same set of bacteria and animals for each stratum using a correlation procedure considering compositionality (Friedman and Alm, 2012). We next measured different bacterial network characteristics (importance/centrality of network components) to approximate the structural importance of the single bacterial taxa within each respective network, including the number of interactions of each bacterium in the network (node degree), the PageRank^TM index (Brin and Page, 1998; Allesina and Pascual, 2009), and the well-known “betweenness”- and “closeness” centralities (Freeman, 1979). In simple terms, PageRank^TM will assign a high importance to genera that interact with other important genera recursively, while betweenness centrality measures the number of interactions in which the respective genus is a mediator (on the shortest path between any pair of taxa). Closeness centrality, on the other hand, measures how distant the node is from all other nodes in the network, whereby high scores indicate a small distance toward the rest of the network members.

When we compare the networks between single time points, patterns resembling community succession and stabilization emerge (Figure 4A, for species networks see Supplememtary Figure S8A). The number of interactions between bacterial genera decreases over time (genera node degree: ρ = -0.497, P < 1.00 × 10^-15; Figure 4B), while the species node degree increases (ρ = 0.144, P < 2.20 × 10^-16; Supplementary Figure S8B). However, not only the number of connections, but also the direction of interactions changes from a higher number of negative interactions to an almost even distribution of positive and negative interactions on the genus level (positive/negative interactions; genus level: TP1: 0.292, TP3: 0.492, TP5: 0.791, TP11: 0.776; species level: TP1: 1.433, TP3: 1.170, TP5: 1.394, TP11: 1.405). Positive interactions are, however, on average stronger than negative ones (Supplementary Figure S9). Furthermore, we could identify an increase in the overall importance of single genera within these assemblages. The average mediator role of bacterial genera and species does not significantly change over time (genera betweenness: ρ = 0.0879, P = 0.0876, species betweenness: ρ = -0.0035, P = 0.8228), and the relative distance of links between community members increases with time due to reduced connectedness (ρ = –0.7643, P < 2.20 × 10^-16; see Figure 4B), while distances among species decrease over time (ρ = 0.1692, P < 2.20 × 10^-16; Supplementary Figure S8B). The PageRank^TM as a generalized importance index reveals an increase of importance with time in genera networks (ρ = 0.148, P = 2.21 × 10^-6; Figure 4B), a pattern also present within species networks (ρ = 0.127, P = 2.31 × 10^-16; Supplementary Figure S8B). All these results hint toward an increasing stabilization of a core set of strongly interacting bacteria over time, with increasing mutualistic- and less competitive interactions between single taxa, which increases modularity (see Supplementary Figure S10 for genus-based network modularity; OTU network modularity: TP1: 0.838, TP3: 0.936, TP5: 0.962, TP11: 0.960) (Clauset et al., 2004; Newman and Girvan, 2004).

To test the robustness of the microbial communities against different disturbance regimes, we simulated random loss of community members or targeted depletion of important keystone members (i.e., network hubs- based on the number of their interactions/node degree) (Albert et al., 2000). For each respective co-occurrence network we measured the resulting changes in the network such as fragmentation (number and size of connected components), size (diameter), and clustering (transitivity, closeness). As expected for ecological communities, which are mainly classified as scale-free networks, we found high resistance against random failures, but fast deterioration of communities after removing small fractions of important, highly connected nodes (see Figure 4C and Supplementary Figure S11, attack vs. random failure). We observed the highest resistance against community collapse shortly after weaning. In this stage many weak interactions seem to prevent the community from rapid deterioration (Figure 4 and Supplementary Figures S10, S11). In later time points bacterial interactions appear to solidify, while others are lost during succession. This differentiation increases the importance of single community members (see the increase of PageRank^TM above) and therefore the vulnerability of the system. Thus, the targeted removal of important mediators (targeted attack) disturbs communities much stronger than random network failures, as shown for example in the stronger network fragmentation (see Figure 4C and Supplementary Figure S8C) and size of the biggest connected component (Supplementary Figures S11A,B). The transitivity/clustering of the networks deteriorates much quicker under targeted removal compared to random removal throughout the time course (Supplementary Figures S11E,F), as does network diameter/size (Supplementary Figures S11C,D) and closeness centrality (increase of distance between nodes; Supplementary Figures S11G,H). In sum, these analyses reveal high community resistance in the early phases of community development and a high susceptibility of climax communities to the removal of key players.

The most influential bacteria (top 10 members for each time point) on network decay were relatively universal and also include indicator taxa for genotype, secretor status, and breeding direction. Several taxa appear to have common large effects on most network parameters (diameter, size of biggest cluster, transitivity, closeness) over more than a single time point, such as Alistipes (OTU-36), Rikenella (OTU-9), and Paludibacter (OTU-120) in species networks, or Alistipes, Enterorhabdus, Lachnospiraceae, Parasutterella, and Ruminococcus in the genus-level networks (Supplementary Tables S11, S12). Several taxa with a large influence on network decay are also indicators for host factors (globally or time point specific), such as Ruminococcus, which is associated to the Fut2^-/- genotype/non-secretor status (Table 2). Indicators for breeding direction are also found among structurally important bacteria including Lachnospiraceae, Barnesiella, and Bacteroides as indicators for the Fut2^+/+ grand dam, or Prevotellaceae, Parasutterella, and Escherichia/Shigella as indicators for the Fut2^-/- grand dam (Table 2 and Supplementary Table S11). Similarly, several of the most important OTUs in the species level networks are indicators of breeding direction such as Paludibacter, Helicobacter, Oscillibacter, and Anaerophaga (Fut2^+/+ grand dam), or Meniscus, Tannerella, and Bacteroides as indicators for the Fut2^-/- grand dam (Supplementary Tables S3, S12). These findings imply a relatively large structural role of bacteria influenced by breeding legacy and transgenerational genotype effects, but also some direct genotype effects (e.g., Ruminococcus association). Furthermore, the behavior of the empirical networks under random error- and targeted attack regimes is similar to simulated networks of a comparable degree distribution (small world, random, degree sequence, power law), but resembles most closely the behavior of “small world”-like networks of similar size (Supplementary Figure S12 and Table S13). Interestingly, the co-occurrence network derived from the last time point appears more resilient than any of the simulated network topologies tested (Supplementary Figure S12).

To investigate the influence of the Fut2 genotype, secretor status, and breeding direction on the robustness of microbial communities, we adapted the heuristics described above and sequentially removed the top 25% of genera and species associated to the respective host factors (as measured by their respective indicator values (De Cáceres et al., 2010) from the network. The removal of bacteria associated to secretor genotypes influences network fragmentation (i.e., number of connected components) earlier and stronger than the removal of non-secretor associated bacteria (Figure 4C and Supplementary Figures S13, S14). Especially in early community development (TP1, TP3) the average number of resulting subnetworks is higher than after removal of non-secretor associated bacteria. Secretor-dependent effects on the network characteristics are most obvious in the species networks, which also show a higher influence of secretor-associated bacteria (Figure 4 and Supplementary Figure S14). We observe similar patterns with respect to Fut2 genotype, with the strongest effect on network characteristics when removing Fut2^+/+-associated bacteria (Supplementary Figures S8D, S15, S16). Further, removal of bacteria associated to a Fut2^+/+ grand dam disturbed the network earlier than the removal of bacteria associated to the Fut2^-/- grand dam (Supplementary Figures S17, S18). However, disturbances introduced by the removal of bacteria (genera/species) associated to the Fut2^-/- grand dam are on average higher than the removal of bacteria associated to the Fut2^+/+ grand dam or random removal. Thus, community stability appears to significantly differ with respect to breeding direction in our experimental cohort of mice.

In summary, our results imply that changes in network characteristics over time resemble patterns of community succession, including an increase in taxon importance, attack vulnerability, and taxonomic/phylogenetic heterogeneity. Furthermore, we find that host genotype and legacy effects can significantly influence network stability. These results illustrate the importance of specialized key members in microbial communities and may suggest new strategies to evaluate community stability and dynamics.