2.1 Coral dataset and experiment

The experimental design and detailed results of the coral dataset we performed LDA on can be found in [19]. The following few sentences summarize the key experimental features and study goals by Rosales et al. [19], who identified bacterial taxa linked to resistance and susceptibility to White Band Disease (WBD) in two species of Acropora coral. In their experiment, diseased tissue was grafted to healthy Acropora cervicornis and Acropora palmata to induce a diseased state. Bacterial communities were sampled prior to grafting (control) and seven days after inoculation (treated). Seven days after inoculation, corals were also examined for visual signs of disease lesions and classified based upon their disease status (diseased or visually unaffected) and severity (low, mid, or high disease susceptibility). Microbial community composition at the ASV level was then compared across treatment and outcome groups using “traditional” community-based methods, including ordinations, differential abundance of taxa, and core microbiome analysis, to identify taxa that 1) may cause WBD, 2) are associated with microbiome changes during disease, and 3) potentially increase Acropora resistance to WBD. Rosales et al. [19] used pairwise comparisons with a Kruskal-Wallis test to assess the significance of Shannon diversity and evenness in host species. ANOVA and PERMANOVA were used to test the significance of dispersion of the samples and to find significant interactions between groups, respectively. The authors also present beta-diversity (PCA with Euclidean distance were used for visualization) and microbial differential abundance plots to identify the difference between treatments.

2.2 Maize dataset and experiment

The maize dataset for this work comes from a serial microbiome propagation studying on how microbiome composition in soil can affect performance of maize and which plant traits are most affected by the microbiome. In this experiment we were interested in understanding the behavior of several plant traits such as stomatal closure point, water use efficiency, maximum rate of photosynthesis, rate of stomatal conductance, drought time (time to permanent stomatal closure from the beginning of the drought treatment), percentage of leaf water content, leaf mass per area, stem diameter, stem height, and root biomass to see if microbiomes that support these traits will improve plant performance under drought and help to find stable microbiomes that help plants to withstand water-limited conditions. While all plant traits are relevant for assessing plant health under drought conditions, we are mainly interested in behavior of stomatal closure point and water use efficiency. Stomatal closure point (SCP) is an important plant desiccation tolerance metric that is linked to plant hydraulics and productivity under drought [26–29]. Water use efficiency (WUEi), on the other hand, is a traditional metric of plant productivity and drought resistance, and it measures the carbon gain of a plant per water lost [30]. Based on the known influence of rhizosphere microbiomes on plant nitrogen availability [31], and the connection between improved nitrogen availability and higher water use efficiency [32, 33], it is plausible to expect that WUEi and SCP could be influenced by the rhizosphere microbiome.

In this work, we present results from data on 119 plant microbiomes from two generations (called generation 0 and 1) of serial microbiome propagation (Fig 1). These analyses on two generations were performed in order to understand the effect of different water treatments on plant-microbiome interactions with the system. In each generation, 64 maize plants were grown in a greenhouse setting between January and May 2020 in individual 2.6 gal (9.8 L) pots filled with 6L of fritted clay (GreensGrade, 20–50 mesh size: Profile (Buffalo Grove, IL, USA)). The plant seed for each pot in each generation was randomly drawn from a stock of an experimental strain of maize, USDA seed bank inbred line “B73”. For generation 0, 48 of the seeds were inoculated using methods described by [34] by microbial communities originating from one of two natural soils: one collected from an agricultural field near Akron, CO, USA and the other collected from a ponderosa pine forest near Los Alamos, NM, USA, and 16 seeds were planted without inoculation to form a control treatment. After seed germination, the plants were divided into two water treatments: full-water (up to 65% volumetric water content 3 times a week) and half-water (up to 45% volumetric water content 3 times a week). After the plants had grown to a stage showing 9 fully grown leaves (∼ 8 weeks), the watering was stopped and the plants were allowed to dry under complete water withdrawal until stomatal closure occurred (terminal drought). Before the terminal drought, the soil in each pot was sampled for microbiome and soil chemistry analyses, plants were sampled for leaf chemistry analysis, and measured for height, stem diameter, maximum photosynthesis rate, and stomatal conductance (infrared gas analyser; Licor 6400, Licor Inc. (Lincoln, NE, USA)) from which water use efficiency was calculated by dividing the maximum photosynthesis rate by the measured stomatal conductance. During terminal drought, stomatal conductance of each plant was measured daily. When stomatal conductance reached zero, a leaf was cut for leaf water potential measurements (pressure chamber method, PMS Instruments (Al bany, OR, USA)) to determine the stomatal closure point (leaf water potential at which the plant closes stomata), and the time since the beginning of the terminal drought was recorded (drought time). For generation 1, new pots with fritted clay were set up, and new seeds were inoculated by serially transferring the microbial communities from the generation 0 pots with no selection with similar microbiome extraction as for generation 0. Therefore, plants in generation 1 were growing with the microbiomes that were direct descendants of those from generation 0. Due to the high density of roots in the pot, all of the soil was considered “rhizosphere”, thus we chose to use soil cores to sample the rhizosphere community in the pots, rather than performing a separate analysis on soil directly adhered to roots, and because they were a more appropriate reflection of the communities we transferred to the next generation. The same water treatments as in generation 0 were imposed after seed germination, but now so that for half of the plants in each treatment the water treatment was switched. This created two additional treatment categories (stable watering vs. switched watering) for microbiomes from each original soil source and the non-inoculated controls. Similarly to generation 0, generation 1 plants were grown to a stage showing 9 fully grown leaves (∼ 8 weeks), and the plant performance measurements and the terminal drought treatment were conducted as for generation 0. For more detailed information on the experimental setup, plant growing conditions, and plant trait measurements see S1 Text.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Experimental design for maize.

