Metagenomic Functional Shifts to Plant Induced Environmental Changes

The Vaccinium angustifolium (wild blueberry) agricultural system involves transformation of the environment surrounding the plant to intensify plant propagation and to improve fruit yield, and therefore is an advantageous model to study the interaction between soil microorganisms and plant–host interactions. We studied this system to address the question of a trade-off between microbial adaptation to a plant-influenced environment and its general metabolic capabilities. We found that many basic metabolic functions were similarly represented in bulk soil and rhizosphere microbiomes overall. However, we identified a niche-specific difference in functions potentially beneficial for microbial survival in the rhizosphere but that might also reduce the ability of microbes to withstand stresses in bulk soils. These functions could provide the microbiome with additional capabilities to respond to environmental fluctuations in the rhizosphere triggered by changes in the composition of root exudates. Based on our analysis we hypothesize that the rhizosphere-specific pathways involved in xenobiotics biodegradation could provide the microbiome with functional flexibility to respond to plant stress status.

Introduction

The rhizosphere is a hotspot of plant–microbiome interactions within soil environments. It is occupied by highly diverse microbial communities which are structurally and functionally influenced by plant and soil type (Philippot et al., 2013). Many previous studies of the metabolic capabilities of plant-associated microbiomes have identified specific functions linked to plant–microbiomes interaction, such as cell motility and root adhesion, metabolism of nitrogen, carbohydrates and vitamins, and xenobiotic degradation (Bulgarelli et al., 2015; Yan et al., 2017; Kamutando et al., 2018). However, it is not clear how these rhizosphere-specific (RS) functions are incorporated into the overall metabolic capabilities of soil microbiomes in general. In particular, it remains unclear whether microbial adaptation to the rhizosphere environment is associated with decreased general metabolic capabilities. The presence of this association would suggest that rhizosphere-associated microbes have some level of plant dependency whereby they rely on plant-produced compounds for their metabolism.

Herein, we conducted a co-occurrence analysis of the pathways found in bulk soil and rhizosphere metagenomes to identify niche-specific sub-networks and their interactions. We hypothesized that if there were a trade-off between rhizosphere and bulk soil metabolic capabilities then the pathway co-occurrence network would have a modular structure with strong positive interactions between niche-specific functions within each sub-network and negative interactions between rhizosphere and bulk soil-specific sub-networks. We were also interested in the identification of rhizosphere specific (RS) functional sub-networks that did not have strong negative interactions with bulk soil-specific (BSS) functional sub-networks, therefore representing an increment in functional repertoire of plant-associated microbiome.

Vaccinium angustifolium (wild blueberry) management is an attractive system to study plant–microbiome interaction since, in contrast to most agricultural systems that involve cultivation on plants in exogenous environments, this system involves transformation of the environments surrounding the plant to intensify plant propagation and to improve fruit yield (Hall et al., 1979; Eaton, 1988; Bell et al., 2009; Drummond et al., 2009). As a result the plant genotype and its native location become the only factors that remaining unchanged in wild blueberry managed habitats. In contrast, the management of the land, plant and soil affects land cover, wild blueberry fruit production and soil properties, which in turn can affect the diversity and structure of the bulk soil and plant-associated microbiome. Consequently, this system is an advantageous model to study the effects of environmental factors on soil and plant-associated microbiomes, the interaction between soil microorganisms and plant–host, and the functional differentiation of rhizosphere microbiomes.

In our previous studies (Yurgel et al., 2017, 2018) the application of both 16S and 18S rRNA amplicon sequencing provided us with some level of understanding of how deterministic factors, such as soil and plant properties shape the wild blueberry microbiome structure and define plant–microbiome interaction. We detected a significant effect of management on the structure of bacterial microbiome in the wild blueberry rhizosphere (Yurgel et al., 2017). Community correlation networks analysis identified several potential hub taxa with important roles in soil fertility and/or plant–microbe interaction and showed that bacterial and eukaryotic interactions became more complex along the soil-endosphere continuum, likely due to the increasing influence of host–plant on microbiome function (Yurgel et al., 2018).

However, the role of environmental factors in the functional assemblage of the wild blueberry microbiome has not been explored. Additionally, the indication that the structure of wild blueberry rhizosphere microbiome was affected by management raised a question of whether the observed taxonomic differences reflected niche-specific microbial functioning. For example, differences in the taxonomic composition of nematode-associated microbiomes did not reflect differences in their functional capabilities (Cheng et al., 2013), likely due to the substantial functional redundancy in microbial communities. Moreover, past work investigating the role of diversity in the selection of bacterial taxa and functions in soil and rhizosphere, suggested that functional traits were a key to the assembly of the rhizosphere microbiome (Yan et al., 2017). To explore this subject further we compared the functional profiles of rhizosphere microbiomes from managed and forest grown plants, which differ significantly in their taxonomic profiles (Yurgel et al., 2017). We were especially interested whether or not the taxonomic differences in rhizosphere microbiomes could be related to their functional capabilities.

Materials and Methods

Sample Collection and DNA Isolation

The wild blueberry root samples, rhizosphere and bulk soil used in this study were collected in August 2015 (Yurgel et al., 2017). The same DNA isolated from bulk and rhizosphere soils from Collingwood (45°36′53″ N 63°56′32″ W) and Debert (45°25′36″ N 63°29′47″ W) sampling sites (Supplementary Table S1) used for 16S and 18S rRNA amplicon analysis (Yurgel et al., 2017, 2018) was used for shotgun metagenome sequencing. The samples collection, processing and DNA isolation is described in Yurgel et al. (2017). In short bulk field soil samples were taken from the top 20 cm of the topsoil layer directly under blueberry plants in an X-shaped pattern with at least 30 m of distance between each collection point. The litter layer was removed and the samples were transported into the laboratory on ice, and immediately stored at −20°C for chemical characterization. After transportation to the lab, 5 g of each soil sample was sieved (2 mm) and immediately stored at −86°C until processing for DNA isolation. DNA was isolated from 0.25 g of the frozen filtered soil samples.

For the collection of rhizosphere samples wild blueberry rhizomes and associated roots were extracted from soil, vigorously shaken, placed in sterile bags, and transported to the laboratory on ice. Root samples were processed immediately after transportation to the laboratory. Blueberry roots were placed in a Falcon tube (50 ml) with 10% glycerol (40 ml) and vortexed until adhering soil had been visibly removed from the root. The roots were removed and the soil suspension was centrifuged at 3,000 × g for 15 min. The supernatant was decanted, and the soil pellets were transferred into 1.5-ml Eppendorf tubes and stored at −86°C. DNA was isolated from 0.25 g (wet weight) of each frozen rhizosphere soil sample.

