Context-dependent dynamics lead to the assembly of functionally distinct pitcher-plant microbiomes

Niche construction through interspecific interactions can condition future community states on past ones. However, the extent to which such history dependency can steer communities towards functionally different states remains a subject of active debate. Using bacterial communities collected from wild pitchers of the carnivorous pitcher plant, Sarracenia purpurea, we tested the effects of history on composition and function across communities assembled in synthetic pitcher plant microcosms. We found that the diversity of assembled communities was determined by the diversity of the system at early, pre-assembly stages. Species composition was also contingent on early community states, not only because of differences in the species pool, but also because the same species had different dynamics in different community contexts. Importantly, compositional differences were proportional to differences in function, as profiles of resource use were strongly correlated with composition, despite convergence in respiration rates. Early differences in community structure can thus propagate to mature communities, conditioning their functional repertoire.


Introduction 33
Microbes profoundly shape our ecosystems, yet we still lack a clear understanding of the 34 processes driving community assembly and related ecosystem functioning 1

. Community 35
assembly is difficult to predict because the dynamics of any particular species can be dependent 36 on the community context; niches are created or destroyed through biotic interactions with other 37 members of an assembling community 2-4 . Microbial communities are complex, with many 38 species and different kinds of interactions among species. For example, microbes can facilitate 39 other species' growth via excretion of metabolic waste products 5-8 , or actively interfere with 40 their growth through the production antimicrobial compounds 9 . Microbes can also engage in 41 strong cooperative interactions, whereby energy transducing metabolic interactions are coupled 42 across species 10 . This diversity of interactions creates many potential contexts for species 43 dynamics, implying that the behavior of one species is dependent on the background of 44 interactions. Stochastic changes in the biotic context-for example, priority effects and random 45 colonization or extinction events-can have long-lasting consequences for community structure 3 . 46 As a result, microbial communities can reach different compositional states due to variation in 47 biotic context in addition to variation in abiotic environmental conditions such as weather events 48 or resource pulses. These history-dependent effects are collectively called "historical 49 contingencies" 2-4 . 50 The extent to which historical contingency leads to alternative community states 2,11 51 remains an active subject of debate. In a strongly selective environment historical contingencies 52 may not matter such that communities converge to the same compositional and functional 53 outcomes 12 . For example, bacterial communities from widely different sources converged 54 reproducibly in single-carbon source synthetic media 7 . Beyond microbes, Mediterranean plant 55 communities in environments with frequent fires were formed by related groups of species that 56 share key traits 13 . In contrast, other studies have found that historical contingency leads to 57 different community compositions and functions, for example: priority effects led to large 58 differences in ecosystem function for wood-decaying fungi 14 and productivity of grassland 59 plants 15 . A third, and perhaps largest, set of studies has found convergence in terms of function 60 but not species (or phylogenetic) composition, for example: grassland plants 16 , the stratified 61 layers of a hypersaline microbial mat 17 , microbial communities colonizing the surface of 62 seaweed 18 , and the bacteria and archaea living in bromeliad tanks 19 . 63 Functional convergence without species convergence is more likely when the functions 64 being measured are performed by many species from different lineages, and thus are redundant 65 within the broader species pool. Here, again, there are contrasting reports in the literature on the 66 prevalence of functional redundancy. A number of studies have found functional redundancy in 67 microbial communities 19,20 , while others emphasize important functional differences that depend 68 on species composition 21 . The degree of redundancy is clearly related to the function and system 69 examined; for example, aerobic respiration is found across many bacteria, while the ability to 70 degrade lignin is rare. Thus, 'narrow' functions, such as the hydrolysis of complex carbon 71 compounds, are generally carried out by rare community members 22 . When relevant functions 72 are variable across genetic backgrounds (not highly redundant) and dependent on interactions, 73 historical contingencies could have major effects on the functional capabilities of a community 74 and on nutrient cycling within ecosystems 14 . Therefore, studies in microbial ecology need to 75 address more specific, relevant functional measurements and examine not just on historically 76 contingent compositional states, but also historically contingent functional states. 77 microcosm communities, ranging from about 6 to 16 (small colored numbers in Figure 1b). In 124 summary, communities tended to first converge due to the fast loss of diversity (primarily along 125 the first NMDS axis), but remain distinct due to historical contingencies (primarily along the 126 second NMDS axis). 127 Figure 1. Microcosm communities approach distinct equilibria. a) Relative abundances of the top 20 ASVs across the ten microcosms during the course of the serial transfer experiment. ASVs are listed one time each on the bar plot, with taxonomic classification in the legend. The DNA concentrations for each timepoint are graphed above as points, with a Loess fit as a solid line. b) Communities change quickly and then stabilize in the Non-metric Multidimensional Scaling (NMDS) plot of Bray-Curtis dissimilarities of the microcosm community compositions. The microcosm name is listed in black next to the Day 0 community of each microcosm, and the lines connect the timepoints. Colored numbers indicate the mean effective number of species for each community post-Day 21. c) The Bray-Curtis dissimilarity of ASV relative abundances between adjacent days decreases over the course of the experiment. The thick black line shows a Loess fit to all data points, and the thin line marks Day 21.
Interestingly, the differences in richness near-equilibrium were seemingly pre-determined 129 at early stages of assembly. Communities lost many ASVs between days 0-3 during initial 130 adjustment to the laboratory environment (Figure 2a), but the richness measured from the first 131 timepoint of the experiment (Day 3) was an excellent predictor of the richness at the end of the 132 experiment (Day 63) with R 2 = 0.9008 and p < 0.0001 (Figure 2b). Notably, richness in samples 133 taken directly from the pitcher plant (Day 0) had no significant correlation with that of Day 63 134 (R 2 = 0.1978, p = 0.1105), consistent with the notion that some ASVs in the pitcher plant fluid 135 that were either metabolically inactive or unable to grow in our experimental conditions. Thus, 136 changes in community composition after only three days of adjusting to the lab environment 137 propagated throughout the assembly process, suggesting that historical contingencies played a 138 significant role in structuring communities. 139  Given the strong correlation between initial and final richness, we wondered whether 141 temporal dynamics of community equilibration were also correlated across microcosms. For example, the temporal dynamics of ASV loss could be driven by external factors such as transfer 143 intervals and the dilution factor, in which case all microcosms should exhibit comparable 144 dynamics. Alternatively, the dynamics might be driven by community context and thus specific 145 to each microcosm. To answer this question, we measured the distribution of ASV extinction 146 times across the microcosm communities (Fig. 2c). As a null model, we tested if the loss of 147 ASVs is a random process where each ASV that is bound to go extinct in a given community 148 context does so with a fixed probability p in each transfer. This assumption implies that the 149 extinction time distribution is described by a geometric distribution, which indeed gave a good 150 fit for all microcosms (Pearson's chi-squared test for differences was not significant: p > 0.05 in 151 all cases). To investigate if microcosms can be described by a common extinction time 152 distribution, we used the Akaike Information Criterion to compare two models: using either one 153 parameter per microcosm or a single parameter describing all microcosms. Surprisingly, the 154 single-parameter model was strongly favored (relative likelihood of ~1000) suggesting that 155 temporal dynamics of ASV extinction in our communities are the same and likely driven by 156 external factors rather than biotic context. This 'universal' dynamic of species loss is also 157 revealed when studying relative richness, i.e. normalized by richness at Day 3 ( Figure 2d). After 158 this normalization, all relative richness curves mapped well to our model's maximum likelihood 159 distribution and approached a common equilibrium relative richness level (approximately ~50% 160 of the initial richness, regardless of its absolute value, Figure 2d). Thus, about half of all initial 161 ASVs were doomed to go extinct at a random point during the assembly process, while the rest 162 persisted indefinitely. 163 We investigated if community context influenced the dynamics of individual ASVs. 164 Although the ten microcosms reached different compositional states, on average ~90% of the communities were composed of ASVs that occurred in more than one microcosm 166 (Supplementary Figure S1), and this overlap allowed us to ask to what extent the same species 167 behave similarly or not when their community context changes. To this end, we first looked at 168 the extinction times and found that the majority of the shared ASVs dropped out of different 169 microcosms at different times (Figure 3a), often with large differences in their persistence time. 170 For example, ASV681 (in Figure 3b) dropped out by Day 12 in microcosm M08, but persisted at 171 high relative abundance through the end of the experiment in microcosms M04, M06 and M10. 172 We focused on ASVs that persisted at least once past the 7 th transfer (Day 21) and found that 173 although the same ASV can act differently in different community contexts (Figure 3b), the 174 temporal trajectories of the same ASVs in different microcosms tended to be somewhat more 175 correlated with each other than randomly chosen ASVs (Fig. 3c). This result implies that, as one 176 might expect, species identity determines species abundance dynamics in different biotic 177 contexts to some degree. However, only 55% of these ASVs had significantly correlated We found that the effects of historical contingency on community assembly were 186 remarkably consistent and reproducible. We performed the same experiment on a second set of 187 microcosms from the same inocula that were not subjected to initial filtering, but otherwise 188 underwent the same transfer protocol and amplicon sequencing. However, despite likely  . Community composition strongly correlates with functional activity. a) Variance in percent CO2 among the ten microcosms over the course of the serial transfer experiment. b) NMDS plot of the Bray-Curtis dissimilarities of functional activity (substrate use) as measured by EcoPlates. c) Procrustes rotation of the composition NMDS plot with the function NMDS plot for Day 63 (Procrustes correlation = 0.8991, p < 0.001). d) Plot of Bray-Curtis dissimilarities for composition by function (Mantel test r = 0.6403, p < 0.001), squared to better illustrate the spread of points. e) Endochitinase activity over time for the five microcosms that strains were cultured from (top row) and frequency of the ASVs over time that map to strains with measurable endochitinase activity (bottom row). Endochitinase activity is shown in units/mL for the strains/ASVs in the bottom row by a gradient from gray (low activity) to red (high activity).
In contrast with CO2 production and the expectation of functional convergence, the 208  To ask to what extent the pattern of substrate utilization of the community could be 251 reduced to the substrate utilization of its members, we applied the same Biolog EcoPlates test to 252 the individual isolates. We found that substrate utilization differences between communities 253 could be attributed to differences in species abundance. Out of the 20 strains measured using 254 Biolog EcoPlates, only one showed high metabolic activity when grown with salicylic acid 255 (strain M09D5GC17 corresponding to ASV863 Burkholderia), and it was a strain present only in 256 M09 which was the only microcosm where growth on salicylic acid was observed. Conversely, 257 the strain with the most growth on itaconic acid (strain M10D5GC19 corresponding to ASV140 258 Achromobacter) was present at the final timepoint in all the microcosms we cultured from except 259 for M09, and the M09 community was the only one not able to grow on itaconic acid 260 (Supplementary Figure S5). Individual strains thus drive key functional differences among 261 microcosms, and we were able to culture and analyze a set of these strains-connecting 262 genotypes with their functional phenotypes. 263