Two generations (G₀, G₁) of maize experiment were grown in the greenhouse with microbiomes extracted from agricultural (M_Ag) and forest (M_Forest) soils. Generation 1 plants were grown with microbiomes that were direct descendants of those from generation 0. full-water (W_w) and half-water (L_w) treatments were imposed to plants during generation 0. The same water treatments were imposed in generation 1 after seed germination, but now so that for half of the plants in each treatment the water treatment was switched.

https://doi.org/10.1371/journal.pcbi.1011075.g001

Of our original 128 microbiome samples, 9 samples were omitted from the final analysis due to insufficient DNA, sequencing failures, or low read quality and quantity, resulting in the analysis of communities from only 119 plants.

In this paper, we analyze microbiome composition with a focus on the links between microbiomes and watering conditions, initial microbiome soil source, and generation of the experiment.

2.3 DNA extraction and microbiome sequencing

Microbial DNA was extracted from homogenized soil core samples using the Qiagen DNeasy PowerSoil kit, with modifications to improve yield from clay soils [35]. 16s rRNA gene amplicon sequencing (V4 region) was used to profile bacterial communities. DNA amplicons were generated using 515F-R806 previously described primers [34] and then sequenced using an Illumina MiSeq (300bp, paired end reads). Raw sequencing reads were demultiplexed using USEARCH11 [36] and dada2 was used to perform quality filtering, primer removal, and read denoising [37]. Default settings were used for dada2 according to the dada2 tutorial (v1.8), except the filterAndTrim function used the following settings to remove primers: truncLen = c(240, 200), truncQ = 2, trimLeft = c(25, 26), maxEE = c(2, 4); only paired reads with a minimum overlap of 100bp were merged, and merged sequences between 250 and 260bp were used to generate the sequence abundance table. Taxonomic classifications were assigned to each unique sequence at the 80% confidence level using the SILVA v138.1 database [38, 39]. Unprocessed sequence data have been deposited to the NCBI Sequence Read Archive under the projects PRJNA780613 (generation 0 samples) and PRJNA 780954 (generation 1 samples).

2.4 Topic modeling

LDA, commonly used in text analysis, is fully described in [4]. Sankaran et al. [5] developed some guidelines of the application of probabilistic latent variable models including LDA to human microbiome data. Before we start to discuss our LDA results, let us introduce the analogy between text and soil microbiome analysis we used: the pot samples, bacterial species (taxa), microbial communities are viewed as the documents, words, and topics, respectively. Therefore, the bacterial abundance counts matrix is viewed as the document-term matrix. That is to say that individual pot microbiomes (documents) are broken down into a distribution of topics, and these topics are distributions of taxa (words). A topic is described as a microbial community that may share similar biological functions.

Fig 2 shows a flow chart of LDA topic modeling at phylum level for K = 3 topics. After soil microbiomes were sequenced, microbiome composition was aggregated to a certain taxonomic level, in this case to the phylum level. The desired number of topics K (we chose the optimal number as described in Results section) was specified in advance and the bacterial counts matrix was used as an input for the LDA algorithm. LDA infers taxa distribution for each topic and topic distribution for each sample. In this analysis, we used the LDA method implemented in the MALLET software [40] to identify topics at different taxonomic levels. To learn the model, we used Gibbs sampling [41]. We ran LDA for 10,000 iterations to allow log-likelihood per word to stabilize.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Flowchart of the study design.

Schematic LDA topic modeling for K = 3 topics.

https://doi.org/10.1371/journal.pcbi.1011075.g002

2.5 Data analysis

Our data is high-dimensional and sparse with a large amount of ASV sequences for which no complete taxonomic classification can be assigned. Due to the differences between text data and microbiome data, some pre-processing approaches and evaluation metrics, such as the removal of common words (stop words) which are assumed not to be useful in describing the document, used for text data were not appropriate in our case. We pre-processed the data by aggregating ASVs to identified taxonomic descriptions associated with the ASV (phylum, class, family, order, and ASV (no aggregation here)), and handled incomplete taxonomic specification for some ASVs by aggregating to the highest known taxonomic level.