Shotgun Library Construction and Sequencing

One nanogram of DNA for each sample was subjected to Nextera XT (Illumina) library preparation. This was done as per the manufacturer’s instructions except the clean-up and normalization stages, which were completed using Just-a-Plate 96 PCR Purification and a Normalization Kit (Charm Biotech). Equal amounts of all barcoded samples were then pooled and sequenced in a shared 150+150 bp paired-end NextSeq run (Illumina High-Output v2 kit). The data generated for this study is available at the short read archive under BioProject PRJNA484230.

Shotgun Metagenomics Pre-processing

The Microbiome Helper (Comeau et al., 2017) metagenomics workflow v2 was followed to process the shotgun metagenomics data with modifications to normalize functional abundance by average genome size within a sample. Briefly, raw paired-end reads were filtered using the KneadData v0.6.1 pipeline parallelized with GNU Parallel v20170322 (Tange, 2011). The filtering step included first trimming reads with Trimmomatic v0.36 (Bolger et al., 2014) in sliding windows of four nucleotides (nt) with a minimum quality score of 30. Trimmed sequences below 50 nt were removed. The next filtering step involved running Bowtie2 v 2.3.4.2 (Langmead and Salzberg, 2012) to remove contaminant reads that mapped human (hg38), PhiX, or Vaccinium corymbosum W8520 (Gupta et al., 2015) genomes with the options “—very-sensitive” and “–dovetail.” HUMAnN2 0.11.2 (Franzosa et al., 2018) was run on these filtered reads to identify the number of reads per kilobase (RPK) of UniRef50 gene families (Suzek et al., 2015) in each sample. UniRef50 gene families were then converted into KEGG orthologs using HUMAnN2’s built in mapping files. KEGG ortholog abundances were then normalized by the number of genome equivalents (average genome size/library size) found in each sample using Microbe Census (Nayfach and Pollard, 2015) to get reads per kilobase per genome equivalent (RPKG). Normalized KEGG orthologs were then mapped to KEGG pathways and KEGG modules using MinPath (Ye and Doak, 2009) and the PICRUSt2 script pathway_pipeline.py (Douglas et al., 2019) using the “–no_regroup” option and mapping files included with PICRUSt2 v2.1.4-b. MetaPhlAn2 v2.7.62 (Truong et al., 2015) was run with default options within HUMAnN2 to identify the relative abundances of taxa within each sample. A custom Python script was used to determine the GC-content of the microbiome overall and for the V. corymbosum reference genome.

HUMAnN2 tables with RPKG of metabolic pathways were transformed to BIOM tables and used for the generation of beta-diversity (weighted Bray–Curtis dissimilarity) (Bray and Curtis, 1957) metrics and analysis of variations in sample groupings explained by weighted Bray–Curtis beta-diversity dissimilarity (Adonis tests, 999 permutations) using QIIME wrapper scripts (Caporaso et al., 2010). An important assumption of this approach is that the read depth per sample is sufficiently high to accurately estimate the abundance of the gene families underlying these metabolic pathways. Analysis of taxonomic and functional profiles was performed using the STAMP software package (Parks et al., 2014). Corrected p-values (q-values) were calculated based on Benjamini–Hochberg FDR multiple-test correction. Features with (Welch’s t-test) q-value < 0.05 were considered significant and were thus retained.

Co-occurrence Network Construction and Analysis

The co-occurrence analysis was performed using the CCREPE (Compositionality Corrected by REnormalization and PErmutation) R package (Schwager and Huttenhower, 2016, Unpublished) with 1000 bootstrap iterations and default setting. This package has previously been used to construct co-occurrence networks from microbial sequencing data (Vazquez-Baeza et al., 2016; Yurgel et al., 2018). This network uses a novel similarity measure, the N-dimensional checkerboard score (NC-score) (Stone and Roberts, 1990), which is particularly appropriate to compositions derived from microbial community sequencing data. First, the co-occurrence and co-exclusion patterns in the samples were scored. The results were filtered to remove non-statistically significant relationships. We generated the network based on strong correlations with p-values < 0.001. The networks were visualized with Cytoscape (Shannon et al., 2003) and were represented as graphs with microbial functions as vertices/nodes and the edges as interaction types.

Results

Sequencing and Data Processing

A total of 44,641,544, 70,129,920, 19,394,324, and 29,124,986 raw reads were obtained from 8 field bulk, 8 field rhizosphere, 4 forest bulk, and 3 forest rhizosphere soil samples, respectively. After assembly, quality control filtering and removal of artificial sequences produced by sequencing artifacts and plant sequences 9,412,563, 9,245,366, 2,864,804, and 3,756,327 high-quality reads were retained in field bulk soil, field rhizosphere, forest bulk, and forest rhizosphere metagenomes respectively (Supplementary Table S2). Interestingly, forest rhizosphere metagenomes exhibited distinct genome properties. It had a significantly lower (Tukey’s test p < 0.05) mean guanine and cytosine (GC) value (∼52%) compared to field rhizosphere (∼57%) and both field (∼59%) and forest (59%) bulk metagenomes.

Taxonomic Composition of Bulk Soil and Rhizosphere Metagenomes

Using a unique clade-specific marker gene approach (Truong et al., 2015) we identified a total of 22 bacterial taxa across all samples. Proteobacteria, Acidobacteria, and Actinobacteria were the predominant bacterial phyla in both bulk soil and rhizosphere samples. A similar taxonomic profile of soil microbiomes obtained by shotgun metagenome analysis was previously reported (Fonseca et al., 2018; Kamutando et al., 2018). A versatile soil bacteria capable of carbon and nitrogen fixation, Rhodopseudomonas palustris, was the most abundant bacterial species represented by >14% of all high-quality reads (Supplementary Table S3) and Rhodococcus erythropolis capable of degrading, toluene, naphthalene, herbicides, and other environmental pollutants (Curragh et al., 1994), was represented by ∼1% of all high-quality reads. We also identified Pseudomonas stutzeri (Parte and Kharat, 2019) a bacterial species capable of denitrifying and potentially degrading insecticides, as well as the nitrogen fixer Bradyrhizobium japonicum, and Granulicella mallensis (Rawat et al., 2012) a versatile heterotroph capable of degrading plant-based carbon polymers (Supplementary Table S3). Interestingly, the human pathogen Mycobacterium intracellulare was also identified, which corroborates confirming previous reports suggesting soil as its potential environmental reservoir (Wallace et al., 2013; Honda et al., 2016).

