COREMIC: a web-tool to search for a niche associated CORE MICrobiome

Datasets used in the study

A diverse set of data composed of 61 samples from two different published datasets and collected from multiple locations (Jesus et al., 2016; Rodrigues et al., 2017) were used for this study. Data were obtained from the NCBI and selected based on the availability of the raw (16S rRNA) sequence data of root-zone bacteria from switchgrass and that for an out-group of reference (native and/or other grasses) plants.

The dataset “Jesus, 2016” (Jesus et al., 2016), PRJEB6704, compared the rhizosphere soil microbial communities associated with restored prairie with three grass crops, namely corn, switchgrass, and mixed prairie grasses. The grasses were grown in fields of Michigan and Wisconsin and were harvested after two and ten years. The V6–V8 region of the 16S rRNA gene was amplified and sequenced using the Roche 454 pyrosequencing. In our study, we used a total of 43 samples (three each from corn, switchgrass, mixed grasses (two years only), and restored prairie grasses grown in Wisconsin and Michigan, and sampled after two and 10 years. Switchgrass grown in Michigan, composed of 4 samples, were collected following 10 years of plant growth.

The dataset “Rodrigues, 2017” (Rodrigues et al., 2017), PRJNA320123, compared the root-zone soil microbial communities associated with switchgrass cultivars: “Alamo” and “Dacotah”. The switchgrass were grown in the greenhouse using soil derived from plots growing Switchgrass (>7 years) near Blacksburg, VA. Switchgrass rhizosphere bacteria were sampled at three different growth stages. The V3–V4 region of the 16S rRNA gene was amplified and sequenced using Illumina MiSeq sequencing. In our study, we used a total of 18 switchgrass samples for Alamo (A) and Dacotah (D) from stages V2 and E3 (4 AV2, 4 DV2, 5 AE3, 5 DE3 = 18).

Overall, these datasets served as a diverse resource (relevant differences are summarized in Fig. 1) to compare the root-zone bacteria and identify core-bacteria associated with switchgrass.

Sequence data analysis and picking of Operational Taxonomic Units (OTU)

For the Rodrigues, 2017 dataset, the OTU table was obtained from previously performed analysis (Rodrigues et al., 2017). For the Jesus, 2016 dataset, quality score (25) and read lengths (150) thresholds were enforced using cutadapt (1.8.1) (Martin, 2011) and an open reference OTU picking (enable_rev_strand_match True) was performed in QIIME v1.8.0 (Caporaso et al., 2010), as previously described (Rodrigues et al., 2015; Rodrigues et al., 2017), to allow comparison with the other dataset. Briefly, uclust (Edgar, 2010) was used to cluster reads into OTUs (97% sequence similarity) and assign taxonomy against the Greengenes reference database version 13.8 (DeSantis et al., 2006; McDonald et al., 2012). Two samples from the Jesus, 2016 dataset were removed from downstream analysis due to very few sequences assigned to OTUs.

Combining two datasets

For each dataset, sequences assigned to OTUs with identical taxonomy were collapsed (summed within each sample) to create a single row for that taxonomy. This method of taxonomic summary is commonly used (e.g., QIIME) and beneficial for working with diverse datasets. Finally, the OTU tables from the two datasets were merged based on identical taxonomy (row-ids). The common OTUs (678 taxa) from the two datasets were selected, converted to biom format and used for further analyses (Fig. 1). The data table was rarefied using a sequence threshold of 1150, and the beta diversity was calculated using Bray-Curtis (Beals, 1984) distance and visualized using Principal Coordinate Analysis (Gower, 2005). Multivariate data analysis methods of MRPP (Mielke, 1984), Permanova (Anderson, 2001) and ANOSIM (Clarke, 1993) were used to identify whether the plant type (switchgrass versus non-switchgrass) were associated with different bacterial communities. Mann–Whitney-U test from QIIME’s group_significance.py was used to identify OTUs that had differential abundance (FDR < 0.05) in switchgrass and non-switchgrass samples.

Core microbiome analysis

To find the set of core OTUs, the samples in the combined OTU table (original data) were first divided into the interest group samples (switchgrass) and out-group samples. The abundance values for each OTU in each sample are then converted to binary (present/absent) values based on whether they are zero or nonzero. For each OTU a one-tailed Fisher’s Exact Test was used to calculate a p-value testing whether an OTU was present in a significantly higher portion in the interest in-group (Switchgrass) compared to the out-group samples (numerous other grass species).

These p-values were corrected for multiple-testing using Benjamini Hochberg. The OTUs with a q-value < 0.05 were then selected to only the OTUs that are present in at least 90% of the interest group samples. Uninformative OTUs (e.g., k_Bacteria;p_;c_;o_;f_;g_;s_) were filtered out and the remaining OTUs were candidates for the core microbiome.

Submitting data to COREMIC

COREMIC needs data and mapping files. Data file(s) provide the abundance (counts) for each OTU (taxonomy) per sample. To allow users to directly use their output of OTU picking methods, COREMIC requires the data file(s) in BIOM format with the taxonomic assignments for each OTU under the “taxonomy” key of the meta-data in the data file(s). Rows with the same taxonomic label will be combined, with their values summed per sample. If multiple data files are uploaded, they will be combined using taxonomic labels, removing any taxonomies that are not present in all the data files.

The mapping file is a tab-delimited file where each sample is represented as a row and contains group (treatment, location, age, etc.) information as columns about the samples from the data file(s). The help section of COREMIC (http://coremic2.appspot.com/help or http://core-mic.com/help) offers details about the input files and provides sample data and group files. The supplementary video (https://oregonstate.box.com/v/coremic-tutorial) offers a tutorial for submitting data to COREMIC.

Implementation of COREMIC

COREMIC receives the input files and parameters from web-tool. The data file is collapsed (sum the sequences of OTUs with identical taxonomy per sample) to have unique taxonomic row ids. Multiple data files, if any, are merged into a single OTU table based on common taxonomic labels. The user preferred normalization is applied and significance is calculated for each OTU followed by multiple testing. Core microbes are identified as those that pass the user-specified criteria (adjusted p-value, presence in interest and out-group). Finally, the results and output files are emailed to the user.

COREMIC and the datasets are available at http://coremic2.appspot.com or http://core-mic.com. Its code is available on github (https://github.com/richrr/coremicro). The web-tool was developed in Python 2.7, and is hosted on Google App Engine. Other requirements include GoogleAppEnginePipeline 1.9.22.1, pyqi 0.3.1, requests 2.10.0, requests-toolbelt 0.6.2, mailjet-rest 1.2.2, biom-format 1.1.2, ete3 3.0.0 (for tree generation–see below for details), webapp2 2.5.2, numpy 1.6.1, matplotlib 1.2.0, jinja2 2.6, ssl 2.7.11. COREMIC is accessible via any internet connected browser and emails the results to the user. The processing times with the default settings after uploading the data are provided in Table S1.

A custom python script generates a phylogenetic tree using the taxonomic labels for each OTU displaying the relationship between the core OTUs obtained from the group of interest and the out-group. This tree is generated using the ete3 3.0.0 library.

Methodological considerations in the use of COREMIC

COREMIC performs a complementary analysis different from that of existing methods by using presence/absence data. For two groups (A and B) it checks whether (pre-determined percentage of) samples from group A have a non-zero value for the OTU. This allows scientists to operate without making assumptions about the PCR-based OTU relative abundances. This is considered a potential advantage of the method because it is unknown whether relative abundance of sequence data is representative of true relative differences between communities. It is well accepted that sequencing depth can affect the occurrence of rare OTUs in both relative and presence/absence data. Relative abundance analysis, however, would bias against rare OTUs, whereas presence absence equally treats abundant and rare OTUs since any count >1 (or as per the user chosen threshold) is treated as present. Using OTUs that are present in at least 90% of the interest group samples along with significance testing to define a core microbiome accounts for possibility that rare OTUs might be missed as being called present in some samples. Therefore, our proposed approach is relatively robust to sequencing depths. Further research, in this regard, will be aimed towards investigating other measures of OTU “presence”, namely the extent of exclusivity, consistency, or abundance of the group that is eventually determined to be a core microbiome.

Sampling plots used in this study were located across a range of diverse environments to help create a backdrop of heterogeneity. While this diversity of habitat conditions ignores the potential for microbe-environment interactions that might be important for the plant-microbial relationship, it has the advantage of being a conservative approach with high veracity for defining a core microbiome regardless of habitat heterogeneity. The locations from which samples were grown (Michigan, Wisconsin, Virginia) were treated as independent to help isolate the overall habitat effect of switchgrass (Werling et al., 2014; Jesus et al., 2016). When the effects of habitat are thought to be habitat specific, researchers can take this into account during the design and analysis using COREMIC.

It is notable that the representation of an outgroup (multiple non-switchgrass species) is an important criteria and choice made by researchers, and is an approach that has both advantages and caveats. By definition, a habitat is defined by its differences from that of other habitats, and therefore the use of the outgroup is an important choice. A counter-argument for the current dataset might argue for exclusion of breeding lines of a cultivated grass (maize) as being unrepresentative of the grass outgroup. In our case, it was thought, a priori, that a diverse set of grasses would provide the best comparison; and no compelling argument was found that supported the exclusion of maize from the analysis. An implicit assumption was also made that the taxonomy of plant species (root-zone habitats) play an important role in determining root-zone microbial communities, an approach supported by extensive findings that different grass species associate with different microbial communities (Kuske et al., 2002; Kennedy et al., 2004; Berendsen, Pieterse & Bakker, 2012; Chaudhary et al., 2012; Turner et al., 2013). So although there is a need for careful consideration of the experimental questions of interest when using COREMIC, this is a common, if not ubiquitous foundation of all experimentation and hypothesis testing. The results provide a statistically valid approach using freely available software to describe and define a core microbiome of switchgrass.

The choice of the outgroup, furthermore, for determining a core microbiome is amenable to choice using deductive reasoning but ultimately limited by available data. This issue almost certainly limits inclusion of many functionally important rhizosphere microbes that could affect the growth of switchgrass. In this study, the proof of concept utilized a conservative approach to highlight the methodology across a diversity of geographies, soil types, and plant ages. The COREMIC tool as well as the multiple methods for defining a core microbiome (e.g., QIIME Caporaso et al., 2010), ISA (Dufrene & Legendre, 1997) will always be defined by the expertise, and the nature of the hypotheses defined and defended by individual researchers.

Core microbes

The individual datasets described in this study had previously focused on identifying abundant microbes and differences due to experimental conditions. The current meta-analysis goes a step further to find common microbiota that are associated with switchgrass across the diverse experimental conditions. In fact, while Lysobacter would be identified as having significantly differential abundance as per the Man Whitney test, it would not be in the top 50 candidates (ranked as per FDR in Table S3) and likely missed from future testing. Clearly, our approach allows us to identify such candidates. The members of the Lysobacter genus, an identified core microbe of switchgrass, are known to live in soil and have been shown to be ecologically important due to their ability to produce exo-enzymes and antibiotics (Reichenbach, 2006). Their antimicrobial activities against bacteria, fungi, unicellular algae, and nematodes have been described (Islam et al., 2005; Jochum, Osborne & Yuen, 2006; Park et al., 2008; Yin, 2010). Strains of this genus, for example, have been used for control of diseases caused by bacteria in rice (Ji et al., 2008) and tall fescue (Kilic-Ekici & Yuen, 2004). Reports of their function thus support the idea that they may play an important role in switchgrass growth and survival. The core microbiome results thus support further research into the role played by this bacterium in the switchgrass rhizosphere.

Similarly, members of the Mesorhizobium genus are well-known diazotrophs (Kaneko et al., 2000) and previously shown to be symbiotically associated with switchgrass (DeAngelis et al., 2010; Bahulikar et al., 2014) and legumes (Weir et al., 2004; Nandasena et al., 2009; Degefu et al., 2013; Guerrouj et al., 2013). Another identified core microbiome taxa, soil-dwelling members of the Chitinophagaceae family are known to have β-glucosidase (Bailey et al., 2013) and Aminocyclopropane-1-carboxylate (ACC) deaminase activities and ability to produce indole-3-acetic acid (IAA) (Madhaiyan et al., 2015). These molecules and enzymes are well known for their effects on plant growth (Zhao, 2010; Van de Poel & Van Der Straeten, 2014). The capacity to degrade cellulose might provide additional and readily available options to aid survival of these bacteria near switchgrass root zones during times of environmental stress. ACC deaminase and IAA production, in contrast, are potent plant growth modulators (Glick, 2014) that could play a role in plant productivity and survival, especially under conditions of plant physiological stress. Though these examples above would need further study, they provide consistent examples describing how a core microorganism could play a role in determining plant function and growth. The power of the approach stems from the ability to identify the core microbes associated with a plant (or other habitat), and that can, with veracity, narrow down potentially important core microbes from otherwise hyperdiverse samples.

From a technological standpoint, it is important to put the current approach into context with research before the metagenomics era. The search and identification of antagonistic plant growth promoting microbes has previously been tedious and labor intensive. Screenings of hundreds of microbes were used to cultivate and identify candidate microbes that might support (or deter) plant growth. In the case of beneficial microbes, even when identified under greenhouse conditions, the beneficial effects rarely translated into plant supportive growth under field growth conditions (Babalola, 2010; Hayat et al., 2010). With the aid of hindsight and new knowledge suggesting the importance of the soil habitat and root-soil interactions in the development of growth promoting plant-microbial relationships, the approach used in this study reverses the focus (from top-down to bottom-up) to search for microbes that appear to already be naturally well-adapted to the root-soil habitats of interest (Trabelsi & Mhamdi, 2013; Souza, Ambrosini & Passaglia, 2015). This process streamlines the search for suitable microbes from a daunting pool of thousands of bacterial taxa. Bacteria and fungi with well-known partnerships with members of the core microbiome, it would be expected, to be more readily adaptable to their native environment. Indeed, the concept of adaptability to an environment has been shown to be true for many types of microbes across the environmental spectrum, and has given rise to the concept of the niche (Lennon et al., 2012). Existing tools (e.g., Corbata, MetaCoMET and QIIME) do not provide statistical significance for the taxa to be a core microbiome compared to the background using presence/absence data. The point of our method is to provide an approach that is definitive, rather than simply stating a set of microbe present in samples. Due to the lack of consensus in the scientific community about what a gold standard protocol/technology is for such studies (mainly due to the relatively early stages of the metagenomics field and technology still in developmental stages), one is limited (in one way or other) to sufficiently utilize the existing data. While our method of combining datasets has its limitation (like the granularity/species level information limitation from 16S, etc.), it still offers a powerful way to (i) understand the taxonomic distribution within samples (ii) mine multiple existing and often-diverse datasets (iii) generate hypothesis for future detailed experiments, and therefore, certainly a preferred alternative than not using existing data. The COREMIC tool provides an alternative and logical approach to help mine available datasets, in the search for core microbiomes associated with habitats that are ecologically and agriculturally important. Finally, each statistical test relies on different assumptions and has different strengths. COREMIC assumes that presence of a microbe, however high or low, can provide meaningful insight into potential host-microbe relationships. Therefore, it provides equal weightage for high and low abundant microbes within a sample. Since microbial abundance can be an important factor in a biological system, we recommend using COREMIC (presence/absence) in complement with other abundance-based methods (e.g., Mann Whitney test, ISA, etc.). Furthermore, while using diverse dataset has its strengths we suggest avoiding datasets where each dataset contains only a single (mutually exclusive) group. The user needs to consider similarities/differences between the datasets and the biological system while choosing an appropriate outgroup for their group of interest. It is up to the user to decide which thresholds make the most sense for their questions and hypothesis, but obviously a stricter threshold will have higher statistical-inferential veracity. Perhaps a combinatorial approach of selecting the top candidates from the different (presence/absence and abundance based) methods, picking microbes that show significant associations by multiple methods, and using the user’s biological expertise might better allow choosing candidates for future testing.