Many bacterial taxa present in soils remain unidentified at some taxonomic levels. Despite advances in sequencing technology that allow deep profiling of the taxonomic composition and functional potential of soil microbial communities, databases containing information about environmental microbes remain sparse. Nevertheless, unidentified taxa, also known as “microbial dark matter” are sometimes in high abundance, and can be essential to microbial community interaction networks where they are often found to be major hubs or keystone taxa in these networks, suggesting that they provide important functions and stability to the community [42, 43]. We elected to include unidentified taxa in our analyses, and aggregated these by the most specific known taxonomic classification. Analysis at the ASV taxonomic level (using ASVs directly) avoids the issue of incomplete taxonomic classification. Using this approach, we could perform the analysis on the complete community, rather than only on the classified subset, which may have changed community network structure [43], and thus the outcomes of the LDA analysis.

We omitted the analysis at genus and species levels because more than half of the taxa were unidentified and thus aggregation was harder to perform at these levels. However, we performed the analysis at the ASV (“deepest”) level which allowed us to determine which ASVs (that belong to certain genera/species) were probable in the LDA topics. When analyzing at certain taxonomic level (except at ASV level), we removed taxa with aggregated count of less than 4. We then performed the same analysis described below for each taxonomic level.

The main parameter of LDA that adjusts the level of detail in the resulting model is the number of topics. If one chooses far too many topics, the main drawback is that the factorization becomes less well-conditioned, resulting in redundant topics and an ambiguous document-topic matrix. With too few topics, not all patterns in the data are captured. We chose the number of topics carefully based on topic selection metrics and also checked that using a different number of topics does not reduce the number of associations significantly. We ran 5-fold cross-validation, measured several metrics, and calculated averaged results as a function of the number of topics. We repeated all calculations five times to account for the method’s stochastic nature. We examined the metrics of perplexity, pairwise cosine similarity, coherence, and exclusivity [40] to help us decide how many topics to use at each level of taxonomic identification (S1 Fig). The perplexity score was calculated to see how well a model performs on the unseen held-out test data (20%). A lower perplexity score indicates better generalization performance. A pairwise cosine similarity was calculated between taxa in topics to determine how distinct were the distributions of taxa between topics. We used coherence [44] to measure whether the most probable taxa in a topic tended to co-occur in other topics. The exclusivity metric [40, 45] we report, finds the number of most unique taxa for each topic.

After selecting the number of topics, we ran LDA once again on the full dataset. We then discussed the weighting of topic abundances in terms of experimental conditions, such as water treatment, watering stability treatment, soil source microbiome inoculation type, and generation. We also calculated the number of topics that had strong or moderate association with experimental conditions at each taxonomic level. A topic was moderately associated with the experimental condition (e.g., half-water treatment) if average topic abundance in all pots treated under this condition (e.g., half-water treatment) was at least twice as large (67–80% abundance) but less than four times larger than the average topic abundance under alternative condition (e.g., full-water treatment). A topic was strongly associated with the experimental condition if average topic abundance in all samples treated under this condition was at least four times larger (80–100% abundance) than the average topic abundance under alternative condition. If there were two alternative conditions (e.g., three soil source microbiome inoculation types) then the strength of topic abundance was determined by comparing it to topic abundances under both alternative conditions. We present topic abundance weighting plots that show topics related to different treatments. The weighting was determined by how much the average topic abundance under one treatment was smaller than that under the alternative treatment on a 0–1 scale. LDA outputs probability distribution over topics for each document (sample). To determine topic weighting in case of two treatments (e.g., half-water vs. full-water) for each topic i and treatment t we first calculated , where r_j is a contribution of the topic in document j and n_t is the number of samples under a certain treatment (e.g., half-water). For each topic i we defined , where K is the number of topics in the model. We then obtained topic weighting w_i,t for each treatment t by normalizing x_i,t values across two treatments so that the sum of the contributions for all topics for each treatment is one. To obtain topic abundance weighting plots we calculated plot coordinates where treatment t₁ is located at the bottom or left half of the plot (e.g., half-water, generation 0) and treatment t₂ is located at the top or right half of the plot (e.g., full-water, generation 1) as follows. If , then topic coordinate on the plot is , otherwise if , then topic coordinate on the plot is .

Each LDA microbial community topic is characterized by the distribution of taxa. This allowed us to detect which topics were more associated with pots that were exposed to water-limited treatment or pots with inoculated microbiomes from different soil sources, etc. To determine which taxa abundances were most amplified (exclusive) in the topic in comparison to their abundances in the overall dataset, we defined the term relative amplification (also known as lift in text analysis [46]) of each taxon within each topic. It is defined as the probability of the taxon given the topic divided by the frequency of that taxon in the overall dataset, and then we normalized each topic across all taxa so that the sum of the relative amplifications for all taxa in each topic is one. In other words, the relative amplification shows which taxon’s abundance is most amplified in the topic in comparison to its abundance in the overall dataset. We also calculated effective number of words for each topic [40], which is computed for each topic as the inverse of the sum of the squared probability of each word in the topic. This metric counts how many taxa are important contributions to each topic, that is the number of taxa that the topic is effectively spread across.

We tested if topics had statistically significant relationships with plant traits. To do this, we calculated a Spearman’s rank correlation coefficient between each topic and each continuous plant trait. Tests were performed at 5% statistical significance level with Holm–Bonferroni correction, which is uniformly more powerful than the Bonferroni correction, to control the family-wise error rate.