Functional Characteristics of Bulk and Rhizosphere Microbiomes

In total, 5321 KEGG orthologs (KOs) comprising 200 KEGG pathways, 304 KEGG modules and 40 KEGG categories were identified in the study. Among them amino acid metabolism, carbohydrate metabolism, metabolism of cofactors and vitamins, metabolism of terpenoids and polyketides, lipid metabolism, and xenobiotics biodegradation and metabolism were the top KEGG functional categories identified representing 21, 14, 10, 9, 8 and 7% of all high-quality reads, respectively. The top 18 functional categories are shown in Supplementary Figure S1.

An analysis of the strength and statistical significance of sample groupings (Adonis test) indicated that niche (bulk soil and rhizosphere) environments were associated with differences in their microbiome’s functional composition. We found differences in the functional characteristics of each microbiome at the KEGG modules, pathways and KO levels (Table 1). When the microbiomes from managed and natural habitats were analyzed together by variations in sample groupings (Adonis test) explained by weighted Bray–Curtis beta-diversity dissimilarities, ∼20% (R² = 0.201, p < 0.001), 17% (R² = 0.174, p < 0.001), and 17% (R² = 0.17460, p < 0.05) of functional variations between bulk soil and rhizosphere communities were detected for the modules, pathways and KO levels respectively. More specifically 51 modules were in differential relative abundance between bulk soil and rhizosphere samples. A number of sugar biosynthetic modules as well as fatty acid and amino acid biosynthetic modules were enriched in bulk soils compared to rhizosphere samples. On the other hand the rhizosphere samples were enriched for amino acid and sugar transport system modules. The 10 most relatively abundant modules overrepresented in bulk soil or in rhizosphere are shown in Supplementary Figure S2.

TABLE 1

Table 1. Variation in sample groupings as explained by Bray–Curtis dissimilarity.

To further refine the functional potential of the bulk soil and rhizosphere microbiomes, the relative abandances of pathways were compared. In total 35 pathways were differentially represented between bulk soil and rhizosphere (Supplementary Table S4). In particular, the pathways constituting carbohydrate, glycan, nucleotide, and vitamin metabolism, as well as folding, sorting and degradation, and replication and repair were overrepresented in the bulk soil compared to rhizosphere samples. In contrast, the pathways involved in degradation of lysine, xenobiotics and plant-derived terpenoids, as well as tyrosine and phenylalanine metabolism and bacterial chemotaxis and flagellar assembly were overrepresented within the field rhizosphere samples (Figure 1).

FIGURE 1

Figure 1. KEGG pathways that were at differential relative abanances between bulk soil and rhizosphere samples. Only pathways with mean relative frequencies > 0.5% are shown. Corrected p-values (q-values) were calculated based on Benjamini–Hochberg FDR multiple test correction. Features with (Welch’s t-test) q-value < 0.05 were considered significant and were thus retained. ^*Pathways overrepresented in rhizosphere.

In total 427 KOs were differentially enriched in field bulk soils and rhizosphere. The 15 most abundant KOs overrepresented within each soil-types are shown in Supplementary Figure S3. The rhizosphere-specific (RS) KOs included K01075, which is found in pathways involved in the biosynthesis of secondary metabolites and benzoate degradation, and K04073, which is found in the pathways involved in the degradation of butanoate, benzoate, xylene, dioxin, and other aromatic compounds. The putrescine transport system permease protein, K11074, and arginine decarboxylase K01584, which is involved in arginine and proline metabolism were also overrepresented in rhizosphere samples. The putative ABC transport system ATP-binding protein, K02003, and the cold shock protein (beta-ribbon, CspA family), K03704 were the most abundant bulk-soil-specific (BSS) KOs.

Managed vs. Natural Habitats

The analysis of strength and statistical significance of sample groupings by management indicated that this parameter did not significantly influence the functional characteristics of rhizosphere and bulk soil microbiomes at the modules, pathways, or KO levels (Table 1). We also did not identify any microbiome functions at any functional level that were differently represented between the rhizosphere of managed and unmanaged plants or between the bulk soils from managed and unmanaged plant habitats.

Functional Network Interaction

We generated a co-occurrence network by correlating the relative abundances between pathways identified in samples from the bulk field and forest soil and the rhizosphere of managed and unmanaged plants. The final co-occurrence network contained 92 pathways forming 49 negative and 180 positive correlations (Supplementary Table S5). This co-occurrence network incorporated the pathways represented by ∼73% of all high-quality reads obtained in the study (Supplementary Table S6). We used the “edge-weighted spring embedded” layout to visualize this interaction. This layout allowed us to visually identify functional clusters by pulling positive correlations (blue) together and pushing negative correlations (red) apart (Figure 2). The functional co-occurrence network was modular and contained four well-defined sub-networks. Two of these sub-networks, Sub-networks 1 and 2, had strong negative correlation with each other, while two others, Sub-networks 3 and 4, were loosely interconnected within the network.

FIGURE 2

Figure 2. Co-occurrence network generated based on pathway relatively abundances in term of reads per kilobase per genome equivalent (RPKG) within bulk soil and rhizosphere samples. The size of the node is proportional to each pathway’s relative abundance across all samples. The lines (i.e., edges) connecting nodes represent a co-occurrence relationship that can be either positive (blue) or negative (red). The intensity of the color and the length of the edges represent the strength of relationship. The positions of the nodes within modules were manually adjusted for better visualization.

Sub-Network 1

It contained 28 pathways with 13 pathways overrepresented in bulk soil (Supplementary Table S6). It had 39 negative interactions with outside nodes and 55 positive interactions within the sub-network with average positive and negative interaction scores of 0.62 and −0.63, respectively. This sub-network was the largest sub-network in the co-occurrence network and it encompassed ∼32% of total high-quality reads identified in the entire microbiome. It contained major functional categories essential for microbiome metabolism and function. The most relatively abundant functional categories found in this sub-network included the metabolism of terpenoids and polyketides (biosynthesis of antimicrobials vancomycin and ansamycins), carbohydrate metabolism, genetic information processing, metabolism of cofactors and vitamins, environmental information processing and amino acid metabolism (Table 2). Since nearly half of the pathways in Sub-network 1 were significantly overrepresented in bulk soil, this sub-network was considered BSS. Many pathways from this sub-network formed strong negative correlations with pathways from Sub-network 2 producing 32 negative interactions.

