The microbiome and metatranscriptome of a panel from the Sarracenia mapping population reveal complex assembly and function involving host influence

Sarracenia provide an optimal system for deciphering the host-microbiome interactions at various levels. We analyzed the pitcher microbiomes and metatranscriptomes of the parental species, and F1 and F2 generations from the mapping population (Sarracenia purpurea X Sarracenia psittacina) utilizing high-throughput sequencing methods. This study aimed to examine the host influences on the microbiome structure and function and to identify the key microbiome traits. Our quality datasets included 8,892,553 full-length bacterial 16s rRNA gene sequences and 65,578 assembled metatranscripts with microbial protein annotations. The correlation network of the bacterial microbiome revealed the presence of 3-7 distinct community clusters, with 8 hub and 19 connector genera. The entire microbiome consisted of viruses, bacterial, archaea, and fungi. The richness and diversity of the microbiome varied among the parental species and offspring genotypes despite being under the same greenhouse environmental conditions. We have discovered certain microbial taxa that are genotype-enriched, including the community hub and connector genera. Nevertheless, there were no significant differences observed in the functional enrichment analysis of the metatranscriptomes across the different genotypes, suggesting a functional convergence of the microbiome. We found that the pitcher microcosm harbors both rhizosphere and phyllosphere microbiomes within its boundaries, resulting in a structurally diverse and functionally complex microbiome community. A total of 50,424 microbial metatranscripts were linked to plant growth-promoting microbial proteins. We show that this complex pitcher microbiome possesses various functions that contribute to plant growth promotion, such as biofertilization, bioremediation, phytohormone signaling, stress regulation, and immune response stimulation. Additionally, the pitcher microbiome exhibits traits related to microbe-microbe interactions, such as colonization of plant systems, biofilm formation, and microbial competitive exclusion. In summary, the demonstrated taxonomical divergence and functionally convergence of the pitcher microbiome are impacted by the host genetics, making it an excellent system for discovering novel beneficial microbiome traits.

Introduction

Plants live in association with a wide range of microorganisms, including bacteria, fungi, protists, nematodes, and viruses, all referred to as the plant microbiome. In natural habitats, through this association, microorganisms play a significant role in enhancing plant growth, productivity, and health (Trivedi et al., 2020). Experimental evidence has demonstrated that plants have the ability to modify the abundance and diversity of microbiomes in their environments (Mahnert et al., 2015). Nevertheless, the precise mechanisms through which this occurs remain unidentified. Plants may have acquired genetic traits that impact their interaction with their microbiomes. The identification of these traits and their mechanisms are essential for the success of the ongoing efforts to harness and manipulate microbiomes in order to optimize their advantages for the host plants (Ke et al., 2021; Afridi et al., 2022; Morales Moreira et al., 2023).

Carnivory evolved independently on 13 occasions in flowering plants, which includes four instances in monocots (three in Poales and one in Alismatales) and nine instances in eudicots (three in Caryophyllales, three in Lamiales, two in Ericales, and one in Oxalidales). There is a total of 810 species belonging to 21 plant genera, which are distributed across 13 families and six orders (reviewed in (Fu et al., 2023)). As a result of convergent evolution, some carnivorous species have developed modified leaves that allow them to capture and kill insect prey. These species also host a complex microbiome that aids in the digestion of the captured prey and support plant growth through an array of direct and indirect traits. These traits might be considered as extensions of the host plant traits, which are equally significant to the inherent plant traits (the concept is reviewed in (Whitham et al., 2006, 2008)). The “pitfall” carnivorous plants evolved to trap insects in three separate lineages (Albert et al., 1992). Using modified leaves called ‘pitchers’, they digest insect prey with the aid of a complex microbiome that develops inside their pitcher’s fluid. In the genus Sarracenia, species differ in their pitcher morphology, which in turn affects the insect species they capture for nutrients (Ellison et al., 2003a, 2004), and also affects the initial assembly of the pitcher microbiomes. The factors involved in shaping the microbiome community likely involve host plant genetics, local environmental conditions, as well as prey’s microbiomes. The interactions between the Sarracenia, its microbiome, and captured prey likely played a role in driving the process of plant speciation. In return, the host species genetics may have evolved traits to impact the composition and function of their microbiome communities. The host impact can be measured by examining the taxonomic divergence and functional convergence of the pitcher microbiome, which are two prevalent phenomena in microbial systems (Louca et al., 2016, 2018).

The bacteria, archaea and eukaryotes that reside in the Sarracenia pitchers have been described by organismal methods (Ellison et al., 2003b; Baiser et al., 2013) and by environmental sequencing (Adlassnig et al., 2011; Boyer and Carter, 2011; Koopman and Carstens, 2011). Multiple studies have explored the impact of the Sarracenia host plant on the microbiome assembly by comparing the microbial communities in natural and artificial pitchers (Koopman et al., 2010; Bittleston et al., 2018; Ellison et al., 2021; Grothjan and Young, 2022), as well as among various Sarracenia species (Heil et al., 2022). These studies concluded that Sarracenia may play an active role in shaping the structure of their microbiome. However, is this impact is regulated in part by the genetic makeup of the host (i.e., deterministic impact)? To address this question, we are using a Sarracenia mapping population in a controlled environment to characterize the microbiome diversity, structure and function among genotypes of same genetic background.

Malmberg et al. developed a Sarracenia genetic linkage map of 437 SNP and SSR markers with a total length of 2017 cM using an F2 generation of 280 plants from a genetic cross between Sarracenia purpurea (Spu) and Sarracenia psittacina (Spa) (Malmberg et al., 2018). They mapped 64 pitcher QTLs traits that were involved in prey attraction, capture and kill, and digestive processes/mechanisms. Despite the fact that the morphological characteristics of these traits are indicative of their adaptation to prey capture, their impact on the structure and function of the microbial community remains uncertain. We hypothesize that Sarracenia evolved genetic features for interacting with their vital microbiomes in a manner similar to the convergent evolution observed in their leaf structure.

Pitcher plants are ideal for studying host-microbiome and microbiome-microbiome interactions at the experimental and molecular levels for several reasons. Firstly, each fluid-containing pitcher is a natural microcosm with well-defined boundaries, making it experimentally amenable. This differs from other symbiotic microbial systems that are less amenable to experimentation, such as gut microbiomes, or those that do not have well-defined boundaries such as plant rhizosphere and phyllosphere microbiomes. Secondly, different species within the Sarracenia genus have varying leaf traits associated with differing microbial communities, that attract and digest different species of insects (Grothjan and Young, 2019). These genetically determined characteristics may have allowed pitcher plants to impose control over the microbiome assembly in their pitchers. Thirdly, the various species of Sarracenia can hybridize readily, offering a unique opportunity to incorporate a genetic approach to examine plant genome influences on microbiome assembly (Furches et al., 2013). Although pitcher plants have long attracted biologists, we know little about the assembly mechanisms involved in the microbiome assembly and structure in their pitchers (Grothjan and Young, 2019). It has been suggested that factors other than prey capture and colonization by eukaryotic species may affect the recruitment of bacteria to their pitchers (Grothjan and Young, 2019).

The objectives of this study are to evaluate the variations in the microbiome structure and function among different genotypes of the Sarracenia mapping populations, namely Spu, Sps, F1 and F2, and to characterize the composition and functions of the Sarracenia microbiome. For a panel of these genotypes, we utilized the single molecule (Oxford Nanopore Technology, ONT) and short reads (Illumina) technologies to sequence the full-length 16S rRNA gene (microbiome) and total mRNA (metatranscriptome), respectively, from their pitcher fluids (Figure 1A). We employed a variety of specialized databases, as well as bioinformatical tools to analyze the data (Figures 1B–D).

Figure 1

Figure 1. Summary of the experiment and analysis pipelines. (A) Photos of the parents (Spu and Sps) and selected F1 and F2 individuals of the Spu X Sps mapping population, and a sample of the collected pitcher fluid. (B) The experimental process of preparing the pitcher fluid samples and sequencing. (C) Modules of the 16S rRNA microbiome analysis pipeline. (D) Modules of the metatranscriptome RNA analysis pipeline.

Materials and Methods

Genotypes

The samples of parental species are true biological replicates because they are clones of the same individuals, but the samples from the F1 and F2 generations are separate individuals. To simplify the study, we treated the F1 individuals as biological replicates of each other and referred to them as “F1”. Similarly, the F2 individuals were treated as biological replicates and referred to them as “F2”. Comprehensive information regarding this mapping population is published in our earlier paper (Malmberg et al., 2018). All plants were grown under the same conditions in the greenhouse. Briefly: plants were grown in 4-6-inch pots on flats without holes and filled with 1/2 to 3/4 inch of water to simulate a tiny bog. The plants are top watered daily in the summer and less often in the autumn, winter, and early spring. The bottom inch of soil was immersed in standing water in the flats, keeping it saturated while the top 3 inches were moist but never saturated.

Collecting Sarracenia pitcher fluid