To compare the LDA approach with traditional microbiome analysis, we ran Spearman’s rank correlation tests and Indicator Species analysis [47–49] to find significant relationships between individual taxa and plant traits. Indicator Species analysis (using the multipatt function in the IndicSpecies v1.7.9 R package) did not identify any taxa significantly associated with any treatment following multiple test correction. Differential abundance analysis was conducted using Corncob R package [50]. It identified numerous taxa that significantly differed in abundance between treatments. Corncob uses a beta-binomial model to compare taxa abundances using un-normalized sequencing data.

To test the robustness of our LDA results that might vary due to the randomness of LDA, we compared results from three different LDA runs at the phylum level. Cosine similarity between taxa in topics obtained from the different runs was calculated to determine the consistency of the topics from run to run. S5 Table shows the consistency of the topics between three different runs. Although topics related to the watering treatments have lowest cosine similarity between runs, the most probable taxa in these topics were the same.

We also fitted different number of topics at the phylum level to see if this significantly affects the results (S6 Table). Even if using more topics can result in more associations between topics and treatments, cosine similarity between topics within a model can be used to avoid overfitting.

3.1 LDA results on coral dataset

To assess the performance of LDA-based microbiome analysis prior to analyzing our own maize dataset, we performed LDA analysis on a previously published dataset [19] that used “traditional” microbiome analysis methods.

By running LDA with 20 topics, we found 7 and 11 topics that were strongly associated with Acropora cervicornis and Acropora palmata species, respectively. Most of the topics were linked to some experimental conditions such as treatment (control and inoculated), outcome (control, visually unaffected, and diseased), and disease susceptibility (low, mid, and high) (Figs 3 and 4). We found that Midichloriaceae was the dominant family in topics relevant to A. cervicornis species, and the families Spirochaetaceae and Endozoicomonadaceae were found to be dominant in topics abundant in samples from A. palmata species (Fig 5).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Topic abundance weighting for the treatment and species types.

All topics to the left of topic 9 are associated with Acropora cervicornis species. All topics to the right of topic 8 are associated with Acropora palmata species. All topics below topic 11 are associated with Control treatment whereas all topics above topic 15 are associated with Inoculated (diseased or visually unaffected) treatment. Topics in color are topics associated with an experimental outcome. Topic 20 is associated with Low disease susceptibility; topics 5 and 14 are associated with Medium disease susceptibility; topics marked with triangle up symbols are associated with High disease susceptibility.

https://doi.org/10.1371/journal.pcbi.1011075.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Ternary plot of topic abundances for the outcome type (control, diseased, visually unaffected).

https://doi.org/10.1371/journal.pcbi.1011075.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Distribution of ASV sequences in each learned LDA topic.

Topics to the left of topic 9 are strongly associated with Acropora cervicornis species, and topics to the right of topic 8 are strongly associated with Acropora palmata species. The labels are written in the phylum_class_order_family format, the names of genera, species, and ASVs are not shown. Only probabilities (color circles) greater than 0.1 are shown. Smaller circles are displayed on the top of larger circles.

https://doi.org/10.1371/journal.pcbi.1011075.g005

Two topics 4 and 16, each representing different Acropora species, were connected to control samples. Topic 4, linked to Acropora cervicornis, was dominated by the presence of ASVs from the families Midichloriaceae, P3OB-42, Endozoicomonadaceae, Nannocystaceae, and several ASVs from to the order Phormidesmiales. On the other hand, topic 16, representing control samples from A. Palmata species, was dominated by the presence of ASVs from Rhodobacteraceae, Colwelliaceae, and Phormidesmiaceae families.

We found three topics related to visually unaffected samples, ASVs from the families Midichloriaceae, Spirochaetaceae, Francisellaceae, Kiloniellaceae, and Endozoicomonadaceae in A. cervicornis were dominant in topic 10, and ASVs from Halobacteroidaceae, Rhodobacteraceae, Spirochaetaceae, Desulfobacteraceae, Clostridiaceae, Lachnospiraceae, and Vibrionaceae families were more probable in topic 2 related to visually unaffected A. palmata species. Topic 20 was the only topic associated with A. cervicornis species and low disease susceptibility. It was mainly dominated by the ASVs from Midichloriaceae and Rhodobacteraceae families and ASVs from the orders Nostocales and Alteromonadales. These results suggest that these taxa may promote resistance to disease in Acropora species. On the contrary, topic 19, associated with A. cervicornis species and high disease susceptibility, was dominated by the presence of ASVs mainly from Endozoicomonadaceae, Midichloriaceae, and P3OB-42 families. Several topics were relevant to A. palmata species and high disease susceptibility and were dominated by the presence of ASVs from Fulvivirga, Arcobacter, Oscillatoria genera, Nostocales order, and Alphaproteobacteria class.