TABLE 2

Table 2. Distribution of KEGG functional categories in the pathway co-occurrence network.

Sub-Network 2

It had 33 negative and 66 positive interactions with average positive and negative interaction scores of 0.62 and −0.62, respectively. It was represented by 17 pathways containing ∼13% of total high-quality reads and was considered as RS since 11 pathways found in this sub-network were significantly overrepresented in rhizosphere (Supplementary Table S6). Sub-network 2 contained a smaller subset of functional categories, some of which might be involved in microbiome responses to environmental factors. For example, the functional category xenobiotics biodegradation was one of the most relatively abundant functional categories in Sub-network 2 and was overrepresented in the sub-network compared to all other sub-networks (Table 2). Additionally, the functional category amino acid metabolism, comprising amino acid degradation and the metabolism of phenolic amino acids, was the most relatively abundant functional category in Sub-network 2. Geraniol degradation was also a part of this sub-network.

The pathways involved in the degradation of aminobenzoate, caprolactam, benzoate, xylene, polycyclic aromatic hydrocarbon (PAH) were significantly overrepresented in the rhizosphere microbiome and were a part of RS Sub-network 2 identified in the co-occurrence network (Figure 3). Bisphenol degradation and the metabolism of xenobiotics by cytochrome P450 were also part of the Sub-network 2 but they could not be attributed to RS pathways by statistical analysis (Supplementary Table S4). However these pathways might be considered as a part of RS functions considering their significant association with several niche-specific pathways. For example, bisphenol degradation (ko00363) was positively correlated with RS metabolism of phenylalanine (ko00360) and aminobenzoate degradation (ko00627), and the metabolism of xenobiotics by cytochrome P450 (ko00980) was positively correlated with RS PAH degradation (ko00624) (Figure 3 and Supplementary Table S3). With the exception of the metabolism of xenobiotics by cytochrome P450 and degradation of bisphenol and PAH, which were represented with 36, 42, and 35% of KOs, respectively, the other pathways involved in xenobiotics, terpenoids, and polyketides metabolism found in Sub-network 2 and were represented with at least 68% of KOs (Supplementary Table S7).

FIGURE 3

Figure 3. Selected KEGG pathways of xenobiotics, terpenoids, and polyketides metabolism in microbial metagenomes from bulk and rhizosphere soils from natural and managed habitats. ^*Pathways overrepresented in rhizosphere.

Sub-Network 3

It had 9 negative and 27 positive interactions with average positive and negative interaction scores of 0.61 and −0.62, respectively. It was represented by ∼7.8% of total high-quality reads obtained in the study and contained 15 pathways from KEGG functional categories. Five pathways belonging to this Sub-network, phenylpropanoid biosynthesis, glycerophospholipid, phosphonate and phosphinate metabolism, bacterial samples (Table 2 and Supplementary Table S6). However, the rest of the pathways found in the Sub-network 3 were similarly represented in bulk soil and the rhizosphere suggesting a boarder role of this sub-networks in both plant-associate and free-living soil microbiome function.

Sub-Network 4

It had 3 negative and 23 positive interactions with average positive and negative interaction scores of 0.64 and −0.81, respectively and contained 11 pathways comprised by ∼7.8% of total high-quality reads (Supplementary Table S6). Genetic information processing and nucleotide metabolism were the most relatively abundant categories in the sub-network (Figure 2). With the exception of pyrimidine metabolism all other pathways comprising Sub-network 4 were similarly represented in bulk and rhizosphere soils suggesting important role of these functional capabilities in both bulk soil and rhizosphere microbiomes.

Discussion

In our previous work we used amplicon-based approaches to study microbial communities associated with wild blueberry production system (Yurgel et al., 2017, 2018). This enabled the investigation of how environmental and plant factors affect microbial community structure, provided an integrative view on inter-kingdom interactions in soil and plant associated microbiomes, and identified microbial taxa with potential importance in plant health and production. However, until now the functional potential of the wild blueberry soil microbiome remained largely unknown. In this study we applied shotgun metagenome sequencing to study the functional characteristics of these communities.

As previously reported (Tessler et al., 2017), whole-metagenome sequencing coupled with clade-based taxonomic algorithm (MetaPhlAn2) approach for the taxonomic profiling of microbial communities was far less comprehensive compared to rRNA gene sequencing analysis. For example, only 22 microbial taxa were identified in our study compared to 996 eukaryotic and 6,802 bacterial OTUs identified by 18S and 16S rRNA sequencing (Yurgel et al., 2018). However, in agreement with our previous reports Proteobacteria, Acidobacteria, and Actinobacteria were the predominant bacterial phyla identified in the wild blueberry-associated soil microbiome (Yurgel et al., 2017, 2018).

Nucleotide Composition

Interestingly, GC-content varied greatly between the forest and field rhizosphere metagenomes. It was previously reported that both phylogeny and the environment affected nucleotide composition of microbial communities from diverse ecological niches (Reichenberger et al., 2015). In our study, the difference in GC-content was only detected between field bulk soil and rhizosphere of managed plants. On the other hand, both plant-proximity and management affected the taxonomic composition of wild blueberry microbial communities (Yurgel et al., 2017, 2018) confirming the importance of environments in shaping microbial nucleotide composition.

It is also important to notice that the raw reads obtained from rhizosphere and bulk soil niches also exhibited distinct GC-signature. Both field and forest bulk soil metagenome had significantly higher (p < 0.05) GC-content (∼60%) compared to rhizosphere metagenomes (∼50%). We calculated that GC-content of V. corymbosum scaffold was ∼38%. Therefore, our initial hypothesis was that residual plant tissue contributed to rhizosphere metagenomes decreasing their GC-content. In our analysis we filtered raw reads against genome of V. corymbosum (Bian et al., 2014) to remove plant derived reads. However, after quality control filtering and removal of contaminant sequences forest rhizosphere metagenome still exhibited lower GC-content compared to the other metagenomes, suggesting that the presence of plant DNA in the raw metagenome was not the only reason for the distinct GC-signature of rhizosphere metagenome.