265
Discussion 266 The effects of historical contingencies are difficult to isolate and detect in the field, 267 because even adjacent sites can experience distinct environmental conditions. The effects can be difficult to capture in a laboratory as well, because historical contingencies may only affect 269 community assembly when disturbance is low and environmental selection is weak 2,12,35 and 270 these conditions are rarely met when natural communities are moved into laboratory settings. 271 Our study demonstrates that historical contingency strongly influences community assembly in a 272 realistic, but controlled, laboratory environment. Moreover, the effects of historical contingency 273 are persistent and reproducible, suggesting that, with enough information about species' 274 functional capabilities, responses to environmental conditions, and interactions with other 275 species, community assembly dynamics might be predictable. natural systems because species that did not grow within the current environment were pruned during transfers, but are more complex than almost all experimental laboratory communities. The 314 ability to culture key community members provides the opportunity for building synthetic 315 communities that retain interactions among species previously established in nature. 316

Sample collection 319
We collected the entire aquatic pools from 10 healthy pitchers of Sarracenia purpurea 320 their native environment, we used cricket media (3 grams cricket powder from farmed Acheta 331 domestica crickets per 1 Liter of milliQ-purified water, acidified to pH 5.6 and then autoclaved). 332 The plate was then placed in a 25ºC incubator. After three days of incubation, each sample was 333 mixed well and 500 µL was transferred to a new plate with 500 µL of sterile cricket media. From 334 our calculations, we added the equivalent of about 1/60 th of a cricket to every well at each 335 transfer. We continued transferring samples and adding cricket media in a 1:1 ratio every three 336 days for a total of 21 plates over 63 days. 337 At the beginning of the experiment (Day 0), we removed a portion of each sample to 338 freeze at -80ºC for later DNA extraction and amplicon sequencing. We removed 100 µL of each 339 filtered sample to first measure optical density (OD) at 600 nm, and then used a Fluorimetric Chitinase Assay Kit (Sigma-Aldrich) to measure the activity of three different types of 341 chitinases: endochitinases, chitobiosidases and β-N-acetylglucosaminidases. We bead-beat each 342 sample for 1 minute and centrifuged it before using a portion of the supernatant in the assay, with 343 two replicates for each sample. Our downstream analyses focused on endochitinases, the 344 enzymes that cleave intramolecular bonds forming new chain ends. 345 We also measured a "functional fingerprint" of the communities with Biolog EcoPlates 37 . 346 We diluted each sample 1:10, combining 1 mL of each filtered sample with 9 mL of purified 347 water that had been acidified to pH 5.6 and then autoclaved. We filled each EcoPlate well with 348 100 uL of sample, and incubated the plates in the 25ºC incubator for three days. At the end of 349 this time the plate was read in a plate reader according to EcoPlate instructions. 350 In addition to measuring chitinase and EcoPlate activity, we measured CO2 production 351 with the MicroResp system 38 . We added 250 µL of each filtered sample to 250 µL of cricket Amplicon Sequence Variants (ASVs) of ~250 base pairs in length. We retained all ASVs with 382 more than two sequences across all samples. We assigned taxonomy using the classify-sklearn 383 method which is a Naive Bayes classifier, and a pre-trained classifier made with the Greengenes 384 database, version 13_8.
Statistical analyses were performed and graphs were generated in R and Mathematica. 386 Reads for each ASV were normalized by the total amount of reads in each sample. Bray-Curtis 387 dissimilarities and NMDS ordinations for Figure 1 were performed using the R vegan package 41 . 388 The effective number of species was calculated as exp(Shannon index). Richness for Figure 2  389 was calculated counting all ASVs present at later time points in each microcosm as present at 390 previous time points to account for ASVs below the sequencing detection limit. Richness of 391 early timepoints (Days 0 and 3)  To describe the dynamics of ASV loss, we employed a geometric model, i.e., the 398 probability ( ) of an ASV that goes extinct eventually to go extinct at transfer is equal to the 399 probability that it did not go extinct in the preceding − 1 transfer. That is, 400 Where is the sole parameter of the model. It can be shown that the maximum likelihood 402 estimator of is given by the mean time to extinction, i.e., ̂= 1 〈 〉 0 . To determine whether the 403 extinction time distributions are best described by a single parameter or whether individual 404 parameters for each microcosm are needed, we computed the Akaike Information Criterion 405 from the likelihood ℒ of the observed extinction under the geometric model, using either a single 406 parameter given by the inverse of the mean extinction time for all ASVs across all 407 microcosms, or for each microcosm individually, i.e.,