The LDA results with the most probable taxa in each topic associated with one of two species, certain treatment or outcomes mostly agreed (Table 1) with the findings in [19]. LDA results agreed with the authors findings of the dominant families related to both species, control samples in Acropora cervicornis species, and low and high disease susceptibilities in Acropora cervicornis species. Rosales et al. [19] found a core bacteria member in both species from the order Myxococcales (P3OB-42) at relatively higher abundances in corals with lower rates of disease development following grafting. Though P3OB-42 was not found among most probable families in topics 2 and 10, it was the most probable family in topic 8 which was related to visually unaffected samples without being ascribed to either species. A single ASV, designated Sphingobium yanoikuyae (family Sphingomonadacea), was significantly detected in both species of disease-exposed and visually unaffected samples [19]. LDA results in agreement with this finding, Sphingobium yanoikuyae was found in topics 17, 19 (diseased) and 8 (visually unaffected). Rosales et al. [19] found Ralstonia genus at higher relative abundances in negative control samples. This genus was the most dominant taxon (22.4%) in LDA topic 9 associated with only control samples and nothing else.

There were some taxa-treatment associations identified by Rosales et al. [19] that LDA failed to identify. In particular, LDA did not detect Spirochaetaceae and Endozoicomonadacea as most dominant families in control samples of A. Palmata species. Rosales et al. [19] found relatively higher abundances of the genus Vibrio in visually unaffected and diseased corals compared to control corals in A. cervicornis. LDA topic associated with visually unaffected corals did not show higher abundances of this genus compared to that in the topic associated with control samples. However, genus Vibrio was the most dominant taxon in topic 5 associated with medium disease susceptibility. For outcome in A. palmata, the families Rhodobacteraceae (genus HIMB11) and Cryomorphaceae (uncultured genus) were significantly abundant and highly associated with visually unaffected and diseased corals [19]. LDA found higher associations of these genera in topic 6 associated with high disease susceptibility. Therefore, the biggest disagreement between LDA results and relative/differential abundances analyses results from [19] was related to diseased samples in both coral species and to control samples in A. palmata species (Table 1). It should be noted that the standard deviations in relative abundance analysis were very high which made the comparison of average relative abundances in samples with LDA probabilities harder.

Though LDA did not detect all associations present in [19], LDA results revealed additional taxa that were linked to certain treatments, outcomes and disease susceptibility that were not previously discussed. Specifically, Rosales et al. [19] did not find any strong associations of taxa and high disease susceptibility in A. palmata species samples, while LDA did.

This suggests that LDA was capable not only to detect most associations as were found by the researchers in the previous study but also complement the results with taxa that were less abundant in data but still may interact with other taxa towards the White Band disease resistance.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Comparison of LDA results and results obtained from the relative abundance analysis (Tables 2 and 3 in [19]).

Table shows mean and standard deviation (SD) relative abundances (RA) in % by taxon of core microbiomes per experimental outcome and coral-host [19] as wells as probability % of these taxa in LDA topics associated with corresponding outcome and coral species. Multiple LDA topics were associated with Diseased outcome in both coral species. Topic 20 and topic 5 were associated with low and medium disease susceptibility, respectively, whereas topics 19, 1, 6, and 18 were associated with high disease susceptibility.

https://doi.org/10.1371/journal.pcbi.1011075.t001

3.2 LDA results from our maize experiment

LDA results provide a window into which taxa might function together as a system and be related to different experimental conditions and plant traits. In this section, we present the analysis results at phylum, class, order, family, and sequence (ASV) taxonomic levels from the microbiome sequential propagation experiment with maize. Based on the topic selection metrics mentioned above, we decided to use K = 6, 8, 10, 20, and 25 topics at the phylum, class, order, family, and ASV levels, respectively, as we did not observe significant improvement in metrics behaviors for larger topic counts. At the phylum level we observed that perplexity dropped at K = 6, cosine similarity was the same for K < 7 (S1(a) Fig), and coherence was not changing much for K = 5, 6, 7 (S1(b) Fig). Therefore, we decided to use K = 6 topics at the phylum level. Detailed plots and tables for each taxonomic level are available in S2–S14 Figs and S1–S19 Tables, including explicit lists of all topic compositions and their associations with treatment conditions and interactions with plant traits. We present summary tables (Tables 2, 3, 4 and 5) that capture the main results of the analyses at all taxonomic levels.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Number of topics associated with experimental conditions at each taxonomic level.

Strong (topic abundance is four times (×4) larger than for other treatments) and moderate (topic abundance is twice (×2) as large but less than four times larger than that for other treatments) topic associations are shown for each experimental condition. The numbers in parentheses represent number of topics associated with at least one experimental condition and overall number of topics used in LDA at taxonomic level.

https://doi.org/10.1371/journal.pcbi.1011075.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Number of topics at each taxonomic level that showed statistically significant relationships between topics and plant traits based on a Spearman’s rank correlation coefficient with Holm–Bonferroni correction.

The numbers in parentheses represent number of different topics associated with significant relationships and overall number of topics used in LDA at taxonomic level. The last four columns are Leaf mass per area, % Leaf water content, Stomatal conductance, and Water use efficiency.

https://doi.org/10.1371/journal.pcbi.1011075.t003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Most probable and amplified taxa in some topics associated with experimental conditions at phylum, class, and order taxonomic levels.