Carbon-Nitrogen Cycling

Microorganisms play a critical role in carbon-nitrogen biogeochemical cycling. In agreement with this, carbohydrate and nitrogen metabolism were the most abundant KEGG functional categories identified in our study. A purple photosynthetic bacterium from the Bradyrhizobiaceae family, Rhodopseudomonas palustris, was the most abundant bacterial species in our samples and is known to have extraordinary metabolic capabilities. It can acquire carbon by catabolism of organic molecules and carbon fixation, use light and organic and inorganic compounds for energy, as well as fix atmospheric di-nitrogen, actively participating in most steps of carbon and nitrogen cycle (Larimer et al., 2004). Several other bacterial taxa involved in nitrogen and carbon cycling from the Bradyrhizobiaceae, Beijerinckiaceae, Burkholderiaceae, and Rhodopseudomonas families were identified in our study. This relatively high abundance of bacteria with potential nitrogen and carbon fixation capabilities might be attributed to the adaptation of the microbiome to low fertility soils typical of wild blueberry habitats (Yurgel et al., 2017). The denitrifying bacteria Pseudomonas stutzeri (Lalucat et al., 2006) was also detected in the microbiome. Although the relative abundance of this bacterium was much lower compared to potential nitrogen fixers, probably because of slow rate of nitrification in acidic soils.

Niche-Specific Functions

The functional characteristics varied significantly between rhizosphere and bulk soil microbiomes. Depending on the level of functional annotation (modules, pathways, or KOs), between 16 and 23% of the variation in community function was niche specific. Carbohydrate metabolism was underrepresented, while the transport of sugars and by product of amino acid degradation, putrescine, were overrepresented in the rhizosphere, potentially due to compensation for decrease in some metabolic capacity of the plant-associated microbiome. Additionally, the degradation of complex organic compounds was overrepresented in the rhizosphere. Previously we showed that the community structure was significantly influenced by the plant, with between 10 and 14% of bacterial community variation being niche-specific (Yurgel et al., 2017, 2018), suggesting that some functional changes between bulk soil and rhizosphere could be attributed to the alteration in microbial community structure.

Functional Stability of Bulk Soil Microbiome

Management was a significant factor influencing the diversity and structure of bulk soil microbiome with 12% of bacterial community variation in bulk soil being attributed to management (Yurgel et al., 2017). However, no significant effect of the management on the bulk soil microbiome function was detected in our study. Furthermore, the functional annotation of bulk soil microbiomes from managed and natural habitats were very similar and no significant differences in the mean relative frequency of functions were detected between them. This data suggested a stronger functional stability of soil microbiome compared to the community structure. Similar functional stability of nematode-associated microbiome was reported previously (Cheng et al., 2013).

Functional Co-occurrence Network Interaction

To extend our view on the functional capabilities within the microbiome the interaction between metabolic pathways was investigated using a co-occurrence network analysis. This approach allowed us to extend the set of niche-specific pathways providing more comprehensive insight on overall microbiome metabolism and its functional specialization. We also performed statistical analysis of differences in the mean relative frequency of pathways between combined rhizosphere and bulk soils samples. As expected, in addition to 35 pathways, which were differentially represented in bulk soil and rhizosphere metagenomes, 14 additional pathways were identified as niche-specific by association. For example, bisphenol degradation was included into RS pathways category by its strong co-occurrence association with phenylalanine and aminobenzoate degradation, which were significantly overrepresented in rhizosphere metagenome.

Xenobiotics Biodegradation and Metabolism

The potential of microbial communities to cycle natural halogens and aromatic compounds in soils is well-documented (Curragh et al., 1994; Weigold et al., 2016; Hlihor et al., 2017; Vergani et al., 2017; Hussain et al., 2018). In agreement with these previous findings, xenobiotics biodegradation and metabolism was among the most relatively abundant functional categories annotated in the metagenome, represented by 7% of all high-quality reads. All but three pathways belonging to this category were identified in our study including degradation of haloalkanes, herbicide atrazine, and aromatic compounds, such as bisphenol, naphthalene, chlorocyclohexane, chlorobenzene, nitrotoluene, styrene, caprolactam, toluene, benzoate, aminobenzoate, PAHs, fluorobenzoate, and dioxin. Herbicide atrazin was traditionally used in management of perennial grasses in wild blueberry fields, but was banned for use in 2013. The samples used in the study were collected in summer 2015, which indicates the long lasting effect of this compound on microbial functional potential.

Around 46% of all high-quality reads annotated as xenobiotics biodegradation and metabolism category comprised pathways which were either significantly overrepresented in rhizosphere or were strongly associated with RS pathways. The enrichment of the pathways involved in degradation of complex organic compounds in the rhizosphere could be attributed to several possible mechanisms. These mechanisms include: an increase in plant-derived complex molecules with the proximity to the plant, the stimulation of microbial biodegradation of organic pollutants by plant rhizo-deposits and inducers, and root exudate facilitated co-metabolic degradation (Haby and Crowley, 1996; Shaw and Burns, 2003; Grandy and Neff, 2008; Toyama et al., 2013).

Overall Functional Structure of Wild Blueberry Microbiome

The pathway co-occurrence network contained four well-defined sub-networks. Sub-network 1 was considered BSS, Sub-network 2 was considered RS, and Sub-networks 3 and 4 was considered niche-independent. The Sub-network 2 comprised many pathways involved in xenobiotics biodegradation and metabolism of terpenoids and polyketides. Many of these functions were negatively correlated with the basic metabolic pathways from Sub-network 1. 32 negative correlations were detected between the sub-networks, giving on average ∼1.9 negative correlation per node in Sub-network 2. Taking into consideration that the majority of pathways in Sub-network 1 were essential for general microbiome function, we can conclude that functions in Sub-network 1 and 2 formed a single unit, where the gain of the specific functions essential for adaptation of the microbiome to the rhizosphere environment were linked to the loss of some basic metabolic functions. For example, glycolysis/gluconeogenesis, pentose phosphate pathway, TCA cycle, carbon fixation, metabolism of cofactors and vitamins and genetic information processing were underrepresented in the rhizosphere. This shift in genome functions might be beneficial for the survival of the microbes in the rhizosphere but might also reduce the ability of these microbes to withstand stresses in bulk soils.

