prey that supplement their nutrient requirements, particularly for nitrogen (Ellison & Gotelli, 2002;Gotelli & Ellison, 2002). However, S. purpurea do not produce several classes of enzymes required to digest prey and instead rely on inquiline macroinvertebrate and microbial communities present in pitcher fluid for the mechanical and hydrolytic enzymatic breakdown of captured insects, respectively (Adlassnig et al., 2011;Bittleston et al., 2018;Gallie & Chang, 1997;Luciano & Newell, 2017;Prankevicius & Cameron, 1991;Young et al., 2018). This breakdown of insect carcasses liberates bioavailable nitrogen and other nutrients for uptake by host plant tissues, while also supporting inquiline communities (Butler & Ellison, 2007;Karagatzides et al., 2009). In this way, both top-down (macroinvertebrate) and bottom-up (microbial) controls shape availability of key resources (e.g. nitrogen) to promote ecosystem stability and function in the S. purpurea ecosystem.

One key macroinvertebrate in the S. purpurea aquatic microecosystem is the pitcher plant mosquito, W. smithii. Unlike most mosquito species that indiscriminately lay their eggs into small bodies of standing water, adult W. smithii females obligately lay their eggs into Sarracenia pitchers, where larvae progress through four aquatic juvenile stages (instars) and an aquatic pupal stage before emerging as terrestrial adults (Bradshaw & Creelman, 1984;Heard, 1994). Pitcher-associated microbes play essential roles in W. smithii larval development by serving as cues for growth and moulting (Arellano & Coon, 2022), as has been documented in other mosquito species (Coon et al., 2014(Coon et al., , 2016a(Coon et al., , 2016b(Coon et al., , 2020;;Valzania et al., 2018). Conversely, W. smithii larvae can also affect the composition and trophic structure of inquiline microbial communities (Peterson et al., 2008), with implications for community function (Nemergut et al., 2013). Changes to microbial metabolism and degradative activity may directly impact nutrient supply to host plants, as they are largely dependent on microbes for enzymatic lysis of prey substrates (Young et al., 2018). The effect of W. smithii on culturable bacterial diversity and abundance has been previously examined, suggesting that larval grazing on microeukaryotes such as rotifers and protozoa may reduce predation pressure on bacteria, increasing bacterial density and diversity (Cochran-Stafira & von Ende, 1998;Hoekman, 2011;Kneitel & Miller, 2002;Peterson et al., 2008). Even so, other studies have shown that direct bacterivory by larvae may also play a role in structuring bacterial communities (Hoekman, 2007;Hoekman et al., 2009). While models or rationale have been proposed for the role of mosquito larvae in pitcher plant inquiline communities (Buckley et al., 2003;Mouquet et al., 2008), the impact of W. smithii on microbial community function in the S. purpurea system has not been directly measured in field populations.

In this study, we manipulated W. smithii larval density and characterized effects on bacterial community composition and function across a natural population of pitchers over time, from pitcher opening to maturation. Bacterial community composition was characterized by high throughput 16S rRNA gene amplicon sequencing, while microbial community function was characterized by profiling carbon substrate use, bulk metabolic rate, hydrolytic enzyme activity and macronutrient pools. We hypothesized that pitcher-associated bacterial communities would change over time, and that the presence and abundance of mosquito larvae would affect community assembly and function in the pitcher microecosystem.

Experimental manipulation of S. purpurea communities was conducted over a 74-day period from June to August 2020 at the Cedarburg Bog in Saukville, WI USA (43°23.2′ N, 88°0.63′ W) (Bott et al., 2008;Grothjan & Young, 2019) with approval from the University of Wisconsin-Milwaukee Field Station advisory committee (Figure 1a).

In brief, 60 unopened pitchers of similar size were tagged and allowed to open naturally. Plants were selected throughout accessible areas to capture population variability. Pitchers were mostly located on different plants, but several pitchers shared a common plant and subsequent analyses accounted for potential non-independence by including pitcher and plant ID as random effects (see "Other statistical analyses"). Newly opened pitchers were then flushed and filled with sterile water prior to being provisioned 10 freeze-dried adult Drosophila melanogaster flies cultured from commercial stocks (Carolina Biologicals) to standardize initial prey availability (= Day 0).

Seven days later, 20 pitchers were assigned to each of three treatments: "Enriched", "Depleted" or "Unmanipulated" (= Day 7). Pitchers were then visited every 10-14 days for a total of 74 days to (i) record the number of mosquito larvae, (ii) maintain a count of 10 larvae in Young, unpublished data), this represented early and late succession, although successional trajectories were not specifically examined.

Enzymatic activities were measured from samples collected on Days 7,18,28,39,and 74. All other nutrient and functional assays were conducted on pitcher fluid sampled on Day 74. All pitcher samples were transported on ice to the laboratory for immediate processing, while any removed volume was replenished with sterile water at each sampling timepoint to maintain pitchers at ~3/4 maximum volume.

Data confirmed that field manipulations were successful in maintaining "high", "medium" and "low" larval abundances in pitchers assigned to the Enriched, Unmanipulated and Depleted treatment groups, respectively, with Enriched pitchers containing significantly higher larval abundances than Depleted pitchers over the entire 74-day experiment (Dunn's test; BH-adj p < .001 for all timepoints) (Figure 1b).

At each DNA sampling timepoint, 1-5 mL of pitcher fluid was centrifuged at 20,000 rcf (× g) for 20 min and supernatant removed prior to storage of pellets at -20°C. Total DNA was extracted from preserved pellets using standard phenol-chloroform methods (Sambrook et al., 1989;Stevenson & Weimer, 2007) prior to one-step PCR amplification of the V4 region of the bacterial 16S rRNA gene using barcoded primers and paired-end sequencing (2 × 250-bp) on an Illumina MiSeq as previously described (Arellano & Coon, 2022). Initial sequence quality control was conducted in QIIME2 (Bolyen et al., 2019), including read-joining using VSEARCH (Rognes et al., 2016). Reads were denoised using qscore based filtering (Bokulich et al., 2013) prior to running Deblur to trim paired-end joined reads (Amir et al., 2017). Taxonomy was assigned using a Naïve-Bayes classifier trained against the Greengenes database (Kozich et al., 2013;McDonald et al., 2012;Pedregosa et al., 2012). Multiple sequence alignment was followed by construction of an unrooted phylogenetic tree using MAFFT and FastTree2 (Katoh & Standley, 2013;Price et al., 2010). All QIIME2 artefacts were converted to phyloseq readable files for data analysis in R (version 4.1.1) (McMurdie & Holmes, 2013). A final phyloseq object was generated following a decontamination procedure using the R package "decontam" (Davis et al., 2018).

Contaminant reads, samples with <100 reads, and any reads classified as "chloroplast" or "mitochondria" were then removed prior to conducting downstream analyses. Species richness and Shannon diversity were estimated using the R packages "breakaway" and "DivNet", respectively (Willis & F I G U R E 1 Study design for field-based manipulations at Cedarburg Bog, Wisconsin, USA (a), and the effects on larval abundance in pitchers assigned to different treatment groups over time (b). Sixty unopened Sarracenia purpurea pitchers were tagged in the field and visited shortly after opening to standardize prey capture by the addition of ten Drosophila melanogaster in 10-15 mL sterile water (= Day 0). Seven days later (= Day 7), 20 pitchers were assigned to each of three treatments: "Enriched" (yellow), "Unmanipulated" (green), and "Depleted" (blue). Water was also collected for bacterial 16S rRNA gene amplicon sequencing and enzymatic assays of protease and chitinase activity (= early timepoint, Day 7). At 10-14-day intervals thereafter, pitchers were visited to maintain treatment manipulations (see Days 18,28,39,53 and 63) and to collect samples for additional assays for enzyme activity (see Days 18,28 and 39). At the late timepoint (= Day 74), all remaining pitcher fluid was collected for bacterial 16S rRNA gene amplicon sequencing and profiling of carbon substrate use, bulk microbial respiration, and macronutrient concentrations. Larvae added to Enriched pitchers were derived from newly hatched eggs laid by adult females in our laboratory Wyeomyia smithii colony (established from mosquitoes collected from natural populations in Wisconsin) and were maintained under sterile conditions prior to field introductions. Boxplots show high, low and median values, with lower and upper edges of each box denoting first and third quartiles, respectively. Dashed lines are "loess" fit lines showing the trend of larval abundance in each treatment group over time. Asterisks (*) indicate timepoints where larval abundance was significantly higher in Enriched pitchers relative to Depleted pitchers (Dunn's test; BH-adj p < .001).

(a) (b) Bunge, 2015;Willis & Martin, 2020), while hypothesis testing for both metrics was conducted using the betta() function in the "breakaway" package (Willis et al., 2017). Beta diversity ordinations were performed on an Aitchison distance (dissimilarity) matrix using principal component analysis (PCoA) after centered log-ratio (clr) transformation of the feature abundance table, as implemented in the R package "microbiome" (Lahti & Shetty, n.d.;McMurdie & Holmes, 2013).

Global dissimilarity was assessed using a permutational multivariate analysis of variance (PERMANOVA) implemented via the adonis() function in the R package "vegan" (Oksanen et al., 2020), with pairwise tests corrected using the Bonferroni method (Martinez Arbizu, 2017). Global and Bonferroni-adjusted pairwise differences in variance were assessed using the betadisper() function in "vegan" (Oksanen et al., 2020). Differentially abundant taxa across different categorical covariate groups were identified via Wilcoxon rank sum tests (two groups) or Kruskal-Wallis tests (three or more groups) using the R package "ALDEx2" (Fernandes et al., 2014), while relationships between different taxa and continuous covariates were identified via ALDEx2-estimated Pearson correlations (Fernandes et al., 2014).

All ALDEx2-generated p-values were corrected using a Benjamini-Hochberg false discovery rate (FDR) adjustment.

At each sampling timepoint for enzymatic activities, 500 μL of pitcher fluid, filtered through 15 μm mesh, was used to assay chitinase and protease activity within 5 h of collection using fluorometric microplate assays to accommodate small volumes as previously described (Young et al., 2018). In brief, chitinase activity was assayed with 4-methylumbelliferyl N-acetylβd-glucosaminide as a substrate (EX 360 nm, EM 420 nm), and protease activity was measured as leucine-aminopeptidase activity using l-leucine-7-amido-4-methylcoumarin hydrochloride as a substrate (EX 460 nm, EM 380 nm) (all Sigma-Aldrich).

Macronutrient concentrations were assayed using a Hach DR900 kit and reagents (Hach Company, Loveland, Colorado), according to the manufacturer's protocols including blanks and standards. Total organic carbon (TOC) and total phosphorus (TP) concentrations were determined using samples diluted 1:5 with organic-free water (Hach Company, Loveland, Colorado) to 10 mL and 25 mL total volumes, respectively, while total nitrogen (TN) was determined using undiluted samples (2 mL pitcher fluid).

A subset of pitcher fluid samples (10 from each treatment) was assayed at the end of the experiment for the capacity of associated bacterial communities to hydrolyze 31 carbon substrates using the Biolog Ecoplate system (Choi & Dobbs, 1999) (Table S1). In brief, each sample was diluted 1:10 with sterile water (pH adjusted to 5.6) and 100 μL of diluted sample was added to each Ecoplate well as previously described (Bittleston et al., 2020). Plates were incubated at 25°C for 3 days and the hydrolysis efficiency for each C substrate was assessed by measuring absorbance at 590 nm (Synergy™ HT; Biotek Instruments, Winooski, Vermont).

The net respiration rate for each pitcher-associated microbial community was also determined at the end of the experiment using a MicroResp™ device (Bittleston et al., 2020;Campbell et al., 2003). In brief, 250 μL of undiluted pitcher fluid was incubated with 100 μL of sterile cricket media slurry prepared as a suspension of food-grade cricket powder (JR Unique Foods) in water at a concentration of 3 g/L (pH adjusted to 5.6) (Bittleston et al., 2020). The deep-well plate was topped with a 96-well plate containing cresol red agar to indicate pH reduction. After 3 days in the dark at 25°C, absorbance of the top plate was measured at 570 nm using a microplate reader. Absorbance values were interpolated to rates of carbon dioxide (CO 2 ) production using a standard curve generated from incubation of known CO 2 concentrations with indicator medium.

In order to examine community composition independently from our a priori timepoints and treatments, samples were clustered based on community composition using a partitioning around medoid (PAM) algorithm implemented using the R package "biotypeR" (Tap & Mende, 2016). This method identifies tentative clusters ("biotypes") using a dissimilarity matrix and assigns centrally located points with minimal within-cluster dissimilarity as "medoids". Cluster computation is then employed to iteratively minimize the dissimilarity between medoid and non-medoid points. To generate the dissimilarity matrix for PAM clustering, sequence data were collapsed to the genus level prior to clr transformation and estimation of Aitchison distances as described above. Validation of the optimal number of clusters and statistical significance of clustering maximized the computed Calinski-Harabasz (CH) indexes (Calinski & Harabasz, 1974) for 1-20 possible clusters and calculated Silhouette coefficients (Rousseeuw, 1987).

Differentially abundant genera among the two and three biotype clusters generated from samples collected at the early ( Day 7) and late (Day 74) timepoints were identified using two independent approaches. First, the "ade4" package in R was used to perform between-class analysis (BCA) of PAM-identified clusters and compute BCA-indices for every taxon in each cluster. High BCA indices for taxa within clusters indicated enrichment relative to other clusters, while low and/or negative values indicated relative paucity. The second approach utilized a Kruskal-Wallis test carried out on the clr-transformed phyloseq object using ALDEx2 in R and identified differentially abundant taxa at the global level (Fernandes et al., 2014). A consensus approach was then used to conservatively identify taxa driving PAM clustering and thus distinguishing biotypes in early and late samples.

Analyses of carbon substrate hydrolysis and macronutrient concentrations by a categorical group, treatment, or biotype, were conducted using one-way ANOVA followed by post-hoc Tukey-Kramer Relationships between different functional assay outputs and larval abundance in individual pitchers were examined using either simple linear regression or linear mixed effects models estimated via the R package "lme4", with pitcher and plant included as nested random effects. Larval abundance was defined in one of two ways, depending on the response variable: (i) as the total number of larvae in a given pitcher at the time of sampling (response variable = protease or chitinase activity), or (ii) the total number of larvae counted in a given pitcher over the entire sampling period (all other response variables; hereafter referred to as "cumulative larval abundance").

Significance of fixed effects in models was calculated using Type III ANOVA with Satterthwaite's method for approximation of degrees of freedom (implemented via the "lmerTest" package) (Kuznetsova et al., 2017), while simulation-based restricted likelihood ratio tests (implemented via the "RLRsim" package) (Scheipl et al., 2008) were used to determine if the variance of random effects significantly differed from zero. The only exception to the above pertained to the relationship between measures of microbial respiration and cumulative larval abundance, which was non-linear and therefore better approximated by a second order polynomial function.

We identified a total of 1103 bacterial ASVs from 16 phyla across all samples (Table S2), with four phyla comprising 87% of total ASVs: Proteobacteria (54%), Bacteroidetes (19%), Actinobacteria (8%) and Firmicutes (5%) (Figure 2a). Alpha diversity, measured by estimates of both species' richness and evenness (Shannon's H index), significantly increased between early ( Day 7) and late (Day 74) samples, independent of any treatment effects (betta test, p < .001) (Figure 2b; Table S2). Among early samples, Enriched pitchers (with high larval abundances were maintained throughout the experiment; Figure 1b) exhibited significantly lower species richness values-but not Shannon's H index values-relative to Depleted pitchers (in which larvae were absent or present at very low abundances) (betta test, p < .001) (Figure 2b; Table S2).

A similar, but more dramatic pattern was observed among late samples, with Enriched pitchers exhibiting significantly lower species richness values (betta test, p = .006) and both Enriched and Unmanipulated pitchers (in which larval abundances naturally varied from 0 to 11 throughout the experiment) exhibiting significantly lower Shannon's H index values than Depleted pitchers (betta test, p < .01) (Figure 2b; Table S2).

Bacterial communities also significantly differed in composition between early and late samples, with 102 ASVs identified as differentially abundant between Day 7 and Day 74 (Table S3). Early samples contained notably greater proportions of taxa within the Proteobacteria (BH-adj p < .001) and Bacteroidetes (BH-adj p < .001), while late samples contained greater proportions of taxa within the Planctomycetes (BH-adj p < .001), Acidobacteria (BH-adj p < .001) and Verrucomicrobia (BH-adj p = .02) (Figure S1). These observations were further supported by ordination analysis using a beta diversity index (Aitchison distance), which revealed significant clustering of samples by timepoint (i.e. early vs. late; Figure S2A). In contrast, no taxa were identified as differentially abundant between Enriched, Unmanipulated, and Depleted pitcher samples at either timepoint, although ordination analysis using the same beta diversity index revealed a marginally significant treatment group effect in late samples driven by differences between Enriched and Depleted pitchers (Figure S2B,C; Table S4).

To address whether specific taxa are associated with the presence and abundance of mosquito larvae in pitchers independent of the pitchers' assigned treatment group, we also conducted differential abundance testing using models that assessed whether individual taxa in pitchers sampled during early and late timepoints covaried with the total number of larvae in each pitcher at the time of sampling. We identified two ASVs that were negatively associated with larval abundance, corresponding to the genera Gordonia (Order: Actinomycetales; BH-adj p = .001) and Sphingomonas (Order: Sphingomonadales; BHadj p = .02). Nine bacterial ASVs were positively correlated with larval abundance, four of which could be identified to genus: Cryocola (Order: Actinomycetales; BH-adj p = .008), Rhodanobacter (Order: Xanthomonadales; BH-adj p = .002), and Novosphingobium (Order: Sphingomonadales; BH-adj p = .001). An additional five positively correlated ASVs were identified to families Acetobacteriaceae (Order: Rhodospirillales; BH-adj p = .04), Methylocystaceae (Order:

Rhizobiales; BH-adj p = .04), Isosphaeraceae (Order: Gemmatales; BHadj p = .03), Sinobacteraceae (Order: Xanthamonadales; BH-adj p = .02) and Comamonadaceae (Order: Burkholderiales; BH-adj p = .001). One other positively correlated ASV had a taxonomic assignment resolved to the order level: Burkholderiales (BH-adj p = .04).

Studies in other aquatic systems support the capacity of environmental alterations to induce changes in microbial community performance that may or may not be associated with changes in community composition (Bier et al., 2015;Comte et al., 2013). In this way, W. smithii larvae have the capacity to shape ecosystem function in the S. purpurea microbial community via impacts on individual taxa and/or community-level processes. We therefore sought to elucidate larval impacts on key measures of pitcher microbial community function, including measures of community hydrolytic enzyme activity, carbon substrate use capacity, and respiration, as well as macronutrient concentrations in pitcher fluid. The activity of two hydrolytic enzymes involved in prey digestion (protease and chitinase) was characterized in fluid sampled from individual pitchers at five timepoints throughout the experiment (Figure 1a), while the capacity of pitcher-associated communities to hydrolyze each of 31 carbon substrates, community bulk respiration rates, and TN, TOC and TP concentrations in pitcher fluid were only characterized for the late timepoint at the end of the experiment (Figure 1a). Linear mixed effects models identified time as well as both discrete and continuous measures of larval abundance (i.e. treatment group and larval abundance counts at the time of sampling) as significant fixed effects in predicting protease activity (Table 1; Table S2),

with Enriched pitchers and pitchers with higher larval abundances

Relative abundance of bacterial orders in fluid collected from "Enriched", "Unmanipulated", and "Depleted" pitchers at the early (upper, Day 7) Alphaproteobacteria Caulobacterales Rhizobiales Rhodospirillales Rickettsiales Sphingomonadales Acidobacteria Acidobacteriales Actinobacteria Acidomicrobiales Actinomycetales uncl. Alphaproteobacteria Betaproteobacteria Burkholderiales Neisseriales SC-I-84 Deltaproteobacteria Bdellobivrionales Myxococcales Gammaproteobacteria Enterobacteriales Legionellales Pseudomonadales Xanthomonadales Chloroflexi JG30-KF-CM45 Bacteroidetes Bacteroidales Cytophagales Flavobacteriales Saprospirales Sphingobacterales Verrucomicrobia Chthoniobacterales Methylacidiphilales Pedosphaerales Firmicutes Clostridiales Planctomycetes Gemmatales

Other Proteobacteria Proteobacteria (cont.) 1365294x, 2024, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17314 by University of Wisconsin Milwaukee, Wiley Online Library on [25/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License exhibiting significantly lower measures of protease activity than Depleted pitchers or pitchers with lower larval abundances (Table 1; Table S2). In contrast, only time was identified as a significant fixed effect in predicting chitinase activity when larval abundance was coded as a discrete variable (Table 1; Table S2), while both time and larval abundance were identified as significant fixed effects in predicting chitinase activity when larval abundance was coded as a continuous variable (Table 1; Table S2).

We also observed significant interactions between time and larval abundance in predicting both protease and chitinase activity, with temporal effects on enzymatic activity being differentiated by larval abundance early in our sampling (Table 1; Figure S3; Table S2).

This was further consistent with the observation that measures of both protease and chitinase activity generally peaked early in our sampling and declined over time (Figure S3; Table S2).

Communities in pitchers assigned to different treatment groups also significantly differed in their capacity to hydrolyze certain carbon substrates, namely glycogen (ANOVA on log-transformed dependent variable, F 2,24 = 3.67, p = .04), with measures of glycogen hydrolysis efficiency being significantly lower in Enriched pitchers than in Depleted pitchers (Tukey-Kramer HSD test, Bonferroniadj p = .03) (Figure 3a; Table S2). Fluid collected from pitchers assigned to different treatment groups further contained different concentrations of the macronutrient nitrogen (ANOVA on square root-transformed dependent variable, F 2,24 = 3.82, p = .04), with TN concentrations being significantly higher in Enriched pitchers than in Depleted pitchers (Tukey-Kramer HSD test, Bonferroni-adj p = .03) (Figure 3a; Table S2). Consistent with both of these treatment group effects, linear regression models revealed a significant negative association between cumulative larval abundance and the efficiency of the associated community's capacity to hydrolyze glycogen (linear regression on log-transformed dependent variable, R 2 = .39, p < .001) (Figure 3b; Table S2). These models also revealed significant positive associations between cumulative larval abundance and both TN concentration (linear regression on square root-transformed dependent variable, R 2 = .36, p < .001) and l-phenylalanine hydrolysis efficiency (linear regression, R 2 = .13, p = .04) (Figure 3b; Table S2), although l-phenylalanine hydrolysis did not significantly differ across pitcher treatment groups.

Enriched and Depleted pitchers (Figure 3a; Table S2). In contrast, measures of community respiration as a proxy for microbial metabolic activity exhibited a more quadratic relationship with larval abundance, with bulk respiration rates being significantly higher in Enriched pitchers than in Unmanipulated pitchers (Dunn's test, BH-adj p = .002) but statistically similar as rates in Depleted pitchers (Figure 3a,b; Table S2). There were no significant differences in TOC, TP or hydrolysis of the other profiled carbon substrates between treatment groups or as a function of larval abundance (Table S2). Furthermore, no particular taxonomic group or ASV was significantly correlated with measures of enzymatic activity, bulk respiration, or carbon substrate hydrolysis (p > .05). TA B L E 1 Model parameters for mixed effects models of enzymatic activity.

Fixed effect Fixed effect coefficient estimate Test statistic and p-value Random effect Test statistic and p-value log(Protease) Treatment Unmanipulated-Depleted: -0.018 Unmanipulated-Enriched: -0.196 t: -0.259; p = .796 t: -2.688; p = .008** Pitcher ID RLRT < 0.001; p = .481 Sampling date -0.006 F 1,243 : 30.405; p < .001*** Plant ID RLRT = 0.922; p = .148 log(Chitinase) Treatment Unmanipulated-Depleted: 0.059 Unmanipulated-Enriched: -0.079 t: 0.758; p = .450 t: -0.992; p = .323 Pitcher ID RLRT < 0.001; p = .481 Sampling date -0.009 F 1,239 : 53.904; p < .001*** Plant ID RLRT = 3.677; p = .024* log(Protease) Larval abundance -0.068 F 1,247 : 18.719; p < .001*** Pitcher ID RLRT <0.001; p = 1 Sampling date Interaction -0.008 0.002 F 1,232 : 27.912; p < .001*** F 1,245 : 7.938; p = .005** Plant ID RLRT <0.001; p = 1 log(Chitinase) Larval abundance -0.049 F 1,254 : 9.326; p = .002** Pitcher ID RLRT <0.001; p = 1 Sampling date Interaction -0.011 0.001 F 1,233 : 50.692; p < .001*** F 1,243 : 8.851; p = .003** Plant ID RLRT <0.001; p = 1 ***0 < p < .001; **.001 < p < .01; *.01 < p < .05.

Owing to known limitations associated with traditional approaches to microbiome data analysis, which can result in overgeneralization (e.g. by comparing observations by alpha diversity measures) and/ or overspecification (e.g. by performing differential abundance tests for every taxon, which incurs a heavy penalty for multiple testing) of the data, our last series of analyses employed an unsupervised clustering method to identify naturally occurring clusters of samples.

To this end, we used the PAM clustering method-a distance-based method originating from studies of host-associated microbiomes (Arumugam et al., 2011) and that has been proposed as a robust approach to cluster microbial samples from a range of environments (Cheng & Ning, 2019).

PAM clustering of communities across the entire dataset identified two distinct biotypes corresponding to samples from the early and late timepoints, with the exception of a single early (Day 7) sample that clustered with late (Day 74) samples (Figure 4a). Reclustering of communities separated by timepoint subsequently identified two and three distinct biotypes among early and late samples, respectively (hereafter referred to as "Early Biotype 1", "Early Biotype 2", "Late Biotype 1", "Late Biotype 2" and "Late Biotype 3"; Figure 4b,c). Early biotypes were distinguished by five taxa from the orders Actinomycetales, Micrococcales and Xanthamonadales that were enriched in Early Biotype 1 relative to Early Biotype 2 (Table S5). Estimates of both species' richness and evenness (Shannon's H index) were also significantly higher for communities assigned to Early Biotype 1 than for communities assigned to Early Biotype 2 (betta, p < .05; Table S2). However, communities assigned to either early biotype were equally likely to have been derived from Enriched, Unmanipulated or Depleted pitchers (χ 2 = 0.37, p = .831) with statistically similar measures of cumulative larval abundance and hydrolytic enzyme activity (Mann-Whitney U test, p > .05). In contrast, the vast majority (~83%) of late communities assigned to Late Biotype 1 were derived from Enriched or Unmanipulated pitchers, while most (~64%) communities assigned to Late Biotype 2 were derived from Depleted pitchers (Figure 4c). Accordingly, measures of cumulative larval abundance were significantly higher in pitchers harbouring communities assigned to Late Biotype 1 than in pitchers harbouring communities assigned to Late Biotype 2 (Dunn's test, p = .04). Communities assigned to Late Biotype 1 were also more likely to be derived from pitchers assigned to Early Biotype 2, although this association was only marginally significant (Fisher's exact test, p = .049).

Consistent with the observation that late timepoint samples from Enriched pitchers harboured communities of lower alpha diversity than samples from Unmanipulated or Depleted pitchers, estimates of alpha diversity were significantly lower for communities assigned to Late Biotype 1 than for communities assigned to Late Biotype 2 or Late Biotype 3 (betta test, p < .001) (Table S2). Fortyfour taxa were also identified as significant drivers of late biotypes, including genera like Novosphingobium and Sphingomonas that were enriched in Late Biotype 1 and Late Biotype 2 and were previously shown to be positively and negatively correlated with larval abundance, respectively (Figure 5; Table S5). However, ASV diversity remained significantly lower for communities assigned to Late Biotype 1 than for communities assigned to Late Biotype 2, regardless of whether ASVs assigned to any biotype-enriched genus or only genera specifically enriched in communities assigned to Late Biotype 1 were considered (Tukey-Kramer HSD test, Bonferroni-adj p < .05;

Figure S4; Table S2).

Finally, fluid collected from pitchers harbouring communities assigned to Late Biotype 1 contained higher TN concentrations than pitchers harbouring communities assigned to Late Biotype 2 (Tukey-Kramer HSD test, Bonferroni-adj p = .02; Table S2), consistent with TN concentrations being significantly higher in fluid collected from pitchers with higher larval abundance than in fluid collected from Depleted pitchers (Figure 3a; Table S2). Communities assigned to Late Biotype 2 also exhibited significantly higher capacities to use Tween40 as a carbon substrate than communities assigned to Biotype 1 (Tukey-Kramer HSD test, Bonferroni-adj p = .003) (Table S2), consistent with hydrolysis efficiencies for other complex carbohydrate polymers (e.g. glycogen) being (i) significantly lower for communities from Enriched pitchers than for communities from Depleted pitchers, and (ii) negatively correlated with cumulative larval abundance in our previous analyses (Figure 3a; Table S2).

In contrast, communities assigned to Late Biotype 2 exhibited protease and chitinase activities that did not statistically differ from levels of communities assigned to Biotype 1 (ANOVA, F 1,26 = 1.863 and 3.681; p = .184 and .070 for protease and chitinase, respectively;

Table S2), despite measures of protease activity being significantly higher in Depleted pitchers than in Enriched pitchers in our treatment group comparisons (Table 1; Table S2).

Here, we used the model aquatic microecosystem hosted by the car- (Day 74) communities assigned to "Late Biotype 1", "Late Biotype 2", and "Late Biotype 3" (left), along with the relative abundances of select genera that were enriched in each biotype (right). Coloured bars in the left panel present the proportion of sequencing reads assigned to a given order, with low abundance orders (<2%) being represented by the "Other" category. Boxplots in the right panel show high, low, and median values, with lower and upper edges of each box denoting first and third quartiles, respectively. Only boxplots for the two genera with the highest BCA coefficients in each biotype are presented. Other taxonomic drivers of biotypes can be found in Table S5. 0.00 0.05 0.10 0.15 0.20 Rel. abundance Burkholderia 0.0 0.2 0.4 0.6 Rel. abundance Novosphingobium 0.00 0.01 0.02 0.03 0.04 Rel. abundance Sphingomonas 0.00 0.01 0.02 0.03 1 2 3 Biotype Rel. abundance Chthoniobacter 0.000 0.005 0.010 0.015 1 2 3 Biotype Rel. abundance Phenylobacterium 0.00 0.01 0.02 0.03 Rel. abundance Edaphobacter Late Biotype 1 0.00 0.25 0.50 0.75 1.00 Proportion of reads 0.00 0.25 0.50 0.75 1.00 Proportion of reads 0.00 0.25 0.50 0.75 1.00 Proportion of reads Late Biotype 2 Late Biotype 3 Chloroflexi JG30-KF-CM45 Bacteroidetes Bacteroidales Cytophagales Flavobacteriales Saprospirales Sphingobacterales Verrucomicrobia Chthoniobacterales Methylacidiphilales Pedosphaerales Firmicutes Clostridiales Planctomycetes Gemmatales Cyanobacteria Other Alphaproteobacteria Caulobacterales Rhizobiales Rhodospirillales Rickettsiales Sphingomonadales Acidobacteria Acidobacteriales Actinobacteria Acidomicrobiales Actinomycetales uncl. Alphaproteobacteria Proteobacteria Betaproteobacteria Burkholderiales Neisseriales SC-I-84 Deltaproteobacteria Bdellobivrionales Myxococcales Gammaproteobacteria Enterobacteriales Legionellales Pseudomonadales Xanthomonadales Proteobacteria (cont.) 1365294x, 2024, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17314 by University of Wisconsin Milwaukee, Wiley Online Library on [25/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Bacterial taxa identified across the field data set were consistent with previous characterization of inquiline bacterial communities in S. purpurea (Bittleston et al., 2018(Bittleston et al., , 2020;;Freedman et al., 2021;Grothjan & Young, 2019, 2022;Northrop et al., 2017;Paisie et al., 2014), sharing broad patterns such as the presence of diverse members of the Alphaproteobacteria (Rhodospirillales, Rhizobiales, Caulobacterales, Sphingomonadales), Betaproteobacteria (Burkholderiales, Neisseriales) and Gammaproteobacteria (Enterobacteriales, Pseudomonadales, Xanthomonadales), as well as members of the Bacteroidetes (Saprospirales, Sphingobacteriales, Flavobacteriales), Actinobacteria (Soilrubribacterales, Actinomycetales) and Firmicutes (Clostridiales). Bacterial community diversity (both species' richness and evenness) also increased from early (7 days after pitcher opening) to late succession (74 days), with significant shifts in bacterial composition observable up to the phylum level. These results are consistent with observations for inquiline communities in other pitcher plant species (Armitage, 2017), but differ from previous studies in S. purpurea reporting either variable shifts in bacterial diversity (Grothjan & Young, 2022) or stable diversity through ecological succession (Miller & TerHorst, 2012). The latter study (Miller & TerHorst, 2012) used culture-based methods to assess bacterial diversity, as opposed to the culture-independent, sequencing-based methods used here, while the former study (Grothjan & Young, 2022) assessed bacterial diversity in pitchers on plants maintained under greenhouse conditions. Greenhouse experiments may incorporate effects of isolation from wetland microbial sources, reiterating the importance of dispersal processes in shaping pitcher community assembly (Boynton et al., 2019;Kneitel & Miller, 2003) and the relative loss of diversity in S. purpurea pitchers outside of a field context (Arellano & Coon, 2022). Diversity of late timepoint samples differed by treatment group, with communities in pitchers with higher mosquito larval abundances having fewer bacterial taxa with more uneven community compositions relative to those in pitchers with few or no larvae. That mosquito larval abundance influenced bacterial diversity from early to late timepoints suggests a role for top predators in mediating successional trajectories via top-down controls on bacterial communities, expanding on former reports of local correlations between W. smithii larval colonization and bacterial diversity (Peterson et al., 2008). Previous studies have suggested that the presence and abundance of other macroinvertebrates-but not mosquito larvae-is associated with higher alpha diversity in pitchers over a continental scale (Freedman et al., 2021), but our observation of lower species richness in pitchers with higher larval abundance within a single wetland site also suggests that the nature of W. smithii-microbiota interactions may vary across the S. purpurea range. Indeed, our results are more consistent with the hypothesis that increasing bacterial diversity within pitcher plant inquiline communities with increasing latitude (Buckley et al., 2003) may be driven by release of filter-feeding pressure by W. smithii, as pitcher mosquito larval density also decreases with increasing latitude (Bradshaw & Holzapfel, 1989). Interestingly, late timepoint communities in Enriched pitchers had significantly lower richness than communities in Unmanipulated and Depleted pitchers, but similar evenness (measured as Shannon's H index) relative to Unmanipulated pitchers.

Elevated richness in Depleted pitchers may suggest a release from predation pressure under field conditions.

While higher larval abundances were associated with significantly lower estimates of alpha diversity for both early and late timepoint communities, only mature communities at the late timepoint were distinguished by treatment group a priori-varying in composition and sample heterogeneity. This may suggest that mosquito-mediated shifts in bacterial composition intensified over time when larval abundance treatment differences were maintained with regular manipulations. Eleven bacterial ASVs exhibited relative abundances that scaled significantly with mosquito larval abundance. This could relate to previous empirical findings that gnotobiotic W. smithii larvae show different developmental outcomes with different bacterial taxa present (Arellano & Coon, 2022). Although specific relationships between each of these taxa and larval developmental outcomes are unclear, this is consistent with our understanding that W. smithii-microbe interactions are more specific than in other mosquito species (Coon et al., 2014(Coon et al., , 2016a(Coon et al., , 2016b(Coon et al., , 2020;;Valzania et al., 2018). This also supports our previous finding that pitcher bacterial communities are curated and actively dispersed by W. smithii mosquitoes (Arellano & Coon, 2022).

In addition to mosquito-mediated impacts on bacterial community composition, mosquito larval abundance was linked to changes in several community-level processes important for ecosystem function, and these associations were generally stronger than those observed between larval abundance and bacterial diversity. This suggests that, in addition to top-down effects of larval predation, the presence of larvae may influence bacterial gene expression related to community functions. Mature pitcher communities differed by treatment in several functional measures, including the capacities of communities to hydrolyze glycogen, bulk community respiration rates and TN concentrations in pitcher fluid. With the exception of community respiration, these relationships seemed to scale linearly with larval abundance; higher larval abundance was correlated with lower community-level glycogen hydrolysis rates and higher TN concentrations, suggesting mosquito larvae influence nutrient cycling and thus the chemical environment, resulting in responses of microbial metabolic capacities. Decline in glycogen hydrolysis rates with more larvae could relate to larval excretion and egestion providing alternative carbon and nitrogen-rich substrates that influence expression of different bacterial catabolic pathways, as occurs in nitrogen salvage by gut bacteria during host hibernation in mammals (Regan et al., 2022). Indeed, higher larval abundances also correlated with higher capacities for communities to hydrolyze l-phenylalanine.

While declines in the capacities of bacterial communities to utilize different substrates could relate to lower bacterial densities due to predator grazing, the positive relationship between larval abundance and L-phenylalanine hydrolysis efficiency also strongly suggests that mosquito-mediated decline in bacterial cell density was not the sole feature dictating metabolic features of late succession microbial communities.

As noted above, W. smithii mosquitoes have a unique association with S. purpurea pitchers, with larvae of most other mosquito species developing in a variety of natural and artificial water container habitats (Day, 2016). Adult females of most mosquito species also must usually feed on blood from a vertebrate host, which, together with teneral reserves acquired during larval feeding, provides essential nutrients for egg production and underlies the capacity of some species to transmit disease-causing organisms to humans (Takken & Verhulst, 2013;Timmermann & Briegel, 1999). However, some mosquito species, including W. smithii, have evolved to be autogenous, with adult W. smithii females in northern populations (including those in Wisconsin) producing eggs without ever blood feeding (Bradshaw, 1980;Bradshaw et al., 2022;Wolff & Riffell, 2018).

Previous studies in other autogenous species implicate enhanced nutrient acquisition as larvae, along with genetic factors, in facilitating the evolution of autogeny (Lounibos et al., 1982;Steffan & Evenhuis, 1981;Tsuji et al., 1990). Our previous work further supports roles for mosquito-associated microbes in shaping teneral reserves and mobilization of key nutrient stores (particularly protein and lipid stores) during oogenesis in adult females (Coon et al., 2016a(Coon et al., , 2020)). In the S. purpurea system, pitcher-associated microbial communities are known to have distinct genomic potential for amino acid metabolism (Bittleston et al., 2018), the expression of which is at least partially related to prey capture (Northrop et al., 2017). A positive correlation between mosquito larval abundance and hydrolysis of l-phenylalanine suggests that W. smithii larvae may affect amino acid turnover in pitcher communities, particularly of a metabolically expensive aromatic essential metabolite for which bacteria are frequently auxotrophic (D 'Souza et al., 2014;Zengler & Zaramela, 2018). In addition to the potential implications for the composition and underlying interactions of pitcherassociated communities, phenylalanine, along with tyrosine, has previously been identified as a key metabolite in the metabolism of non-blood-feeding relative to blood-feeding subspecies of the common house mosquito, Culex pipiens (Beintema et al., 1994;Siperstein et al., 2022;Wheeler & Buck, 1996). Whether phenylalanine hydrolysis is also important for the maintenance of a non-blood-feeding life history in W. smithii warrants further study, especially in light of results from a previous study that failed to detect differences in phenylalanine metabolism between blood-feeding and non-bloodfeeding W. smithii populations (Bradshaw et al., 2018). Whether the symbiotic association of W. smithii with S. purpurea and inquiline microbial communities may have resulted in an "out-sourced" route to a non-blood-feeding life history also warrants further study, given the broader implications for our understanding of the factors driving mosquito evolution away from life histories important for disease transmission.

Higher TN in the presence of W. smithii larvae is consistent with other studies reporting higher ammonia release in pitcher fluid related to presence of inquiline larvae (W. smithii and/or the pitcher plant midge, Metriocnemus knabi) (Bradshaw & Creelman, 1984).

However, prior research suggested the presence of a "complete" pitcher-associated food web (i.e. predatory pitcher plant flesh fly, Fletcherimyia spp., larvae; W. smithii larvae; M. knabi larvae; and other pitcher-associated eukaryotes) may not translate to greater nitrogen pools within host plant tissue (Bradshaw & Creelman, 1984;Butler et al., 2008). The influence of W. smithii larval density on nitrogen concentrations may also depend on interactions with intermediate trophic levels and vary as a function of prey (bacterivorous protozoa or rotifers versus direct bacterivory) (Mouquet et al., 2008). This study offers evidence that W. smithii can influence key macronutrient pools in the S. purpurea inquiline ecosystem, although some nitrogen may be used within the community. While explicit consideration of other food web components could help, our results are consistent with the idea that nitrogen release from captured insect prey is facilitated first by mechanical break-down by inquilines with a shredding feeding modality (i.e. Fletcherimyia spp. and M. knabi), followed by enzymatic hydrolysis by bacteria. Filter feeding W. smithii may then act to promote mineralization of nitrogen by preying upon microorganisms and releasing bioavailable nitrogen as a bottom-up effect in addition to the top-down effects of grazing (Bradshaw & Creelman, 1984;Heard, 1994).

Activity of the hydrolytic enzymes chitinase and protease changed over time, consistent with past experiments characterizing enzymatic activity in response to insect prey additions (Young et al., 2018). Despite strong temporal trends, a mixed effects modelling approach supported a trajectory of pitcher enzymatic activity as partially attributable to W. smithii larval abundance, whereby higher larval numbers correlated with lower chitinase and protease activity.

Notably, for chitinase activity, plant identity was important, potentially due to pitcher differences in natural prey capture over time after the initial standardized feeding (Grothjan, 2021). However, while we did not explicitly control for variability in prey capture rates over the season, our selected pitchers were within an area of similar prey availability and both plant and pitcher ID were included as random effects to at least partially account for potential prey capture biases. Protease production can be stimulated by relatively limiting N availability (Young et al., 2018); higher TN with more larvae may alleviate N limitation, reducing protease expression. Degradative activity in S. purpurea pitchers is largely attributable to microbial enzymes (Adlassnig et al., 2011;Young et al., 2018); thus, lower activity with higher larval abundance may also result from larval predation decreasing microbial biomass, although several functional assays scaled non-linearly or positively with larval abundance at the late timepoint. For example, significantly higher bulk microbial community respiration was observed in pitchers with higher mosquito larval abundance. The parabolic relationship between larval abundance and respiration suggests that microbial respiration is highest under

1365294x, 2024, 7, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17314 by University of Wisconsin Milwaukee, Wiley Online Library on [25/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 13 of 17 ARELLANO et al. low and high predation, but is lowest under intermediate predation,

or larval abundance, so community respiration rate may relate to an interaction of physiological status related to nutrient resources as well as cell density, related to predation. A similar relationship has been proposed for protozoan predation on bacteria under N-limited conditions (Hunt et al., 1977), analogous to the relatively low nutrient pitcher environment (Grothjan & Young, 2022). Furthermore, enhanced nutrient mineralization and stimulation of microbial respiration have been observed in response to higher predation pressure (Barsdate et al., 1974;Stout, 1973).

Overall, the more pronounced shifts in function with varying mosquito larval abundance despite more marginal shifts in bacterial community composition suggest differential expression of functions in response to larval abundance. However, while we were able to observe marked differences in bacterial composition between Days 7 and 74, more frequent sampling of community composition is necessary to demonstrate successional patterns and contextualize the higher temporal resolution of our enzymatic activity data.

Experiments are also warranted to better understand how changes in abiotic conditions interplay with mosquito larval abundance to influence pitcher-associated microbial communities in the field, although such changes will be difficult to resolve at the short time scales relevant for most microbes without disrupting community succession. We further recognize that, while our functional assays reflected the activity of complete pitcher-associated microbial communities, our DNA sequencing focused only on bacteria and non-bacterial components may contribute to metabolic activity (Young et al., 2018). Hydrolytic enzyme activity has been associated with bacteria (Young et al., 2018), but respiration and carbon substrate use could have contributions from fungal or other microeukaryotes, which could also be affected by top-down pressures of W. smithii (Boynton et al., 2019;Paisie et al., 2014). As noted above, our results could also reflect release from bacterivore predation by other microeukaryotes such as rotifers and protozoa that, like bacteria, are directly grazed upon by W. smithii (Cochran-Stafira & von Ende, 1998;Hoekman, 2011;Kneitel & Miller, 2002;Peterson et al., 2008). Indeed, future research should explore the varying degrees of W. smithii omnivory and its role in mediating potential trophic cascades in this system (sensu Paine, 1980). Varied intra-and inter-trophic processes may indirectly influence mosquito-bacterial interactions and thus pitcher-associated microbial community composition and function. et al., 2018). Another mosquito-associated taxon, Burkholderia, is also comprised of members that are well known to engage in symbiotic associations with diverse insect species (Kaltenpoth & Flórez, 2020). Indeed, field-derived isolates of the Burkholderia-Paraburkholderia complex have also been shown to support robust development of gnotobiotic W. smithii larvae in a laboratory setting (Arellano & Coon, 2022), suggesting they may also be functionally significant in this symbiosis. Functionally, the mosquito-depleted Late Biotype 2 also demonstrated low TN relative to the other biotypes, reinforcing the positive correlation between larval abundance and TN. In contrast, biotypes identified via unsupervised clustering failed to recapitulate associations between mosquito larval abundance and hydrolytic enzyme activity, despite previous analyses identifying a significant negative association between mosquito larval abundance and protease activity. These differences underscore how myriad factors in the field are likely to additionally shape late succession communities, as has been explored in previous studies investigating temporal trajectories of pitcher or pitcher-derived bacterial communities in field and greenhouse settings, respectively (Bittleston et al., 2021;Grothjan & Young, 2022). That early biotypes non-randomly associated with late biotypes in our study is also consistent with previous work in pitcher systems characterizing the importance of historical contingency in community progression (Bittleston et al., 2020). Future work is necessary to address how different variables in the field shape pitcher microbial communities independent of, or in concert with, top-down effects of W. smithii

to affect community-level processes underlying key functions in the S. purpurea model microecosystem.

This study supports a role for W. smithii mosquitoes in top-down regulation of the S. purpurea system through modulation of bacterial diversity, community composition and ecosystem function. Our novel findings link functional traits of the pitcher-associated microbial community to larval abundance over time, even in the presence of strong temporal structuring and field conditions. Functional