Left: average percentages of most abundant taxa in pots for each taxonomic level. Middle: topic association with experimental conditions (HW: half-water, FW: full-water, G0: generation 0, G1: generation 1, AG: agricultural soil source). Right: most probable taxa in topic with probabilities are displayed in each upper row; most amplified taxa with its value (see Methods) are displayed in each lower row, some taxa were unidentified at lower levels. Note that unidentified taxa were aggregated to the highest known taxonomic level, meaning that they may represent one or several classes, orders, etc. Not all topics are listed.

https://doi.org/10.1371/journal.pcbi.1011075.t004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Most probable and amplified taxa in topics strongly associated with experimental conditions at family and ASV taxonomic levels.

Left: average percentages of most abundant taxa in pots for each taxonomic level. Middle: topic association with experimental conditions (HW: half-water, FW: full-water, G0: generation 0, G1: generation 1, AG: agricultural soil source, FR: forest soil source, ST: stable watering, SW: switched watering). Right: most probable taxa in topic with probabilities are displayed in each upper row; most amplified taxa with its value (see Methods) are displayed in each lower row, some taxa were unidentified at lower levels. Note that unidentified taxa were aggregated to the highest known taxonomic level, meaning that they may represent one or several classes, orders, etc. Not all topics are listed.

https://doi.org/10.1371/journal.pcbi.1011075.t005

When analyzed at different taxonomic levels, LDA allowed us to identify topics moderately and strongly related to the experimental conditions at each level (Table 2). We found that microbiome topics were more associated with experimental conditions at lower taxonomic levels than at higher taxonomic levels as expected based on high functional diversity that can be present in groups belonging to the same high taxonomic level (Table 2). Topics were associated mostly with water treatment and plant generation, and less with soil microbiome inoculation source type and stability of the watering treatment. We found that at each taxonomic level there was at least one topic that was associated with the full- or half-water treatments, however, strong associations were only observed starting from the order taxonomic level. Strong topics associations with the generation and particular soil microbiome inoculation source type were detected only at order level and lower. We observed significant skewness of topics toward the stability of watering treatment only at the lowest ASV taxonomic level. At the order level and lower, several topics that were simultaneously strongly associated with multiple experimental conditions appeared. For example, one of the topics was strongly associated with half-water treatment and was more prevalent in generation 1 plants with agricultural microbiome inoculation. At both family and ASV taxonomic levels, all topics were moderately or strongly related to at least one experimental condition. We also detected many moderate associations of topics with experimental conditions at each level.

We found statistically significant relationships between learned topics and plant structural traits, including stem diameter, height, and root biomass at each taxonomic level (Table 3). We did not observe any statistically significant relationships between learned topics and functional plant traits (stomatal closure point (SCP) and water use efficiency (WUEi)) except at the ASV taxonomic level where we found a significant relationship of one topic with WUEi. Overall at each taxonomic level, some relationships were significant (p-value (SCP) < 0.001 and p-value (WUEi) < 0.05), but the corrections for multiple testing precluded their significance.

3.3 Identifying taxa related to experimental conditions

We present taxa along with their probabilities and relative amplifications in several selected topics strongly and moderately associated with at least one experimental condition in Tables 4 and 5. For each topic we calculated the effective number of taxa, and we highlighted the taxa whose ranks are under this threshold in bold (S1, S7 and S10 Tables), indicating which taxa are the most important contributions to the content of this topic. For phylum-level analysis, this was the top 2–4 taxa for each phylum in the topic. For the ASV level, each topic was spread across several dozen taxa.

According to the previous studies conducted on different host plants across different soils, Actinobacteriota and Firmicutes tend to increase in response to drought, Bacteroidota and Verrucomicrobiota decrease, and Proteobacteria, Planctomycetota, and Acidobacteriota compositions change [21, 25]. LDA results at phylum taxonomic level mostly agreed with these trends. In particular, Proteobacteria and Verrucomicrobiota were dominant in the topic linked to full-water treatment, whereas Actinobacteriota was dominant phyla in the topic related to half-water treatment. Bacteroidota was dominant in both half and full-water associated topics. Firmicutes phylum was present only in 1% of data, as expected it did not showed up among most probable phyla in any topic. However, Firmicutes was the one of the most amplified phyla in the topic related to half-water treatment (Table 4).

Microbial communities detected in class-level topics related to water treatment followed the same trends as at the phylum level. In particular, Actinobacteria was the most probable class at half-water associated topics, and Verrucomicrobiae was the most dominant class at full-water associated topics. Classes Alphaproteobacteria and Bacteroidia were dominant in both half- and full-water associated topics. Bacilli in Firmicutes phylum was the most amplified class in the topic related to half-water treatment.

Starting from order taxonomic level, most topics have become associated with multiple experimental conditions, and therefore it became more difficult to highlight orders that were only associated with one of the conditions. We note that several most amplified orders were unidentified ones. This result highlights an importance of keeping unidentified taxa in the analysis. Orders Micrococcales in Actinobacteria class, Sphingobacteriales in Bacteroidia class, and Rhizobiales in Alphaproteobacteria class formed a microbial community associated with water-limited soil. On the other hand, orders Burkholderiales (Gammaproteobacteria) and Rhizobiales (Alphaproteobacteria) were dominant in the topic linked to full-water treatment. Another Verrucomicrobiae order, Verrucomicrobiales, was the most probable in the topic linked to full-water treatment and to plants without soil microbiome inoculation.