Many basic metabolic functions were similarly represented in bulk soil and rhizosphere. For example, within the carbohydrate metabolism category, starch, sucrose, pyruvate, galactose, ascorbate, glyoxylate, dicarboxylate, aldarate, propanoate, butanoate, and C₅-branched dibasic acid metabolism and pentose and glucuronate interconversions were not niche specific, indicating functional capability of both in bulk soils and rhizosphere microbiome to metabolize a divers set of sugars and organic acids. In general, primarily metabolites (sugars, organic acids, and amino acids) are predominant compounds found in root exudates, while secondary metabolites (polycyclic aromatic compounds and phytohormones) are less abundant (van Dam and Bouwmeester, 2016). However, the composition of root exudates is affected by plant genotype, age, and biotic and abiotic factors (Chaparro et al., 2013; Mhlongo et al., 2018). It was proposed that plant root exudates are derived from the most abundant metabolites available at the time. For example, even small changes in plant water status can impact the composition of root exudates. It was shown that Quercus ilex (holm oak) root exudates contained 71% of secondary metabolites and 81% of primary metabolites under drought stress and the recovery phase, respectively (Gargallo-Garriga et al., 2018). The functional structure of the wild blueberry soil microbiome indicated that the basic metabolic functions necessary for both bulk soil and rhizosphere microbiomes are complemented with rhizosphere-specific functions, which could provide the microbiome with additional functional capabilities to respond to environmental fluctuations in rhizosphere triggered by changes in the composition of root exudates. Interestingly, the majority of the pathways from Sub-network 2 were involved in the degradation of secondary metabolites, such as polycyclic aromatic compounds and terpenoids suggesting the importance of these pathways in microbiome responses to plant stress status.

On the other hand, Sub-networks 3 and 4 were less connected in the network. We hypothesize that this low integration of the pathways comprising Sub-network 3 into general microbiome metabolic functions (Sub-networks 1 and 2) could provide the microbiome with additional functional capabilities for interaction with host–plant, including chemotaxis, plant–pathogen interaction, and lipid and glycan metabolism.

Conclusion

We investigated the functional characteristics of the wild blueberry soil microbiome with shotgun metagenomics sequencing. The metagenome GC-content varied greatly between rhizosphere microbiomes from managed and unmanaged habitats indicating the importance of environments in shaping microbial nucleotide composition. A relatively high abundances of microorganisms with potential nitrogen and carbon fixation capabilities were detected suggesting an adaptation of the microbiome to low fertility soils typical for wild blueberry habitats. In contrast to the rhizosphere metagenome, the bulk soil metagenome exhibited functional stability indicating similar functional repertoire in microbiomes associated with soils from managed and natural habitats. Lower relative abundance in general metabolic functions, such as carbohydrate and amino acid metabolism and higher relative abundance in sugar and putrescine transport, as well as degradation of complex organic compounds were detected in the rhizosphere microbiome. The analysis of the pathway co-occurrence network extended the set of niche-specific functions and provided a better understanding of overall microbiome metabolism and its functional specialization. The strong modular structure of metagenome indicated a potential trade-off between functional adaptations of microorganisms to the rhizosphere environment and its basic metabolic function. Based on our analysis we hypothesize that the rhizosphere-specific pathways involved in biodegradation xenobiotics and terpenoids could provide the microbiome with functional flexibility to respond to plant stress status.

Data Availability

Publicly available datasets were analyzed in this study. This data can be found here: https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA484230.

Author Contributions

SY and ML obtained funding and analyzed and discussed the results. SY designed the study, performed data pre-possessing, and wrote the manuscript. SY and JN performed the bioinformatics analyses. GD assisted with the bioinformatics analyses. SY, JN, GD, and ML participated in the production and the final version of the manuscript.

Funding

This work was supported by Nova Scotia Research Acceleration Grant RA2016-005, Discovery Grant from Natural Sciences and Engineering Research Council of Canada and Wild Blueberry Producers Association of Nova Scotia to SY. Additional support was provided to ML through the Strategic Research Initiatives Fund from Dalhousie University and a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada. GD is funded by a Natural Sciences and Engineering Research Council (NSERC) Alexander Graham Bell Graduate Scholarship (Doctoral). JN is a trainee in the Cancer Research Training Program of the Beatrice Hunter Cancer Research Institute, with funds provided by the Terry Fox Research Institute (TFRI).

Conflict of Interest Statement

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.

Supplementary Material

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

References

Bell, D. J., Rowland, L. J., Zhang, D., and Drummonda, F. A. (2009). Spatial genetic structure of lowbush blueberry, Vaccinium angustifolium, in four fields in Maine. Botany 87, 932–946.

PubMed Abstract | Google Scholar

Bian, Y., Ballington, J., Raja, A., Brouwer, C., Reid, R., Burke, M., et al. (2014). Patterns of simple sequence repeats in cultivated blueberries (Vaccinium section Cyanococcus spp.) and their use in revealing genetic diversity and population structure. Mol. Breed. 34, 675–689. doi: 10.1007/s11032-014-0066-7

CrossRef Full Text | Google Scholar

Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Bray, J. R., and 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

Bulgarelli, D., Garrido-Oter, R., Munch, P. C., Weiman, A., Droge, J., Pan, Y., et al. (2015). Structure and function of the bacterial root microbiota in wild and domesticated barley. Cell Host Microbe 17, 392–403. doi: 10.1016/j.chom.2015.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Caporaso, J. G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F. D., Costello, E. K., et al. (2010). QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336.

Google Scholar

Chaparro, J. M., Badri, D. V., Bakker, M. G., Sugiyama, A., Manter, D. K., and Vivanco, J. M. (2013). Root exudation of phytochemicals in Arabidopsis follows specific patterns that are developmentally programmed and correlate with soil microbial functions. PLoS One 8:e55731. doi: 10.1371/journal.pone.0055731

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, X. Y., Tian, X. L., Wang, Y. S., Lin, R. M., Mao, Z. C., Chen, N., et al. (2013). Metagenomic analysis of the pinewood nematode microbiome reveals a symbiotic relationship critical for xenobiotics degradation. Sci. Rep. 3:1869. doi: 10.1038/srep01869

PubMed Abstract | CrossRef Full Text | Google Scholar