Correlation analysis 414
For the correlation analysis in Fig. 3c- shared between microcosms. We found that f = 80% of strains occurring in two microcosms had 426 the same fate in both, suggesting that a given strain has an 80% chance of having the same fate in 427 a new microcosm as it had in its current microcosm. We thus estimate the probability P(n) of a 428 similar strain occurring in n microcosms to have to same fate in all n microcosms as ( ) = 429 K*+ , which is shown as the red line in Fig. 3e. 430

Strain isolation and identification 432
Individual strains were isolated from five of the ten microcosms (M03, M05, M07, M09 433 and M10) by plating the culture fluid and picking around 100 colonies per microcosm. We 434 cultured on petri plates at 1:10 5 and 1:10 6 dilutions using both cricket media with the addition of 435 agar and a second medium containing 11.28 g/L M9 salts, vitamin solution 42 , trace metals 43 , 436 agar, and 2.5 g/L N-acetylglucosamine (GlcNAc) as the sole carbon source. Plates were put in a 437 25ºC incubator for one week, after which single colonies were picked and then re-streaked at 438 least two additional times before being grown up in liquid media and frozen in glycerol at -80ºC. 439 To identify and barcode our strains, we added 2 µL of liquid culture to 20 µL of sterile, 440 nuclease free water and, after a freeze-thaw cycle, did direct PCR using primers 27F and 1492R 441 to amplify the 16S ribosomal RNA gene and the Q5 High-Fidelity kit (New England Biolabs) 442 with an initial 5-minute incubation at 100ºC. Before Sanger sequencing, we tested for successful 443 amplification with gel electrophoresis and cleaned the PCR products with SPRI beads according 444 to the Agencourt AMPure XP protocol. Sanger sequences were trimmed and filtered with 445 Geneious, and assigned taxonomy using the RDP classifier. 446 447

Measurements of strain functional activity 448
We measured the functional activity of approximately 50 strains with 100% matches in 449 their Sanger-sequenced 16S rRNA gene to ASVs from the amplicon sequencing. When multiple 450 strains mapped to the same ASV, their enzyme activities were averaged. Each strain was 451 streaked out on cricket-M9 media plates from the frozen glycerol stock, and then a single colony 452 was grown in liquid cricket media. Chitinase activity was measured as described for the