Analysis at family level showed that family Micrococcaceae (Actinobacteria) together with a few families from Sphingobacteriales order formed a community linked to half-water treatment. Topics mainly consisted of families Verrucomicrobiaceae (Verrucomicrobiae) and Comamonadaceae (Gammaproteobacteria) were linked to full-water treatment.

All topics at ASV taxonomic level were moderately or strongly related to at least one experimental condition (Fig 6). An ASV from Pseudarthrobacter genus (from 5 to 18% in LDA topics) along with an ASV from Methylobacterium-Methylorubrum genus (Alphaproteobacteria class) (3.5–9%) were found in topics linked to half-water treatment (Fig 7). Several ASVs from Verrucomicrobiales order (from 5 to 18% in LDA topics) along with an ASV in Comamonadaceae family (Gammaproteobacteria class) (3.5–11%) formed a community linked to full-water treatment. The microbial community consisted mainly from three ASVs in Streptomyces genus (18.2%), autotrophicum/hiltneri/lusitanum species (19.2%), and Opitutaceae family (16.2%), respectively, was associated with half-water treatment and generation 1 plants with agriculture soil microbiome inoculation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Topic abundance weighting for the water treatment and generation on a 0–1 scale.

All topics to the left of topic 5 are associated with generation 0. All topics to the right of topic 25 are associated with generation 1. All topics below topic 25 are associated with half-water treatment whereas all topics above topic 15 are associated with full-water treatment. Topics are colored according to the soil microbiome inoculation source. The shapes of the markers corresponding to the topic association with stability of watering treatment. Note: the weighting was calculated by how much average topic abundance under one treatment was different from that under alternative one (see Methods).

https://doi.org/10.1371/journal.pcbi.1011075.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Distribution of ASVs in each learned LDA topic.

Topics to the left of topic 10 are associated with half-water treatment, and topics to the right of topic 15 are associated with full-water treatment. The labels are written in the phylum_class_order_family_genus_species format, the names of genera, species, and ASV are not shown. Only probabilities (color circles) greater than 0.025 are shown. The sizes of the circles representing probabilities are multiplied by 4 for visualization purposes. Smaller circles are displayed on the top of larger circles.

https://doi.org/10.1371/journal.pcbi.1011075.g007

The analyses on different taxonomic levels gave insights into the homogeneity of the function of the lower level microbial communities that compose higher level taxonomic communities. Certain taxa were highly probable in topics associated with the same treatments at each taxonomic level. This happened either because a certain ASV was highly abundant in the data and was linked to a certain treatment which caused the same associations at higher taxonomic levels, or because multiple ASVs had the same treatment associations at higher taxonomic levels.

In some cases, single ASVs were associated with the same experimental conditions at different taxonomic levels. An ASV belonging to Ralstonia genus (Gammaproteobacteria class) was highly probable in one single topic that did not show any associations with treatments at phylum, class, and order levels. However, at family and ASV levels this ASV was the most probable (45.94% at family and 49% at ASV levels) in the topic linked to half-water treatment and generation 1 plants with forest soil microbiome inoculation. An ASV from autotrophicum/hiltneri/lusitanum species (Gammaproteobacteria class) was detected in topics related to half-water treatment at ASV (19.2%), family (8.6%) and order (16.3%) levels. An ASV from Rhodobacter genus (Alphaproteobacteria class) was prevalent in topics related to generation 0 and full-water treatment at ASV (7.1%), family (19.5%) and order (8.4%) levels. An ASV from Vampirovibrio genus (Cyanobacteria phylum) was detected in topics related to plants with agricultural soil microbiome inoculation at ASV (5.3 and 7.1%), family (37.5%), order (7.01%), and class (18.7%) levels. An ASV from Pseudarthrobacter genus (Actinobacteria class) was highly probable in topics associated with half-water treatment. There were several other ASVs from Micrococcales order that contributed to half-water associated topics but their probabilities were much smaller than that of Pseudarthrobacter.

Multiple ASVs appeared in the same topics associated with certain experimental conditions. Two ASVs from Deinococcus genus were highly probable in only one topic associated with generation 0 at each taxonomic level. Three ASVs from Streptomyces genus (Actinobacteria class) were abundant in the topic related to stable half-watering treatment and generation 1 plants with agricultural soil microbiome inoculation at order taxonomic level and lower. Four ASVs from Sphingobacteriales order (Bacteroidia class) were highly abundant in topics linked to half-water treatment and plants without soil microbiome inoculation at family and ASV taxonomic levels. Starting from order level, five ASVs in Chitinophagacea family (Bacteroidia class) were found in topics associated with half-water treatment, whereas two ASVs in Dyadobacter genus (Bacteroidia class) were found in topics associated with full-water treatment. At family and ASV taxonomic levels, two ASVs in Microscillaceae family (Bacteroidia class) were linked to topics related to plants with agricultural soil microbiome inoculation.