Comeau, A. M., Douglas, G. M., and Langille, M. G. (2017). Microbiome helper: a custom and streamlined workflow for microbiome research. mSystems 2:e0127-16. doi: 10.1128/mSystems.00127-16

PubMed Abstract | CrossRef Full Text | Google Scholar

Curragh, H., Flynn, O., Larkin, M. J., Stafford, T. M., Hamilton, J. T., and Harper, D. B. (1994). Haloalkane degradation and assimilation by Rhodococcus rhodochrous NCIMB 13064. Microbiology 140(Pt 6), 1433–1442. doi: 10.1099/00221287-140-6-1433

PubMed Abstract | CrossRef Full Text | Google Scholar

Douglas, G. M., Maffei, V. J., Zaneveld, J., Yurgel, S. N., Brown, J. R., Taylor, C. M., et al. (2019). PICRUSt2: an improved and extensible approach for metagenome inference. BioRxiv

Google Scholar

Drummond, F., Smagula, J., Annis, S., and Yarborough, D. (2009). Organic Wild Blueberry Production. Orono, ME: Maine Agricultural And Forest Experiment Station.

Google Scholar

Eaton, L. J. (1988). Nutrient Cycling in Lowbush Blueberry Stands. Ph. D. thesis, Dalhousie University, Halifax, NS.

Google Scholar

Fonseca, J. P., Hoffmann, L., Cabral, B. C. A., Dias, V. H. G., Miranda, M. R., de Azevedo Martins, A. C., et al. (2018). Contrasting the microbiomes from forest rhizosphere and deeper bulk soil from an Amazon rainforest reserve. Gene 642, 389–397. doi: 10.1016/j.gene.2017.11.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Franzosa, E. A., McIver, L. J., Rahnavard, G., Thompson, L. R., Schirmer, M., Weingart, G., et al. (2018). Species-level functional profiling of metagenomes and metatranscriptomes. Nat. Methods 15, 962–968. doi: 10.1038/s41592-018-0176-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Gargallo-Garriga, A., Preece, C., Sardans, J., Oravec, M., Urban, O., and Penuelas, J. (2018). Root exudate metabolomes change under drought and show limited capacity for recovery. Sci. Rep. 8:12696. doi: 10.1038/s41598-018-30150-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Grandy, A. S., and Neff, J. C. (2008). Molecular C dynamics downstream: the biochemical decomposition sequence and its impact on soil organic matter structure and function. Sci. Total Environ. 404, 297–307. doi: 10.1016/j.scitotenv.2007.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, V., Estrada, A. D., Blakley, I., Reid, R., Patel, K., Meyer, M. D., et al. (2015). RNA-Seq analysis and annotation of a draft blueberry genome assembly identifies candidate genes involved in fruit ripening, biosynthesis of bioactive compounds, and stage-specific alternative splicing. Gigascience 4:5. doi: 10.1186/s13742-015-0046-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Haby, P. A., and Crowley, D. E. (1996). Biodegradation of 3-chlorobenzoate as affected by rhizodeposition and selected carbon substrates. J. Environ. Qual. 25, 304–310.

Google Scholar

Hall, I. V., Aalders, L. E., and Nickerson, N. L. (1979). The biological flora of Canada. I. Vaccinium angustifolium Ait., sweet lowbush blueberry. Can. Field Nat. 93, 415–430.

Google Scholar

Hlihor, R. M., Rosca, M., Tavares, T., and Gavrilescu, M. (2017). The role of Arthrobacter viscosus in the removal of Pb(II) from aqueous solutions. Water Sci. Technol. 76, 1726–1738. doi: 10.2166/wst.2017.360

PubMed Abstract | CrossRef Full Text | Google Scholar

Honda, J. R., Hasan, N. A., Davidson, R. M., Williams, M. D., Epperson, L. E., Reynolds, P. R., et al. (2016). Environmental nontuberculous mycobacteria in the Hawaiian Islands. PLoS Negl. Trop. Dis. 10:e0005068. doi: 10.1371/journal.pntd.0005068

PubMed Abstract | CrossRef Full Text | Google Scholar

Hussain, I., Aleti, G., Naidu, R., Puschenreiter, M., Mahmood, Q., Rahman, M. M., et al. (2018). Microbe and plant assisted-remediation of organic xenobiotics and its enhancement by genetically modified organisms and recombinant technology: a review. Sci. Total Environ. 628-629, 1582–1599. doi: 10.1016/j.scitotenv.2018.02.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamutando, C. N., Vikram, S., Kamgan-Nkuekam, G., Makhalanyane, T. P., Greve, M., Le Roux, J. J., et al. (2018). The functional potential of the rhizospheric microbiome of an invasive tree species, Acacia dealbata. Microb. Ecol. 77, 191–200. doi: 10.1007/s00248-018-1214-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Lalucat, J., Bennasar, A., Bosch, R., Garcia-Valdes, E., and Palleroni, N. J. (2006). Biology of Pseudomonas stutzeri. Microbiol. Mol. Biol. Rev. 70, 510–547. doi: 10.1128/mmbr.00047-05

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

Larimer, F. W., Chain, P., Hauser, L., Lamerdin, J., Malfatti, S., Do, L., et al. (2004). Complete genome sequence of the metabolically versatile photosynthetic bacterium Rhodopseudomonas palustris. Nat. Biotechnol. 22, 55–61. doi: 10.1038/nbt923

PubMed Abstract | CrossRef Full Text | Google Scholar

Mhlongo, M. I., Piater, L. A., Madala, N. E., Labuschagne, N., and Dubery, I. A. (2018). The chemistry of plant-microbe interactions in the rhizosphere and the potential for metabolomics to reveal signaling related to defense priming and induced systemic resistance. Front. Plant Sci. 9:112. doi: 10.3389/fpls.2018.00112

PubMed Abstract | CrossRef Full Text | Google Scholar

Nayfach, S., and Pollard, K. S. (2015). Average genome size estimation improves comparative metagenomics and sheds light on the functional ecology of the human microbiome. Genome Biol. 16:51. doi: 10.1186/s13059-015-0611-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Parks, D. H., Tyson, G. W., Hugenholtz, P., and Beiko, R. G. (2014). STAMP: statistical analysis of taxonomic and functional profiles. Bioinformatics 30, 3123–3124. doi: 10.1093/bioinformatics/btu494

PubMed Abstract | CrossRef Full Text | Google Scholar