We collected the fluid from the pitchers of Spu, Sps, F1 and F2 plants (Figure 1A). All pitchers were at least 8 weeks old to ensure stable microbiome communities. The pitcher fluid was collected using sterile glass pipettes. For Spu, F1 and F2 plants, the glass pipettes were easily inserted into the pitchers as the flap is open. The Sps plants have a narrow body and a closed flap, which required a small incision in the flap to insert the glass pipette into the pitcher (Figure 1A). Due to the morphology of Sps, some pitchers were dry or contained a small amount of fluid. In these cases, sterilized water was added to elute the microbiome that is attached to pitcher’s internal surface. The fluid samples were collected into 50 ml conical sterilized tubes and stored in a -20 °C freezer.

Sample preparation

We thawed the pitcher fluid samples on ice for 30-40 minutes, followed by filtration using the Steriflip Sterile Centrifuge Tube Top Filter Unit (Cat# SCNY00020) to remove any plant material and dead insect matter. A schematic diagram of the sample preparation is shown in Figure 1B. The filtered fluid was a cloudy and homogenous mixture, indicating the presence of microbial growth. We centrifuged the samples at 3000 rpm at 4°C for 20 minutes to precipitate the microbiota, discarded the supernatant, suspended the pellet in 450 µl of PBS buffer, and mixed by pipetting. We added 750 µl of DNA/RNA Shield (Cat# R2002) to the suspension and mixed by pipetting and inverting. The material was transferred into a bead-bashing lysis tube and battered at high frequency for 2 minutes. Next, the samples were then centrifuged for 30 seconds at 16,000 rpm, divided into 400 µl aliquots in 2 ml tubes with 2 volumes of the DNA lysis buffer added to each tube. Finally, we used the ZymoBIOMICS DNA/RNA Miniprep Kit (Cat# R2002) to isolate both DNA and RNA.

Microbiome DNA isolation

We transferred the entire sample into a SpinAway Yellow Filter unit, followed by centrifugation at 16,000 rpm for 30 seconds. The flow-through, which contained the RNA, was saved in a 15 ml tube. An equal volume of 100% ethanol was added, and the tube was placed on ice. After processing the entire sample, a new collecting tube was used and 400 µl of DNA/RNA Prep Buffer was added to the column, followed by centrifugation at 16,000 rpm for 30 seconds, adding 700 µl of DNA/RNA wash buffer to the column, and centrifugation at 16,000 rpm for 30 seconds. Another 400 µl of the wash buffer was added to the column followed by centrifugation at 16,000 rpm for 2 minutes. We then carefully transferred the yellow filter into a nuclease-free 1.5 ml tube and added 50 µl of DNase/RNase-free water directly in the center of the yellow column matrix and incubated it at room temperature for 5 minutes. While waiting, we assembled a Zymo-Spin III-HRC Filter unit and added 600 µl of ZymoBIOMICS HRC Prep Solution into it, followed by centrifugation at 8000 rpm for 3 minutes. The elution column was then centrifuged at 16,000 rpm for 30 seconds, and the eluted DNA was transferred into the prepared HRC filter unit and centrifuged for 16,000 rpm for 3 minutes to elute the final DNA into1.5 ml tube.

Microbiome RNA extraction

We transferred the RNA sample into a SpinAway Green Filter unit, followed by centrifugation at 16,000 rpm for 30 seconds, replacing the collection tube, adding 400 µl of DNA/RNA buffer to the column, mixing by pipetting, and centrifugation at 16,000 rpm for 30 seconds. We added 700 µl of DNA/RNA wash buffer to the column and centrifuged at 16,000 rpm for 30 seconds, followed by adding another 400 µl of the wash buffer to the column, centrifugation at 16,000 rpm for 2 minutes. We then carefully transferred the green filter into a nuclease-free 1.5 ml tube, added 30 µl of DNase/RNase-free water directly in the center of the yellow column matrix, and incubated it at room temperature for 5 minutes. While waiting, we assembled the Zymo-Spin III-HRC unit, added 600 µl of HRC prep solution, and centrifuged at 8000 rpm for 3 minutes. The green column matrix was centrifuged at 16,000 rpm for 30 seconds to elute the RNA into the nuclease-free tube. The eluted RNA was transferred into the prepared HRC filter in another 1.5 ml nuclease-free tube and centrifuged at 16,000 rpm for 3 minutes to elute the final RNA in 1.5 ml tube. The RNA tubes were labeled and stored at -80˚C.

DNA and RNA quality assessment

We used the Fragment Analyzer Long Fragment (LF) assay to examine the DNA quality and the Bioanalyzer Pico RNA Assay to assess the RNA integrity. The concentrations were assessed using the Qubit 4 fluorometer high-sensitivity assays.

Full-length 16S rRNA gene libraries and sequencing

we used the Oxford Nanopore Technology (ONT) amplicon sequencing method to build the 16s rRNA gene libraries and sequence them on the MinIon Mk IC platform. First, we built the libraries using the 16S barcoding kit (SQK-16S024) following the manufacturer’s protocol. Briefly, the ONT primers (F: TTTCTGTTGGTGCTGATATTGCAGRGTTYGATYMTGGCTCAG, and R: ACTTGCCTGTCGCTCTATCTTCRGYTACCTTGTTACGACTT) were used to amplify the full-length rRNA gene from the isolated microbiome DNA samples. The amplicons (~ 1500 bp) were barcoded by PCR using the ONT barcoding primers (SQK-16S024), followed by multiplexing, and ONT adapter ligation. The Final library pool was sequenced on the Mk IC platform using FLO-MIN106 flow cells. The final pool of 16S full-length ONT libraries included 13 samples (4 Spu, 3 Sps, 3 F1, and 3 F2). One of the three F2 samples failed in the sequencing, leaving only two F2 samples.

Metatranscriptomic library preparation and sequencing

We treated the metatranscriptome total RNA with the FastSelect 5S/16S/23S reagent (Qiagen, cat# 335921) to deplete the rRNA following the manufacturer protocol. We started with 12.5 µl total RNA, fragmented for 2 min at 89 ˚C, and eluted the fragment RNA in 18 µl RNase-free water. The rRNA depleted was used immediately to prepare the metatranscriptome libraries using the SEQuoia complete RNAseq kit (Bio-Rad cat# 17005726) following the manufacturer’s procedures without any modifications. The final sequencing libraries were assessed on the Fragment analyzer, pooled equimolarly, and sequenced on the NextSeq 500 sequencer using a mid-output kit. The sequencing run produced over 171 million reads.

Microbiome data analysis

Taxonomical classification and abundance

We carried out the primary analysis using the EPI2ME 16S workflow (ONT 16S workflow 2019) to generate the taxonomic classification and abundance of the 16S genes for all samples. The 16S workflow uses the Centrifuge algorithm (Kim et al., 2016) to map the full-length 16s reads to the NCBI 16S-18S database and classify them according to the taxonomy of their best matches. The obtained taxonomical matrix for all samples was then used for more in-depth analysis using R as described in the results section. The overall analysis pipeline (Figure 1C) is divided into three modules: data preprocessing, microbiome comparison, and core and hub microbiome analyses. The R packages and analysis tools utilized in this study pipeline are listed in Figure 1C and will be elaborated on in the results section.

Statistical analysis of alpha and beta diversity

We utilized the Alpha diversity test to evaluate microbial diversity within each group, that is, among replicates of the same genotype (4 Spu, 3 Sps, 3 F1, and 2 F2). The alpha diversity was assessed using three methods: Chao1 (Chao, 1984), Pielou’s Evenness (Volvenko, 2014), and Shannon (Shannon, 2001). The pairwise Wilcoxon test (Wilcoxon, 1992) was then used to test the significance of the variance of the alpha diversity among genotypes. The difference in the microbiome composition among genotypes was calculated using beta diversity analysis, followed by visualization using PCoA (Principal Coordinate Analysis) (Gower, 1966). The sample distances were calculated using Bray-Curtis distances (Bray and Curtis, 1957; Kers and Saccenti, 2022).

Ubiquitous, genotype-enriched, and core microbiomes

We used an abundance-based approach to identify ubiquitous microbiomes, which is the microbiome existing in all samples above the background noise. This was followed by an overlap analysis using the ggvenn R package (Yan, 2021) to identify the genotype-enriched microbiomes. As for the core microbiome, we used the microbiome R package (Lahti and Shetty, 2018) to identify the core pitcher microbiome after optimizing the prevalence and limit of detection factors. We then employed the Linear Discriminant Analysis Effect Size (LEfSe) (Segata et al., 2011) method to identify the most distinctive microbial genera showing significant differential abundance between each pair of genotypes. The LefSe method is powerful because it combines statistical significance with the evaluation of biological consistency (effect size) using several tools including the Kruskal-Wallis test (Kruskal and Wallis, 1952) and Linear Discriminant Analysis (LDA) (Fisher, 1936).

Network and hub microbiome

We implemented network analysis to gain insights into the structure of the pitcher microbiome community. We deployed ggClusterNet (Wen et al., 2022) to infer the correlation network followed by visualization using Gephi (Bastian et al., 2009). Subsequently, we computed three scores for each node: hub score, within-module connectivity Z-score, and participation coefficient scores. The scores reflect the node’s connectivity to all other nodes in the network, to nodes within the same community, and the evenness of edges distribution within each community, respectively. The network structure, and hence all three node scores, will change when the number of nodes, the threshold for correlations, and the p-values are altered. Therefore, we examined the changes in the network structure across various correlation criteria, p-values, and node counts (Supplementary Figure S1). As a result, we set the correlation criterion to 0.7 and the p-value cutoff to 0.05.

Metatranscriptome data analysis

Metatranscriptome assembly and annotation

The metatranscriptome analysis pipeline is shown in Figure 1D. First, we trimmed the sequencing reads to remove adapter traces as well as short and low-quality reads using Trim_galore. Then we used the rRNA database v4 and sortMeRNA (Kopylova et al., 2012) to remove rRNA reads from each sample. The remaining metatranscriptomic mRNA reads from all samples were combined and assembled using rnaSPAdes (Bushmanova et al., 2019). The assembled metatranscriptome was deduplicated to remove identical metatranscripts followed by annotation against a customized uniprot_sprot protein database, which included bacterial, fungal, viral and archaea sprot sequences (uniprot_sprot_microbial) using diamond blastx (–top 1 –evalue 0.05). Also, we mapped the metatranscriptome to the NCBI nonredundant and GTDB protein databases using diamond blastx to obtain accurate taxonomical lineages. Subsequently, we mapped the cleaned sequence reads of each sample to the annotated metatranscriptome using STAR (Dobin et al., 2013) and compiled the count matrix for all samples with taxonomical information. Lastly, we converted the count matrix to CPM matrix.

Abundance of genera and metaproteins

In order to determine the abundance of taxa, the count matrix was transformed from reads per protein to genus abundance. The genus abundance was determined by counting all metatranscripts mapped to proteins from the same genus (according to the uniprot_sprot annotation) for the Spu, Sps, F1, and F2 samples. Next, we calculated the z-core and p-value for each genus count in each sample and filtered out all genera with p-values > 0.05. The TMM method was used to normalize the abundances of the filtered genera. Similarly, in order to identify prevalent proteins, the metatranscripts count matrix was transformed into a metaprotein abundance matrix. The metaprotein abundance was determined by counting the metatranscripts mapped to the same protein, based on the uniprot_sprot annotation, from multiple genera. Next, we calculated the z-core and p-value for each metaprotein count in each sample and filtered out all metaproteins with p-values > 0.01, followed by a TMM normalization.

Functional core microbiome

We used two criteria to identify the functional core microbiome: prevalence > 0.4 and limit of detection > 50 metatranscripts. Unlike the structural core microbiome, which was identified based on the 16S rRNA gene data, the functional core microbiome relies on the prevalence and abundance of protein-coding metatranscripts, hence their functions.

Genotype functional-enrichment analysis

We conducted function category enrichment analysis on the proteins that were enriched in each genotype (p-value < 0.05) using ShinyGO 0.80 (Ge et al., 2020) and Acidobacteria Bacteriam 13_1_20CM_53_8STRINGdb as a reference. The pathway databases of the enrichment analysis included GO biological processes, GO cellular components, GO molecular functions, annotated UniProt keywords, and local network clusters (STRING). We then manually filtered the list of enriched functions to remove redundancy, which existed as a result of using multiple databases.

Overall microbiome functional analysis

We used the intensively curated Plant Growth-Promoting Traits (PGPT) (Patz et al., 2024) from the Plant-associated Bacteria (PLaBAse) project (Patz et al., 2021) to gain insights into the major functions of the pitcher microbiome community. The metatranscriptome was mapped to the mgPGPT protein database using the Diamond blastx (mini. E-value ≤ 0.05) and used the Megan software (Gautam et al., 2023) to further analyze and visualize the results.

Results

Microbiome sequencing data

The ONT sequencing run produced 9,772,540 full-length 16S rRNA gene sequences of which 9,742,538 reads were taxonomically classified (~ 99%) with an average classification accuracy of 93%. The read length average and mode are 1,244bp, and 1,780bp, respectively. The quality score average and mode are 13.91, and 14.05 respectively. Reads without valid barcodes were removed leaving a total of 8,999,370 sequences. The EPI2ME 16S workflow assigned a lowest common ancestor (LCA) score (Munch et al., 2008) to each full-length 16S read (Supplementary Table S1). The LCA score is used instead of the classical “num_genus_taxid” and is used to store records of the classification status. It has three values: 1 for a species-level classification without accurate genus classification; 0 for both species and genus classification; -1 for unsuccessful species classification. We removed 106,817 reads that had a -1 LCA score. The remaining 8,892,553 reads (minimum classification accuracy is 77%) were used in the subsequent analyses.

Microbiome diversity and intergenotype variance

The use of full-length 16S rRNA genes in this analysis offers a heightened degree of precise classification analysis, hence eliminating the necessity for data rarefaction. The non-rarefied count matrix was further normalized as Counts Per Million (CPM) and was utilized to create the phyloseq (McMurdie and Holmes, 2013) object for further examination in the subsequent analysis.

Alpha diversity

The alpha diversity is measured using Chao1 (Chao, 1984), Pielou’s Evenness (Volvenko, 2014), and Shannon (2001) methods (Figure 2A). Each method measures particular aspects of microbial diversity. Chao1 assesses the species richness, Pielou’s Evenness focuses on the evenness of species distribution, and Shannon measures both richness and evenness. All three measures reveal noticeable variations in microbial richness and evenness between the genotypes (intergenotype) and among the replicates of the same group (intragenotype). When both richness and evenness are assessed simultaneously using Shannon’s approach, the median alpha diversities of all samples were more closely aligned. The observed variations may be attributed to the limited number of replicates, variations in sequence depth among samples, genetic differences among the replicates of the F1 and F2 samples, as well as the differences in microbial diversity. We then used the pairwise Wilcoxon method (Wilcoxon, 1992) to test if the alpha diversity among the parents, F1 and F2 is significantly different. All of the p-values were close to 1 (Supplementary Table S1), indicating that none of the pairs’ diversity showed a significant difference. This suggests that the low replication and sequencing depth variation did not result in any statistically significant difference in the diversity and species evenness observed across the samples. In the subsequent analysis, we assume there is no significant difference in the alpha diversity among the parents, F1, and F2 generations.

Figure 2

Figure 2. Microbiome diversity assessment. (A) Comparison of alpha diversity of different genotypes using Chao, Pielou, and Shannon tests to assess richness, evenness, and both, respectively. (B) Visualization of the relationship between every sample using Principal Coordinate Analysis. The color represents the genotype of each sample. The ellipses mark the variation of each genotype.

Beta diversity

In Figure 2B, we illustrate the beta diversity among genotypes using the PCoA (Principle Coordinate Analysis) (Gower, 1966), where the colored ellipse depicts the dispersion of all samples from each genotype, encompassing 95% of the diversity, assuming a normal distribution. Since there are only two replicates in the F2, it is infeasible to construct an ellipse for it, and hence it uses a straight line to indicate the dispersion. The results show that the microbiome compositions of the two parent species, Spu and Sps, are different from the descendant generations, F1 and F2. To further explore the inter-genotype microbiome differences, we implemented PERMANOVA (Anderson, 2001) to test the difference in the distribution between genotypes. The p-value for the test is 0.011, which indicates significant differences in microbial community structure between genotypes. We also implemented the PERMDISP test (Anderson, 2006) to assess inter-genotypes microbiome dispersion. The p-value of the test is 0.825, indicating the microbiome dispersion is not significantly different between genotypes. Combining the results from the two tests, we can conclude that there is a significant difference in the microbial community structure between the parental species and their F1, and F2 genotypes.

The genotype-enriched and core microbiomes

We identified the genotype-enriched and core microbiomes in the pitcher fluids of the parents, F1, and F2 groups. Specifically, a genus is deemed genotype-enriched only if it is represented in every replicate of that genotype by at least one full-length 16S rRNA read. Accordingly, out of 1,960 identified bacterial genera, the pitcher microbiomes of Spu, Sps, F1, and F2 contained 180, 177, 213, and 155 bacterial genera, respectively (Table 1; Supplementary Table S1). There is a total of 342 non-redundant bacterial genera constituting both the ubiquitous (62), and genotype-specific (145) bacterial taxa, as well as taxa existed in two to three genotypes (135). We utilized the ggvenn R package (Yan, 2021) to identify the microbiome overlap among genotypes (Figures 3A, C). Specifically, 62 (18.1%) bacterial genera exist in all samples of Spu, Sps, F1, and F2, among which the top 5 abundant genera are Aeromonas, Rhodopseudomonas, Achromobacter, Paraburkholderia, and Azospirillum.

Table 1

Table 1. summary of the results.

Figure 3

Figure 3. The core and genotype enriched microbiome analysis. (A) A Venn diagram showing the number of overlapped genera among Spu, Sps, F1, and F2. (B) Heatmap of the prevalence of different microbiomes under different detection thresholds. The detection value and prevalence cut-off limits are represented by the horizontal red line and transparent area, respectively. (C) Four Venn diagrams comparing the overlapping genera among every three genotypes to highlight the genotype specific numbers. (D) Genotype-enriched genes detected by Linear discriminant analysis Effect Size.

For the core microbiome, we considered a genus to be a member of the core microbiome if its least prevalence under different detection thresholds is greater than 40%, as described in the microbiome analysis R package (Lahti and Shetty, 2018). Unlike the ubiquitous microbiome, the core microbiome requires a specific abundance (the detection threshold) in a specific proportion of samples (prevalence) instead of in every sample. The heatmap in Figure 3B illustrates the prevalence of the core microbiome under different detection thresholds (100 to 1,000). With a prevalence equal to 0.4 and a detection limit equal to 500, a total of 58 bacterial genera constituted the core microbiome of Spu, Sps, F1, and F2 (Supplementary Table S1). The top 5 genera with the highest CPM normalized count in the core microbiome are Azospirillum, Achromobacter, Rhodopseudomonas, Mucilaginibacter, and Aeromonas. There are 42 genera that are shared by the ubiquitous and the core microbiomes, representing 67.7% and 72.4% of each, respectively (Supplementary Table S1).

In the genotype-enriched microbiome, 35 genera (10.2%) are unique to the Spu genotype, 23 genera (6.7%) to the Sps genotype, 68 genera (19.9%) to the F1 genotype, and 19 genera (5.6%) to the F2 genotype (Figure 3A). The results show that the microbiome of the F1 genotype differs significantly from that of the other genotypes. A similar observation can be obtained from the comparisons presented in Figure 3C. The microbiomes of the two parents (Sps and Spu) exhibit a greater similarity to each other than to either the F1 or F2 microbiomes (Figure 3A). The proportion of shared genera between Sps and Spu microbiomes is 37.8% when compared to F1 and 44.5% when compared to F2 microbiomes, which is much higher than the overlap between F1 and F2 microbiomes (Figure 3C). The overlaps between F1 and F2 microbiomes are 35.4 relative to the Spu and 36.8 relative to the Sps microbiomes. The results confirm that the Sps microbiome is similar to that of Spu, while the F1 microbiome is more closely related to the F2 microbiome (Figure 2B).

We employed LEfSe analysis (Segata et al., 2011) to examine the microbiomes of Spu, Sps, F1, and F2 to identify the most distinctive microbial genera showing significant differences in abundance across all genotype. Figure 3D displays the LEfSe analysis results for the comparisons Sps vs. F1, Spu vs. F1, and Spu vs. Sps. The F2 microbiome was excluded from the analysis due to having only two replicates, as previously explained. At a log10 LDA score of 2, there are more distinct genera in the microbiomes of F1 compared to Spu and Sps, than between Spu and Sps (Figure 3C). These results are consistent with the PCoA analysis (Figure 2B) and microbiome overlapping analysis (Figures 3A, C).

The microbiome communities networks and hubs

This analysis included 582 genera (out of 1,960) with an accumulative normalized abundance of more than 100 across all samples. In the correlational network (Figure 4A), nodes with correlation ≥ 0.7 and p-value ≤ 0.05 are connected with edges. The greedy algorithm (Clauset et al., 2004) assigned the network nodes to 5 microbiome communities (clusters), which are indicated by the node and edge colors (Figure 4A-left). The edge colors represent the correlation values of connected nodes, where blue and red edges signify positive and negative correlations, respectively. Two communities contain less than 10 genera (Cluster 1, 2), and the majority of the genera (97%) are assigned to the other three communities (Cluster 3, 4, 5). Cluster 1 (dark green) is located in the negatively correlated part of the network and serves as the hub node for the other three major microbiome communities. Cluster 2 (orange) has very few connections to the rest of the network. Cluster 3 (pink), Cluster 4 (light green), and Cluster 5 (blue) construct the main part of the microbiome network. The top 5 genera with the highest normalized count in these five communities are: Sphingobacterium, Pedobacter, Tissierella, Epilithonimonas, and Diaphorobacter for Cluster 1; Gracilibacter, Desulfonatronobacter, Desulfococcus, Maritalea, and Desulfoconvexum for Cluster 2; Azospirillum, Mucilaginibacter, Legionella, Terriglobus, and Methylacidimicrobium for Cluster 3; Achromobacter, Rhodopseudomonas, Aeromonas, Comamonas, and Burkholderia for Cluster 4; Sphingomonas, Luteibacter, Granulicella, Reyranella, and Terrimonas for Cluster 5.

Figure 4

Figure 4. Network analysis. (A) the whole network of 582 filtered genera (left) and connectivity scores (Right). (B-E left) the subnetworks of top 100 genera ranked by their hubscore (B), z-score (C), normalized abundance (D), and participant coefficient (E). The connectivity scores of the four subnetworks are in their corresponding right panels.

We also computed the within-model connectivity z-score and participation coefficient (Guimera and Nunes Amaral, 2005) for each node in the whole network, which indicates the connectivity within and between communities, respectively (Figure 4A-right). Nodes are classified into four categories based on their scores: peripheral, connections, module hubs, and network hubs (Guimera and Nunes Amaral, 2005). A network hub node has a stronger connection to the network than the ordinary node. Connector and module hub nodes have more connectedness across and within communities than the average of community nodes. The remaining genera are classified as peripheral. As shown in Figure 4A-right, the whole network has 2 module hubs and 10 connectors. The module hubs are Schlesneria and Methylophaga, whereas the connectors are Dendronalium, Calothris, Lacrimispora, Sphingobium, Mesorhizbobium, Parabacteroides, Legionella, Leifsonia, Lactococcus, and Desulfococcus.

To gain a closer insight into the microbiome communities, we zoomed in on the top 100 genera ranked by their hubscore, within-model connectivity z-score, normalized abundance, and participation coefficient (Figures 4B–E). In the co-occurrence subnetworks (Figures 4B–E, left), node colors denote the greedy algorithm-detected clusters and node sizes represent the accumulated abundance in all 12 samples. The connectivity scores of the four subnetworks are illustrated in Figures 4B–E (right panel), respectively. The inferred networks of the top 100 genera ranked by hubscore and within-model z-score have similar topology including 3 highly connected community clusters each, two module hubs (Schlesneria and Methylophaga), and no connectors or network hubs (Figures 5B, C). The high connectivity in these two subnetworks can be explained by the fact that both the hubscore and within-model connectivity z-score prioritize the connections of each node, resulting in a highly interconnected network with fewer distinct clustering patterns. Similarly, the inferred subnetworks of the top 100 genera ranked by normalized abundance and participation coefficient (Figure 4E, E) have similar structures with 6 and 7 clusters, and Phenylobacterium and Schlesneria as module hubs, respectively. The abundance-based subnetwork has the following 11 connector nodes: Phenylobacterium, Schlesneria, Chitinophaga, Rhodoferax, Rhizobium, Flavobacterium, Dokdonella, Caballeronia, Gluconacetobacter, Calothrix, Erythobacter, Novosphingobium, Arcicella, and Entrobacter (Figure 4E-right). The abundance and participation coefficient subnetworks seem to capture more of the network topological information and identify more community connectors and hubs.

Figure 5

Figure 5. Genera overlapping among subnetworks and abundance of hub genera. (A-C) Venn diagrams showing the number of overlapped nodes (genera) among the four subnetworks. Overlapping between the hubscore and z-score top most 100 genera (A), normalized abundance and participant coefficients (B), and all four subnetworks (C). (D) A heatmap showing the normalized abundance of all hub and connector genera in the Spu, Sps, F1, and F2 samples.

The overlap among the four subnetworks is presented in Figure 5. A total of 67 genera (50.4%), including a module hub (Schlesneria), appeared in both hubscore and z-score subnetworks (Figure 5A). On the other hand, only 22 genera (12.4%) are common to both the abundance and participation coefficient subnetworks (Figure 5B). Interestingly, no module hubs or connectors were shared between these two subnetworks. Only two genera overlapped among the four subnetworks (Figure 5C). The participation coefficient quantifies the overall connectivity of each node, making it a valuable tool for capturing the structure of the network of the top 100 genera. While there is minimal overlap between the genera in the abundance and participation coefficient subnetworks (Figure 5B), both networks have a similar topology, suggesting that these two parameters represent comparable community clusters. The heatmap in Figure 5D illustrates the TMM normalized abundance of the discovered hubs and connector genera in the parents, F1, and F2 data. It reveals three clusters and six subclusters. Also, it reveals the presence of a cluster of extremely prevalent hub and connector nodes in each of the genotypes, indicating that the genotype of the host may impact the formation of the hubs and connectors within the microbiome community.

Characterization of the pitcher metatranscriptome

The assembled pitcher metatranscriptome included a total of 132,710 metatranscripts, with an N50, average and maximum length of 211 bp, 219 bp, and 8048 bp, respectively. About 65,578 (49.4%), 62,382 (47%), and 63,105 (47.6%) metatranscripts have blastx protein hits (e-value ≤ 0.05) in the uniprot_sprot_microbial, NCBI nonredundant, and GTDB protein databases, respectively. Figure 6A shows a random example of the blastx alignment of metatranscripts to reference protein (min. coverage ≥ 70%). The total number of non-redundant proteins is 9693. These best protein hits come from 2157 microbial species representing 121 viral, 1067 bacterial, 11 archaeal and 270 fungal lineages (Table 1; Supplementary Table S2). At the microbial level, the highest represented viral clade is Riboviria (71 metatranscripts), the highest bacterial taxon is proteobacteria (17,517 metatranscripts), and an unclassified Archaea (3975 metatranscripts). About 9021 metatranscripts mapped to proteins from eukaryotic organisms, the highest among them were Cryptophyceae (432), Discoba (302), Opisthokonta (503), and Sar (4765) (Figure 6B). The metatranscriptome data revealed a functional core microbiome consisting of 90 microbial species, including 68 bacterial taxa and 22 fungal taxa (Supplementary Table S2). Interestingly, sixteen bacterial taxa are part of both the structural and functional microbiomes (Supplementary Table S2). The taxonomy tree generated from the metatranscripts (Figure 6B) offers a visual representation of the general classification of the microbial taxa in the pitcher microbiome. It is worth noting that there can be heterogeneity in the microbiome structure from one individual to another. However, as a whole, the pitcher microbiome appears to have a very intricate structure with multiple levels of interactions.

Figure 6

Figure 6. Taxonomy tree of the identified taxa based on the pitcher metatranscriptome alignment to proteins. (A) Metatranscripts of variable lengths aligned to the protein MAG TPA: carbamoyl-phosphate synthase large subunit (Sphingobacterium sp.). The vertical gray lines represent gaps. (B) The taxonomy tree of the identified taxa. The circles represent the log scale of number of metatranscripts assigned to the tree nodes and leaves. The numbers listed after each taxon name is the number of metatranscripts mapped to the proteins of this taxon.

Genotype-based differences in the pitcher metatranscriptomes

We examined the differences in the pitcher metatranscriptome among the genotypes (Spu, Sps, F1, and F2) using the genera and metaproteins enrichment analysis. For each genotype, we selected the genera with p-value < 0.05 and metaproteins with p-values < 0.01. There were 25 genera and 212 metaproteins meeting these criteria. The normalized abundance of these genera and proteins is illustrated in the heatmaps in Figure 7. It is evident from the heatmap clustering that there are differences in the genera abundance (Figure 7A) and metaprotein accumulation (Figure 7B) among the genotypes. Note that in Figure 7B, the heatmap includes only 25 proteins, but the complete list of proteins is provided in Supplementary Table S2. Despite the intergenotypic difference in the genera diversity and metaproteins, the enriched functional categories in all genotypes seem to be very similar, with small variations in the fold changes and number of genes in each category (Figure 8). For each of genotype, proteins with p-values < 0.01 (see materials and methods) were used in the functional category enrichment analysis. These sets included 199, 237, 236, and 201 proteins in Spu, Sps, F1, and F2 respectively. The analysis revealed a significant (FDR < 0.01) enrichment (more than two folds) of multiple functional categories in all genotypes. Nevertheless, there were no discernible differences in the enriched categories across the different genotypes (Figure 8), suggesting a microbial function convergence in the pitcher microcosm despite variations in structure and diversity.

Figure 7

Figure 7. The pitcher metatranscriptome: the abundant genera and proteins showing species-related patterns. (A) heatmap showing the TMM normalized abundance of the genera in the samples of Spu, Sps, F1, and F2 (p-values < 0.05). (B) similar to “A” but showing the abundance of the 25 protein with lowest P-values.

Figure 8

Figure 8. Function enrichment of the genotype-enriched metaproteins against the Go functional categories, uniprot keywords, STRING functional lists of Acidobacteria. The number of genes, FDR and fold enrichment of the significantly enriched functional categories in Spu, Sps, F1, and F2 metatranscriptomes are shown in (A-D), respectively.

Phyllosphere and rhizosphere microbiomes

We compared the identified pitcher microbial taxa with previously reported phyllosphere and rhizosphere microbiomes (reviewed in (Saeed et al., 2021; Bashir et al., 2022)). As presented in Table 2, groups of 30, 23, and 20 microbial taxa that were reported to be performing different functions in the phyllosphere, rhizosphere, and both, respectively, were also present in the pitcher microbiomes in this study. Furthermore, many of these microbes were community connector hubs, indicating that they play a major role in microbiome assembly and functions (Table 2). The richness of the pitcher microbiome with both phyllosphere and rhizosphere microbes indicates its intricate nature and suggests that its functions extend beyond the traditional role of supporting plant nutrition.

Table 2

Table 2. The phyllosphere (PS) and Rhizosphere (RS) bacterial (B) and fungal (F) taxa in the pitcher microbiome and metatranscriptome.

The pitcher metatranscriptome functions

Our analysis of the pitcher metatranscriptome using the PGPT database revealed major function categories (Table 1; Figure 9). About 50,424 metatranscripts were mapped to PGPT proteins that are associated with different plant growth-promoting functional categories. Globally, 13,519 and 29,028 metatranscripts are mapped to functional groups with direct and indirect effects on plant growth traits, respectively (Figure 9). Among the microbiome traits that have direct effects on plant growth are biofertilization, bioremediation, and phytohormonal plant signals. The microbiome traits that have indirect effect on plant growth include stress control, immune response stimulation, colonizing plant systems and microbial competitive exclusion. As shown in Figure 9, the size of the green circle represents the log scale of the number of metatranscripts assigned to the corresponding function. Altogether, there are 44 major microbiome traits that have direct effect on the host-microbiome and microbiome-microbiome interactions. Traits that are represented by over a thousand metatranscripts include nitrogen acquisition (1003), phosphate solubilization (1751), heavy metals detoxifications (2413), plant vitamin production (1069), neutralizing biotic stress (1578), neutralizing abiotic stress (6230), universal stress response (1062), chemotaxis (1033), plant-derived substrate usage (7067), colonization surface attachment (1485), cell envelop remodeling (1618), biofilm formation (2313), bacterial fitness (2526), and bacterial secretion (1188). These findings reveal the complex functions of the pitcher microbiome at several levels, providing a chance to investigate a wide range of microbiome traits and their host interactions.

Figure 9

Figure 9. Classification of the major function categories of the pitcher microbiome based on the number of metatranscripts mapped to PGPT ontology. The circles represent the log scale of number of metatranscripts assigned to each function. The number following each function is the number of metatranscripts mapped to this function.

Discussion

The Sarracenia genus relies on its microbiome to facilitate the digestion of trapped prey and to provide important nutrients such as nitrogen (N) for protein synthesis, phosphorus (P) for nucleic acid synthesis, and magnesium (Mg) and iron (Fe) for chlorophyll synthesis. Previous research suggested that factors other than prey capture and colonization by eukaryotic species may affect the recruitment of bacteria to the pitchers (Grothjan and Young, 2019). Also, It has been demonstrated that plants have the ability to modify the composition of the microbiomes in their surroundings (Mahnert et al., 2015) and host genetic variations have an impact on the structure and functions of the microbiome in humans (Turnbaugh et al., 2009; Blekhman et al., 2015). Extensive research has been conducted on the convergent evolution of carnivorous plant leaves, particularly those belonging to the Sarracenia genus (Thorogood et al., 2018; Wheeler and Carstens, 2018; Chomicki et al., 2024). The evolution of the composite trait of the insect-trapping pitchers in carnivorous plants is believed to have occurred through separate modifications of shape, coloration, and biosynthetic pathways (reviewed in (Chomicki et al., 2024)). We are exploring the concept of considering the pitcher microbiome properties as an extension of plant traits that have evolved alongside carnivory as a part of the insect-trapping composite trait. Our hypothesis is that carnivorous plants have gained genetic traits during the convergent evolution of their pitcher trait, enabling them to interact with their vital microbiomes. Our results show that pitcher plants impact the structure and function of their microbiomes, which in turn have a broader effect on their growth, extending beyond mere nutritional support.

The host genotype significantly influence the microbiome assembly

In their natural habitats, the microbiome diversity diverged among the Sarracenia host species (Heil et al., 2022). In this study, we detected statistically significant differences in the structure of the pitcher microbiome communities among the parental species and their F1 and F2 genotypes under the same greenhouse conditions. The finding is supported by the PERMANOVA analysis, which indicated that the variations in distribution among genotypes are statistically significant (p-value = 0.011). The PERMDISP test indicated that there is no significant difference in the dispersion (p-value = 0.825) between genotypes, ruling out the possibility that difference in sequencing depth between samples could have produced the observed distribution differences. Both Alpha-diversity and PCoA (Figure 2) analyses unveiled differences in the microbiome diversity among individuals of the same group and between the genotypes. Although the majority of the microbial taxa existed in the microbiomes of all species and genotypes, there were genotype-enriched microbial taxa (Figures 3A–D). The microbiome of the F1 generation differed significantly from all other samples, which could be attributed to its dominant genetics that are expressed in more distinctive traits. The variation in the microbiota between the parental species is smaller than the variation observed between either of them and their F1 and F2 genotypes. This observation is reinforced by the number of distinctive genera, as identified by the LEfSe analysis, between the parental species and their F1 and F2 generation (Figure 3C). The F2 generation exhibits increased genotype complexity and genetic variation, resulting in a wide range of traits, including those that influence the assembly and structure of their microbiome. The abundance patterns of the overall microbiome community hub and connector genera varied significantly among the different genotypes (Figures 4, 5). This suggests that the genotypes may influence the development of their respective local community hubs and connectors. By utilizing the Sarracenia mapping population, our results suggest that host genetic factors are impacting the observed significant variance in the microbiome structure.

The core microbiomes and plant growth-promoting traits

The richness and diversity of the Sarracenia core microbiome suggests that it possesses a wide range of complex traits and levels of interactions. Over 1900 bacterial genera were identified using 16S the rRNA gene, including 62 ubiquitous genera, 145 genotype-enriched genera, and 58 genera constituting the structural core microbiome. More than 2100 microbial genera were found using the metatranscriptomic data, including 68 bacterial and 22 fungal genera constituting the functional core microbiome. There are 16 genera existed in both structural and functional core microbiomes, including Acidovorax, Bradyrhiobium, Bukholderia, Caulobacter, Chromobacterium, Clostridium, Cupriavidus, Magnetospirilllum, Mesorhizobium, Novosphingobium, Paraburkholderia, Pseudomonas, Rhizobium, Rhodopesudomonas, Rickettsia, and Stenotrophomonas. The functional core microbiome relies on the prevalence and abundance of protein-coding metatranscripts, while the identification of the structural microbiome relies on the prevalence and abundance of 16S rRNA gene reads. Many of the core genera especially the prevalent ones are among the group of plant-associated microbes known as plant growth-promoting Rhizobacteria (PGPR) (Beneduzi et al., 2012; Agarwal et al., 2020). As shown in Table 2, we identified the fungal and bacterial taxa in the pitcher microbiome, which have previously been discovered in rhizospheres and/or phyllosphere of other plant species (reviewed in (Saeed et al., 2021; Bashir et al., 2022)). While PGPR has traditionally been found in the rhizosphere, it is not unexpected to discover them in the pitcher microbiome (phyllosphere) due to the fact that this plant species depends on pitcher leaves for acquiring essential nutrients. This finding prompts inquiries into the functional and structural interactions between the rhizosphere and phyllosphere microbiomes while coexisting in the confined microenvironment of the pitchers. The functional traits of these PGPR genera include the ability to enhance plant growth through biofertilization, act as biocontrol agents to combat pests, and serve as biological fungicides to protect their specific host plants.

Microbiome community clusters and hubs

The network analysis (Figures 4, 5) unveiled many clusters of microbial communities (Supplementary Table S1), resulting in the identification of community hubs and connector genera. The number of community clusters varied between 3 and 7, depending on the network factors used to identify the clusters. For example, the analysis revealed the presence of five distinct community clusters in the entire network. Additionally, the networks of the top 100 abundant and top 100 participation coefficient nodes (genera) exhibited the identification of six and seven community clusters, respectively. But regardless of the method, there are always distinctive community clusters in the inferred microbiome network. The key characteristic of these communities is the discovery of a total of 27 genera playing the hub and connector roles in these communities. It is crucial to note that the hub and connector genera, based on hubscore, within-model connectivity z-score, and participation coefficient, have a low abundance. This indicates that the genera that perform important roles in the microbiome community are not always the most prevalent genera. Interestingly, the analysis of the normalized abundance of the 27 hub and connector genera in the parents, F1, and F2 genotypes uncovers distinct clusters of genera that vary significantly in their abundance among the genotypes (Figure 5D). This suggests that the host genetic makeup influences the selection of the microbiome hubs and connectors. For instance, the genus Schlesneria appeared as a hub genus in the whole network, as well as in the subnetworks, despite not being among the highly abundant genera compared to the other hub and connectors genera (Figure 5D; Supplementary Table S1). As shown in Table 2, the majority of the hub and connector genera are commonly reported as plant-associated genera living in the plant’s phyllosphere and rhizosphere (Mahnert et al., 2018). Overall, the inferred microbiome networks offered insights into the assembly and organization of the pitcher microbiome. This involved the identification of many subcommunities, each with their hub genera, as well as connector genera that facilitate communications between these subcommunities. Interestingly, the community networks also show that the hub and connector genera exhibit varying levels of abundance across different genotypes (Figure 5D), suggesting that the host genotypes may have an impact on their development.

The complex function of pitcher microbiome

Although our findings demonstrate a considerable influence of the host on the diversity and organization of the microbiome, we did not discover any notable genotype-related variations in the enriched functional categories in the metatranscriptome (Figure 8). This analysis, however, lacks the level of detail necessary to discover any fine functional differences. Nevertheless, it suggests that even though there may be fine difference in the microbiome functions among genotypes, the overall microbiome traits and functions remain consistent. This could be explained by functional redundancy among the taxa of the pitcher microbiome and a microbiome functional convergence to achieve similar functional traits (Louca et al., 2016, 2018). A similar conclusion was reached by analyzing the microbial communities and functions in Spu populations in natural habitats (Grothjan and Young, 2019). The functional analysis, however, revealed complex and multifaceted functional traits of the pitcher microbiome (Figure 9). At the host-microbiome interaction level, approximately 13.5% of the metatranscripts that have been classified functionally are related to genes involved in biofertilization activities. These activities encompass nitrogen acquisition, phosphate solubilization, and iron acquisition. The presence of a significant portion of the pitcher microbiome that contributes to plant biofertilization is expected, given that acquiring nutrients from leaves is a key characteristic associated with the carnivorous behavior. Nevertheless, the discovery of the pitcher microbiome’s bioremediation (6%) and phytohormonal plant signal generation (7%) activities are novel additions to our understanding of the Sarracenia-microbiome interaction. The microbiome has a significant role in the detoxification of heavy metals (2413 metatranscripts) and the degradation of xenobiotics (607 metatranscripts) in the pitcher fluid, which serves to safeguard the host plant. These functions directly impact the fitness and defense of the pitcher plant and can be considered extended plant traits. Approximately 58% of the pitcher metatranscriptome is assigned to a set of functions categorized as traits that have an indirect impact on the host plan (Figure 9). The most prevalent functions, as indicated by the abundance of metatranscripts, are those associated with stress regulation, promotion of plant immune response, and colonization of plant systems. On the microbiome-microbiome interaction level, the analysis revealed that major portions of the metatranscriptome are involved in microbial competitive exclusion activities (8376 metatranscripts) and usage of plant-derived substrates (7067 metatranscripts).

These results offer a comprehensive understanding of the functions of the microbiome in connection to the host plant, microbial community, and environmental conditions inside the pitcher microcosm. Furthermore, they highlight some strategies utilized by the host plant to regulate the assembly of the microbiome, such as plant immunity and the release of extracellular substrates. The pitcher contains a varied microbial population with complicated multiscale functions, which can be used for metagenomic mining to discover new beneficial microbial taxa.

Conclusions and prospects

The variations in the assembly and structure of the microbial communities within the pitchers of the Sarracenia mapping population can be attributed, at least in part, to the genetic characteristics of the host. The pitcher microbiome performs multiple functions that have both direct and indirect effects on the host. Nevertheless, the precise genetic and biochemical pathways via which the host plant interacts with the microbiome are still unidentified. Our goal is to link the extended traits of the Sarracenia plant (namely, microbiome traits and functions) to the genetic elements of Sarracenia using cross-species QTL linkage mapping. This will enable us to gain insights into the mechanisms underpinning the interactions between the host and its microbiome, which is an essential requirement for the pursuits of microbiome engineering.

Data availability statement

The raw sequencing data are deposited in NCBI RSA under bioporjct PRJNA1163405. The RSA Ids of the datasets are 43849136 through 43849147.

Author contributions

JC: Data curation, Formal analysis, Software, Visualization, Writing – original draft, Writing – review & editing. IM: Investigation, Methodology, Writing – review & editing. WR: Investigation, Methodology, Writing – review & editing, Resources. MZ: Writing – review & editing, Formal analysis. LJ: Writing – review & editing, Resources. RM: Resources, Writing – review & editing. MA: Conceptualization, Funding acquisition, Resources, Writing –review & editing, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Supervision, Validation, Visualization, Writing – original draft.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2024.1445713/full#supplementary-material

Supplementary Figure 1 | Examining the changes in the network structure under various correlation, p-value, and node counts criteria.

Supplementary Table S1 | microbiome data.

Supplementary Table S2 | metatranscriptome data.

References

Adlassnig, W., Peroutka, M., Lendi, T. (2011). Traps of carnivorous pitcher plants as a habitat: composition of the fluid, biodiversity and mutualistic activities. Ann. Bot. 107, 181–194. doi: 10.1093/aob/mcq238

PubMed Abstract | Crossref Full Text | Google Scholar

Afridi, M. S., Ali, S., Salam, A., César Terra, W., Hafeez, A., Sumaira, et al. (2022). Plant microbiome engineering: hopes or hypes. Biol. (Basel) 11, 1782. doi: 10.3390/biology11121782

Crossref Full Text | Google Scholar

Agarwal, P., Giri, B. S., Rani, R. (2020). Unravelling the role of rhizospheric plant-microbe synergy in phytoremediation: A genomic perspective. Curr. Genomics 21, 334–342. doi: 10.2174/1389202921999200623133240

PubMed Abstract | Crossref Full Text | Google Scholar

Albert, V. A., Williams, S. E., Chase, M. W. (1992). Carnivorous plants: phylogeny and structural evolution. Science 257, 1491–1495. doi: 10.1126/science.1523408

PubMed Abstract | Crossref Full Text | Google Scholar

Anderson, M. J. (2001). A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26, 32–46. doi: 10.1111/j.1442-9993.2001.01070.pp.x

Crossref Full Text | Google Scholar

Anderson, M. J. (2006). Distance-based tests for homogeneity of multivariate dispersions. Biometrics. 62. doi: 10.1111/j.1541-0420.2005.00440.x

Crossref Full Text | Google Scholar

Baiser, B., Buckley, H. L., Gotelli, N. J., Ellison, A. M. (2013). Predicting food-web structure with metacommunity models. Oikos 122, 492–506. doi: 10.1111/j.1600-0706.2012.00005.x

Crossref Full Text | Google Scholar

Bashir, I., War, A. F., Rafiq, I., Reshi, Z. A., Rashid, I., Shouche, Y. S. (2022). Phyllosphere microbiome: Diversity and functions. Microbiol. Res. 254, 126888. doi: 10.1016/j.micres.2021.126888

PubMed Abstract | Crossref Full Text | Google Scholar

Bastian, M., Heymann, S., Jacomy, M. (2009). “Gephi: An Open Source Software for Exploring and Manipulating Networks.” in Proceedings of the International AAAI Conference on Web and Social Media. 3, 361–362. doi: 10.1609/icwsm.v3i1.13937

Crossref Full Text | Google Scholar

Beneduzi, A., Ambrosini, A., Passaglia, L. M. P. (2012). Plant growth-promoting rhizobacteria (PGPR): Their potential as antagonists and biocontrol agents. Genet. Mol. Biol. 35, 1044–1051. doi: 10.1590/S1415-47572012000600020

PubMed Abstract | Crossref Full Text | Google Scholar

Bittleston, L. S., Wolock, C. J., Yahya, B. E., Chan, X. Y., Chan, K. G., Pierce, N. E., et al. (2018). Convergence between the microcosms of Southeast Asian and North American pitcher plants. eLife 7, e36741. doi: 10.7554/eLife.36741.023

PubMed Abstract | Crossref Full Text | Google Scholar

Blekhman, R., Goodrich, J. K., Huang, K., Sun, Q., Bukowski, R., Bell, J. T., et al. (2015). Host genetic variation impacts microbiome composition across human body sites. Genome Biol. 16, 191. doi: 10.1186/s13059-015-0759-1

PubMed Abstract | Crossref Full Text | Google Scholar

Boyer, T., Carter, R. (2011). Community analysis of green pitcher plant (Sarracenia oreophila) bogs in Alabama. Castanea 76, 364–376. doi: 10.2179/10-048.1

Crossref Full Text | Google Scholar

Bray, J. R., Curtis, J. T. (1957). An ordination of the upland forest communities of Southern Wisconsin. Ecol. Monogr. 27, 325–349. doi: 10.2307/1942268

Crossref Full Text | Google Scholar

Bushmanova, E., Antipov, D., Lapidus, A., Prjibelski, A. D. (2019). rnaSPAdes: a de novo transcriptome assembler and its application to RNA-Seq data. GigaScience 8, giz100. doi: 10.1093/gigascience/giz100

PubMed Abstract | Crossref Full Text | Google Scholar

Chao, A. (1984). Nonparametric estimation of the number of classes in a population. Scandinavian J. Stat 11, 265–270. Available online at: https://www.jstor.org/stable/4615964.

Google Scholar

Chomicki, G., Burin, G., Busta, L., Gozdzik, J., Jetter, R., Mortimer, B., et al. (2024). Convergence in carnivorous pitcher plants reveals a mechanism for composite trait evolution. Science 383, 108–113. doi: 10.1126/science.ade0529

PubMed Abstract | Crossref Full Text | Google Scholar

Clauset, A., Newman, M. E. J., Moore, C. (2004). Finding community structure in very large networks. Phys. Rev. E 70, 066111. doi: 10.1103/PhysRevE.70.066111

Crossref Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | Crossref Full Text | Google Scholar

Ellison, A. M., Buckley, H. L., Miller, T. E., Gotelli, N. J. (2004). Morphological variation in Sarracenia purpurea (Sarraceniaceae): geographic, environmental, and taxonomic correlates. Am. J. Bot. 91, 1930–1935. doi: 10.3732/ajb.91.11.1930

PubMed Abstract | Crossref Full Text | Google Scholar

Ellison, A. M., Gotelli, N. J., Błędzki, L. A., Butler, J. L. (2021). Regulation by the Pitcher Plant Sarracenia purpurea of the Structure of its Inquiline Food Web. amid 186, 1–15. doi: 10.1674/0003-0031-186.1.1

Crossref Full Text | Google Scholar

Ellison, A. M., Gotelli, N. J., Brewer, J. S., Cochran-Stafira, D. L., Kneitel, J. M., Miller, T. E., et al. (2003a).The evolutionary ecology of carnivorous plants. In: Advances in Ecological Research (Academic Press). Available online at: https://www.sciencedirect.com/science/article/pii/S0065250403330090 (Accessed March 16, 2024).

Google Scholar

Ellison, A. M., Gotelli, N. J., Brewer, J. S., Cochran-Stafira, D. L., Kneitel, J. M., Miller, T. E., et al. (2003b). The evolutionary ecology of carnivorous plants. Adv. Ecol. Res. 33, 1–74. doi: 10.1016/S0065-2504(03)33009-0

Crossref Full Text | Google Scholar

Fisher, R. A. (1936). The use of multiple measurements in taxonomic problems. Ann. Eugenics 7, 179–188. doi: 10.1111/j.1469-1809.1936.tb02137.x

Crossref Full Text | Google Scholar

Fu, C.-N., Wicke, S., Zhu, A.-D., Li, D.-Z., Gao, L.-M. (2023). Distinctive plastome evolution in carnivorous angiosperms. BMC Plant Biol. 23, 660. doi: 10.1186/s12870-023-04682-1

PubMed Abstract | Crossref Full Text | Google Scholar

Furches, M. S., Small, R. L., Furches, A. (2013). Hybridization leads to interspecific gene flow in Sarracenia (Sarraceniaceae). Am. J. Bot. 100, 2085–2091. doi: 10.3732/ajb.1300038

PubMed Abstract | Crossref Full Text | Google Scholar

Gautam, A., Zeng, W., Huson, D. H. (2023).DIAMOND +  MEGAN Microbiome Analysis. In: Metagenomic Data Analysis (New York, NY: Springer US) (Accessed May 6, 2024).

Google Scholar

Ge, S. X., Jung, D., Yao, R. (2020). ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics 36, 2628–2629. doi: 10.1093/bioinformatics/btz931

PubMed Abstract | Crossref Full Text | Google Scholar

Gower, J. C. (1966). Some distance properties of latent root and vector methods used in multivariate analysis. Biometrika 53, 325–338. doi: 10.1093/biomet/53.3-4.325

Crossref Full Text | Google Scholar

Grothjan, J. J., Young, E. B. (2019). Diverse microbial communities hosted by the model carnivorous pitcher plant Sarracenia purpurea: analysis of both bacterial and eukaryotic composition across distinct host plant populations. PeerJ 7, e6392. doi: 10.7717/peerj.6392

PubMed Abstract | Crossref Full Text | Google Scholar

Grothjan, J. J., Young, E. B. (2022). Bacterial recruitment to carnivorous pitcher plant communities: identifying sources influencing plant microbiome composition and function. Front. Microbiol. 13, 791079. doi: 10.3389/fmicb.2022.791079

PubMed Abstract | Crossref Full Text | Google Scholar

Guimera, R., Nunes Amaral, L. A. (2005). Functional cartography of complex metabolic networks. nature 433, 895–900. doi: 10.1038/nature03288

PubMed Abstract | Crossref Full Text | Google Scholar

Heil, J. A., Wolock, C. J., Pierce, N. E., Pringle, A., Bittleston, L. S. (2022). Sarracenia pitcher plant-associated microbial communities differ primarily by host species across a longitudinal gradient. Environ. Microbiol. 24, 3500–3516. doi: 10.1111/1462-2920.15993

PubMed Abstract | Crossref Full Text | Google Scholar

Ke, J., Wang, B., Yoshikuni, Y. (2021). Microbiome engineering: synthetic biology of plant-associated microbiomes in sustainable agriculture. Trends Biotechnol. 39, 244–261. doi: 10.1016/j.tibtech.2020.07.008

PubMed Abstract | Crossref Full Text | Google Scholar

Kers, J. G., Saccenti, E. (2022). The power of microbiome studies: some considerations on which alpha and beta metrics to use and how to report results. Front. Microbiol. 12, 796025. doi: 10.3389/fmicb.2021.796025

PubMed Abstract | Crossref Full Text | Google Scholar

Kim, D., Song, L., Breitwieser, F. P., Salzberg, S. L. (2016). Centrifuge: rapid and sensitive classification of metagenomic sequences. Genome Res. 26, 1721–1729. doi: 10.1101/gr.210641.116

Crossref Full Text | Google Scholar

Koopman, M. M., Carstens, B. C. (2011). The microbial phyllogeography of the carnivorous plant sarracenia alata. Microb. Ecol. 61, 750–758. doi: 10.1007/s00248-011-9832-9

PubMed Abstract | Crossref Full Text | Google Scholar

Koopman, M. M., Fuselier, D. M., Hird, S., Carstens, B. C. (2010). The carnivorous pale pitcher plant harbors diverse, distinct, and time-dependent bacterial communities. Appl. Environ. Microbiol. 76, 1851–1860. doi: 10.1128/AEM.02440-09

PubMed Abstract | Crossref Full Text | Google Scholar

Kopylova, E., Noé, L., Touzet, H. (2012). SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics 28, 3211–3217. doi: 10.1093/bioinformatics/bts611

PubMed Abstract | Crossref Full Text | Google Scholar

Kruskal, W. H., Wallis, W. A. (1952). Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 47, 583–621. doi: 10.1080/01621459.1952.10483441

Crossref Full Text | Google Scholar

Lahti, L., Shetty, S. (2018). Introduction to the microbiome R package. Available online at: https://microbiomegithubio/tutorialshttps://s3.jcloud.sjtu.edu.cn/899a892efef34b1b944a19981040f55b-oss01/bioconductor/3.10/bioc/vignettes/microbiome/inst/doc/vignette.html (Accessed December 18, 2023).

Google Scholar

Louca, S., Jacques, S. M. S., Pires, A. P. F., Leal, J. S., Srivastava, D. S., Parfrey, L. W., et al. (2016). High taxonomic variability despite stable functional structure across microbial communities. Nat. Ecol. Evol. 1, 1–12. doi: 10.1038/s41559-016-0015

Crossref Full Text | Google Scholar

Louca, S., Polz, M. F., Mazel, F., Albright, M. B. N., Huber, J. A., O’Connor, M. I., et al. (2018). Function and functional redundancy in microbial systems. Nat. Ecol. Evol. 2, 936–943. doi: 10.1038/s41559-018-0519-1

PubMed Abstract | Crossref Full Text | Google Scholar

Mahnert, A., Haratani, M., Schmuck, M., Berg, G. (2018). Enriching beneficial microbial diversity of indoor plants and their surrounding built environment with biostimulants. Front. Microbiol. 9, 2985. doi: 10.3389/fmicb.2018.02985

PubMed Abstract | Crossref Full Text | Google Scholar

Mahnert, A., Moissl-Eichinger, C., Berg, G. (2015). Microbiome interplay: plants alter microbial abundance and diversity within the built environment. Front. Microbiol. 6. Available at: https://www.frontiersin.org/journals/microbiology/articles/10.3389/fmicb.2015.00887/full.

PubMed Abstract | Google Scholar

Malmberg, R. L., Rogers, W. L., Alabady, M. S. (2018). A carnivorous plant genetic map: pitcher/insect-capture QTL on a genetic linkage map of Sarracenia. Life Sci. Alliance 1. doi: 10.26508/lsa.201800146

Crossref Full Text | Google Scholar

McMurdie, P. J., Holmes, S. (2013). phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PloS One 8, e61217. doi: 10.1371/journal.pone.0061217

PubMed Abstract | Crossref Full Text | Google Scholar

Morales Moreira, Z. P., Chen, M. Y., Yanez Ortuno, D. L., Haney, C. H. (2023). Engineering plant microbiomes by integrating eco-evolutionary principles into current strategies. Curr. Opin. Plant Biol. 71, 102316. doi: 10.1016/j.pbi.2022.102316

PubMed Abstract | Crossref Full Text | Google Scholar

Munch, K., Boomsma, W., Willerslev, E., Nielsen, R. (2008). Fast phylogenetic DNA barcoding. Phil Trans. R Soc. B 363, 3997–4002. doi: 10.1098/rstb.2008.0169

Crossref Full Text | Google Scholar

ONT 16S workflow 16S Workflow (2019). EPI2ME WIMP workflow: quantitative, real-time species identification from metagenomic samples (Oxford Nanopore Technologies). Available online at: https://nanoporetech.com/resource-centre/epi2me-wimp-workflow-quantitative-real-time-species-identification-metagenomic (Accessed February 23, 2024).

Google Scholar

Patz, S., Gautam, A., Becker, M., Ruppel, S., Rodríguez-Palenzuela, P., Huson, D. (2021). PLaBAse: A comprehensive web resource for analyzing the plant growth-promoting potential of plant-associated bacteria. 2021.12.13.472471. doi: 10.1101/2021.12.13.472471

Crossref Full Text | Google Scholar

Patz, S., Rauh, M., Gautam, A., Huson, D. H. (2024). mgPGPT: Metagenomic analysis of plant growth-promoting traits. bioRxiv 2021.12.13.472471. doi: 10.1101/2024.02.17.580828

Crossref Full Text | Google Scholar

Saeed, Q., Xiukang, W., Haider, F. U., Kučerik, J., Mumtaz, M. Z., Holatko, J., et al. (2021). Rhizosphere bacteria in plant growth promotion, biocontrol, and bioremediation of contaminated sites: A comprehensive review of effects and mechanisms. Int. J. Mol. Sci. 22, 10529. doi: 10.3390/ijms221910529

PubMed Abstract | Crossref Full Text | Google Scholar

Segata, N., Izard, J., Waldron, L., Gevers, D., Miropolsky, L., Garrett, W. S., et al. (2011). Metagenomic biomarker discovery and explanation. Genome Biol. 12, R60. doi: 10.1186/gb-2011-12-6-r60

PubMed Abstract | Crossref Full Text | Google Scholar

Shannon, C. E. (2001). A mathematical theory of communication. SIGMOBILE Mob Comput. Commun. Rev. 5, 3–55. doi: 10.1145/584091.584093

Crossref Full Text | Google Scholar

Thorogood, C. J., Bauer, U., Hiscock, S. J. (2018). Convergent and divergent evolution in carnivorous pitcher plant traps. New Phytol. 217, 1035–1041. doi: 10.1111/nph.14879

PubMed Abstract | Crossref Full Text | Google Scholar

Trivedi, P., Leach, J. E., Tringe, S. G., Sa, T., Singh, B. K. (2020). Plant-microbiome interactions: from community assembly to plant health. Nat. Rev. Microbiol. 18, 607–621. doi: 10.1038/s41579-020-0412-1

PubMed Abstract | Crossref Full Text | Google Scholar

Turnbaugh, P. J., Hamady, M., Yatsunenko, T., Cantarel, B. L., Duncan, A., Ley, R. E., et al. (2009). A core gut microbiome in obese and lean twins. Nature 457, 480–484. doi: 10.1038/nature07540

PubMed Abstract | Crossref Full Text | Google Scholar

Volvenko, I. V. (2014). General patterns of spatial distribution of the integral characteristics of benthic macrofauna of the northwestern pacific and biological structure of ocean. Open J. Ecol. 4, 196–213. doi: 10.4236/oje.2014.44020

Crossref Full Text | Google Scholar

Wen, T., Xie, P., Yang, S., Niu, G., Liu, X., Ding, Z., et al. (2022). ggClusterNet: An R package for microbiome network analysis and modularity-based multiple network layouts. iMeta 1, e32. doi: 10.1002/imt2.32

PubMed Abstract | Crossref Full Text | Google Scholar

Wheeler, G. L., Carstens, B. C. (2018). Evaluating the adaptive evolutionary convergence of carnivorous plant taxa through functional genomics. PeerJ 6, e4322. doi: 10.7717/peerj.4322

PubMed Abstract | Crossref Full Text | Google Scholar

Whitham, T. G., Bailey, J. K., Schweitzer, J. A., Shuster, S. M., Bangert, R. K., LeRoy, C. J., et al. (2006). A framework for community and ecosystem genetics: from genes to ecosystems. Nat. Rev. Genet. 7, 510. doi: 10.1038/nrg1877

PubMed Abstract | Crossref Full Text | Google Scholar

Whitham, T. G., DiFazio, S. P., Schweitzer, J. A., Shuster, S. M., Allan, G. J., Bailey, J. K., et al. (2008). Extending genomics to natural communities and ecosystems. Science 320, 492–495. doi: 10.1126/science.1153918

PubMed Abstract | Crossref Full Text | Google Scholar

Wilcoxon, F. (1992).Individual Comparisons by Ranking Methods. In: Breakthroughs in Statistics (New York, NY: Springer New York). Available online at: http://link.springer.com/10.1007/978-1-4612-4380-9_16 (Accessed December 18, 2023).

Google Scholar

Yan, L. (2021). ggvenn: Draw Venn Diagram by ‘ggplot2.’ R Package Version 19.

Google Scholar