3.4 Correlation and differential abundance analysis

We saw significant Spearman’s correlations between learned topics and some plant morphological traits at each taxonomic level (Table 3). However, we did not find any statistically significant correlations between topics and plant functional traits, such as SCP and WUEi, except for one topic at ASV taxonomic level that was correlated with WUEi.

To compare the LDA approach with traditional microbiome analysis, we ran Spearman’s rank correlation tests between taxa and plant traits at both phylum and ASV levels. We used an indicator species analysis to find significant relationships between individual taxa and plant traits. We also conducted differential abundance analysis to identify taxa that significantly differed in abundance between treatments. Indicator species analysis did not identify any taxa significantly associated with any treatment following multiple test correction. Both correlation analyses done for individual phyla and on LDA topics mostly agreed with each other. Out of 3946 ASVs, 108 of them were found to have significant correlations with root biomass, drought time, stem height, stem diameter, and leaf mass per area (S19 and S20 Tables). At the phylum level, out of 27 phyla, 11 were found to have significant correlations with root biomass, drought time, stem height, and stem diameter (S3 Table). We expect to see taller plants, plants with increased leaf mass per area, and plants with the smaller values of the drought time and stem diameter for the plants grown under half-water treatment than under full-water treatment. Therefore, taxa abundant in water-limited soil were expected to be positively correlated with the stem height and negatively correlated with the drought time, stem diameter, and root biomass, and show the opposite behavior under the full-water treatment. On average, plants in generation 0 were taller, had thicker stems and larger root biomasses than generation 1 plants. Correlation analysis at the phylum level showed that Actinobacteriota was negatively correlated with the drought time. Both Crenarchaeota and Bdellovibrionota were negatively correlated with the stem diameter, Crenarchaeota was also negatively correlated with the root biomass, and Bdellovibrionota was positively correlated with the stem height. Thermoplasmatota was positively correlated with the stem diameter and root biomass and negatively correlated with the stem height. Gemmatimonadota and Verrucomicrobiota were positively correlated with the stem diameter and the drought time. WPS-2 was positively correlated with the stem diameter and Deinococcota was positively correlated with the root biomass. These results obtained from the correlation analysis performed on individual taxa were in agreement with the results from LDA analysis. Actinobacteriota, Bdellovibrionota, Crenarchaeota amplified in half-water associated topic were negatively correlated with the stem diameter and positively correlated with the leaf mass per area. Gemmatimonadota, Verrucomicrobiota, and WPS-2 were amplified only in the full-water treatment associated topic and this topic was negatively correlated with the leaf mass per area. Deinococcota was found probable in generation 0 associated topic which was positively correlated with root biomass and stem diameter, and negatively correlated with the stem height. LDA topic results mostly agreed with the results obtained from the correlation analysis performed on individual taxa and agreed with the bacteria abundance trends observed in other studies. Instead of analyzing individual taxa, LDA learned topics (groups of microbial communities) that may be associated with experimental conditions and correlated with plant traits. LDA analysis directly highlights communities in the microbiome that can be generically expected to be related to the conditions of the plant and the response of the plant to those conditions.

Differential abundance analysis showed plenty of significant ASVs that differ in abundance between treatments (Figs 8 and S12, S13 and S14). Both LDA and differential abundance analyses showed the same ASVs that were related to certain treatments. ASVs from genera such as Pseudarthrobacter, Methylobacterium-Methylorubrum, Noviherbaspirillum, Pedobacter, Niastella, Ralstonia, SH-PL14 and ASV from Microscillaceae family were abundant in the half-water treatment relative to the full-water treatment in both LDA and differential abundance analyses. On the other hand, ASVs from genera such as Rhodobacter, Sphingopyxis, Luteimonas, Caulobacteraceae, ASV from Comamonadaceae, Verrucomicrobiaceae families, and ASV from Burkholderiales, Cyanobacteria orders were abundant in the full-water treatment relative to the half-water treatment in both LDA and differential abundance analyses. When comparing differential abundance analysis results with that from LDA we only considered ASVs with the probabilities larger than 0.025.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Difference in abundances of ASVs in the half-water treatment relative to the full-water (dashed middle line).

Dots represent the differential abundance coefficient and the error bars are standard errors. The taxa shown are only those that are significant after a p-value correction with the FDR set to 0.05. Plots were produced using corncob R package [50]. The dots colored in red show taxa that were abundant in half-water related LDA topics. The dots colored in blue show taxa that were abundant in full water related LDA topics. Note that the comparison was made with the Fig 7 where only ASVs with p > 0.025 are shown. The second half of the plot with other significant taxa is shown in S12 Fig.

https://doi.org/10.1371/journal.pcbi.1011075.g008

Research data for this article

Data used in this analysis come from our greenhouse experiment. Unprocessed sequence data have been deposited to the NCBI Sequence Read Archive under the projects PRJNA 780613 (generation 0 samples) and PRJNA 780954 (generation 1 samples). The code as wells as data files containing the distribution of taxa in LDA topics are publicly available at https://github.com/lanl/ML-for-microbiomes.