Parte, S. G., and Kharat, A. S. (2019). Aerobic degradation of clothianidin to 2-Chloro-methyl Thiazole and Methyl 3-(Thiazole-yl) methyl guanidine produced by Pseudomonas stutzeri smk. J. Environ. Public Health 2019:4807913. doi: 10.1155/2019/4807913

PubMed Abstract | CrossRef Full Text | Google Scholar

Philippot, L., Raaijmakers, J. M., Lemanceau, P., and van der Putten, W. H. (2013). Going back to the roots: the microbial ecology of the rhizosphere. Nat. Rev. Microbiol. 11, 789–799. doi: 10.1038/nrmicro3109

PubMed Abstract | CrossRef Full Text | Google Scholar

Rawat, S. R., Mannisto, M. K., Bromberg, Y., and Haggblom, M. M. (2012). Comparative genomic and physiological analysis provides insights into the role of Acidobacteria in organic carbon utilization in Arctic tundra soils. FEMS Microbiol. Ecol. 82, 341–355. doi: 10.1111/j.1574-6941.2012.01381.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Reichenberger, E. R., Rosen, G., Hershberg, U., and Hershberg, R. (2015). Prokaryotic nucleotide composition is shaped by both phylogeny and the environment. Genome Biol. Evol. 7, 1380–1389. doi: 10.1093/gbe/evv063

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaw, L. J., and Burns, R. G. (2003). Biodegradation of organic pollutants in the rhizosphere. Adv. Appl. Microbiol. 53, 1–60. doi: 10.1016/s0065-2164(03)53001-5

CrossRef Full Text | Google Scholar

Stone, L., and Roberts, A. (1990). The checkerboard score and species distributions. Oecologia 85, 74–79. doi: 10.1007/BF00317345

PubMed Abstract | CrossRef Full Text | Google Scholar

Suzek, B. E., Wang, Y., Huang, H., McGarvey, P. B., Wu, C. H., and UniProt, C. (2015). UniRef clusters: a comprehensive and scalable alternative for improving sequence similarity searches. Bioinformatics 31, 926–932. doi: 10.1093/bioinformatics/btu739

PubMed Abstract | CrossRef Full Text | Google Scholar

Tange, O. (2011). GNU parallel – the command-line power tool. USENIX Mag. 36, 42–47.

Google Scholar

Tessler, M., Neumann, J. S., Afshinnekoo, E., Pineda, M., Hersch, R., Velho, L. F. M., et al. (2017). Large-scale differences in microbial biodiversity discovery between 16S amplicon and shotgun sequencing. Sci. Rep. 7:6589. doi: 10.1038/s41598-017-06665-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Toyama, T., Ojima, T., Tanaka, Y., Mori, K., and Morikawa, M. (2013). Sustainable biodegradation of phenolic endocrine-disrupting chemicals by Phragmites australis-rhizosphere bacteria association. Water Sci. Technol. 68, 522–529. doi: 10.2166/wst.2013.234

PubMed Abstract | CrossRef Full Text | Google Scholar

Truong, D. T., Franzosa, E. A., Tickle, T. L., Scholz, M., Weingart, G., Pasolli, E., et al. (2015). MetaPhlAn2 for enhanced metagenomic taxonomic profiling. Nat. Methods 12, 902–903. doi: 10.1038/nmeth.3589

PubMed Abstract | CrossRef Full Text | Google Scholar

van Dam, N. M., and Bouwmeester, H. J. (2016). Metabolomics in the rhizosphere: tapping into belowground chemical communication. Trends Plant Sci. 21, 256–265. doi: 10.1016/j.tplants.2016.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Vazquez-Baeza, Y., Hyde, E. R., Suchodolski, J. S., and Knight, R. (2016). Dog and human inflammatory bowel disease rely on overlapping yet distinct dysbiosis networks. Nat. Microbiol. 1:16177. doi: 10.1038/nmicrobiol.2016.177

PubMed Abstract | CrossRef Full Text | Google Scholar

Vergani, L., Mapelli, F., Zanardini, E., Terzaghi, E., Di Guardo, A., Morosini, C., et al. (2017). Phyto-rhizoremediation of polychlorinated biphenyl contaminated soils: an outlook on plant-microbe beneficial interactions. Sci. Total Environ. 575, 1395–1406. doi: 10.1016/j.scitotenv.2016.09.218

PubMed Abstract | CrossRef Full Text | Google Scholar

Wallace, R. J.Jr., Iakhiaeva, E., Williams, M. D., Brown-Elliott, B. A., Vasireddy, S., Vasireddy, R., et al. (2013). Absence of Mycobacterium intracellulare and presence of Mycobacterium chimaera in household water and biofilm samples of patients in the United States with Mycobacterium avium complex respiratory disease. J. Clin. Microbiol. 51, 1747–1752. doi: 10.1128/JCM.00186-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Weigold, P., El-Hadidi, M., Ruecker, A., Huson, D. H., Scholten, T., Jochmann, M., et al. (2016). A metagenomic-based survey of microbial (de)halogenation potential in a German forest soil. Sci. Rep. 6:28958. doi: 10.1038/srep28958

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, Y., Kuramae, E. E., de Hollander, M., Klinkhamer, P. G., and van Veen, J. A. (2017). Functional traits dominate the diversity-related selection of bacterial communities in the rhizosphere. ISME J. 11, 56–66. doi: 10.1038/ismej.2016.108

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, Y., and Doak, T. G. (2009). A parsimony approach to biological pathway reconstruction/inference for genomes and metagenomes. PLoS Comput. Biol. 5:e1000465. doi: 10.1371/journal.pcbi.1000465

PubMed Abstract | CrossRef Full Text | Google Scholar

Yurgel, S. N., Douglas, G. M., Comeau, A. M., Mammoliti, M., Dusault, A., Percival, D., et al. (2017). Variation in bacterial and eukaryotic communities associated with natural and managed wild blueberry habitats. Phytobioms J. 1, 102–113. doi: 10.1094/pbiomes-03-17-0012-r

CrossRef Full Text | Google Scholar

Yurgel, S. N., Douglas, G. M., Dusault, A., Percival, D., and Langille, M. G. I. (2018). Dissecting community structure in wild blueberry root and soil microbiome. Front. Microbiol. 9:1187. doi: 10.3389/fmicb.2018.01187

PubMed Abstract | CrossRef Full Text | Google Scholar