Divergence of Root Microbiota in Different Habitats based on Weighted Correlation Networks

Summary

Network analysis was applied to evaluate the association of various ecological microbial communities, such as soil, water and rhizosphere. Presented here is a protocol on how to use the WGCNA algorithm to analyze different co-occurrence networks that may occur in the microbial communities due to different ecological environments.

Abstract

The root microbiome plays an important role in plant growth and environmental adaptation. Network analysis is an important tool for studying communities, which can effectively explore the interaction relationship or co-occurrence model of different microbial species in different environments. The purpose of this manuscript is to provide details on how to use the weighted correlation network algorithm to analyze different co-occurrence networks that may occur in microbial communities due to different ecological environments. All analysis of the experiment is performed in the WGCNA package. WGCNA is an R package for weighted correlation network analysis. The experimental data used to demonstrate these methods were microbial community data from the NCBI (National Center for Biotechnology Information) database for three niches of the rice (Oryza sativa) root system. We used the weighted correlation network algorithm to construct co-abundance networks of microbial community in each of the three niches. Then, differential co-abundance networks among endosphere, rhizoplane and rhizosphere soil were identified. In addition, the core genera in network were obtained by the "WGCNA" package, which plays an important regulated role in network functions. These methods enable researchers to analyze the response of microbial network to environmental disturbance and verify different microbial ecological response theories. The results of these methods show that the significant differential microbial networks identified in the endosphere, rhizoplane and rhizosphere soil of rice.

Introduction

Microbiome research has important implications for understanding and manipulating ecosystem processes^1,2. Microbial populations are interconnected by interacting ecological networks, whose characteristics can affect the response of microorganisms to environmental changes^3,4. Furthermore, the properties of these networks affect the stability of microbial communities, and are closely associated with soil function⁵. Weighted gene correlation network analysis has now been widely applied for research on the relationship between genes and microbial communities⁶. Previous studies have focused mainly on the associations between networks of different genes or populations and the outside world⁷. However, the differences in correlation networks formed by microbial populations under different environmental conditions have been scarcely investigated. The purpose of the research presented in this paper is to provide insights and details on the rapid implementation of the WGCNA algorithm to construct a co-occurrence network of microbiome samples collected under different environmental conditions. Based on the analysis results, we assessed the composition and differences of the population and further discussed the relationship between different microbial populations. The following basic flow of weighted correlation network algorithm⁸ was applied. First, a similarity matrix needed to be constructed by calculating the Pearson correlation coefficient between the Operational Taxonomic Units (OTU) expression profiles. Then, the parameters of the adjacency functions (the power or the sigmoid adjacency functions) were adopted with a scale-free topology criterion, the similarity matrix was transformed into an adjacency matrix, and each co-occurrence network corresponded to an adjacency matrix. We used average linkage hierarchical clustering coupled with the TOM-based dissimilarity to group OTUs with coherent expression profiles into modules. Further, we calculated the relationship between conservative statistics and the related parameter analysis modules, finally identifying the hub OTU in the module. These methods are particularly suitable for analysis of the differences in network structures among various microbial populations under divergent environmental conditions. In this manuscript, we have described in detail the method of co-expression network development, the analysis of the dissimilarities between the modules, and have provided a brief overview of the steps in the procedure applied to obtain the core species in different module networks.

Protocol

1. Data Download

1. Download the data of the accession PRJNA386367 form the NCBI database. From the data of the accession PRJNA386367, select the rhizosphere, rhizoplane, and endosphere microbiome data from rice plants grown for 14 weeks in a submerged rice field in Arbuckle, California in 2014.
   ​NOTE: The rhizosphere, rhizoplane, and endosphere microbiome data were presented by the OTUs table in accession PRJNA386367.

2. Optimal power value determination

NOTE: The WGCNA package contains all of the following functional parameters. WGCNA is an R package for weighted correlation network analysis. The key command lines refer to the Supplement S1.

1. In the R language environment, open the Rstudio software and install the WGCNA package.
2. Load the data and use the goodSamplesGenes function to check the correctness of the data. Execute the command lines:
   "gsg = goodSamplesGenes(datExpr0, verbose = 3)
   gsg$allOK "
   Click Run.
3. Check for outliers and store samples that meet the requirements. When the check result is TRUE, continue to the next step. Save the result.
4. Use the PickSoftThreshold function to calculate the scale-free index R² of the two groups of the data under different power values. Execute the command line:
   "sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)"
   Click Run.
5. Visualize the results (Figure 1). Execute the command line:
   "plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
      xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
      main = paste("ES_Scale independence"));
   text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
      labels=powers,cex=cex1,col="red");
   abline(h=0.9,col="red")
   plot(sft$fitIndices[,1], sft$fitIndices[,5],
      xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
      main = paste("ES_Mean connectivity"))
   text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")"
   Click Run.
   NOTE: The premise of the weighted correlation network algorithm is that the established co-expression network structure conforms to the standards of the scale-free topology criterion, increasing its robustness. A scale-free index closer to 1 indicates a network structure that is closer to the scale-free network.
6. Select the power value when the scale-free index R² squared greater than 0.9 and proceed to the next step of analysis.
   ​NOTE: When the scale-free index is close to 1, the network structure is closer to the scale-free network. When analyzing two or more networks, it is necessary to choose to make each network close to the power value of the scale-free network to satisfy the comparability between the co-expressed networks.

3. Construction of a co-expression network and module identification

NOTE: Based on the above calculated power value, the co-occurrence network is constructed. The key command lines refer to the Supplement S2.

1. Use the adjacency function in the WGCNA package to add signed parameters for the construction of a symbolic co-occurrence network. Execute the command line:
   "adjacency = adjacency(datExpr0, power = softPower)"
   Click Run.
2. Apply the TOM-similarity function to develop a topological overlapping network and calculate the dissimilarity network. Execute the command line:
   "TOM = TOMsimilarity(adjacency);
   dissTOM = 1-TOM"
   Click Run.
   NOTE: The signed parameter was added to set the topology overlap network type.
3. Use the hclust function to select the average linkage hierarchical clustering method for hierarchical clustering. Execute the command line:
   "geneTree = hclust(as.dist(dissTOM), method = "average");"
   Click Run.
4. Use the cutreeDynamic function to perform dynamic branch cutting and set the minClusterSize parameter to 30. Obtain the module recognition result. Execute the command line:
   "dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize);"
   Click Run.
   NOTE: The minimum module size could not be lower than 30.
5. Calculate the module eigen of each OTUs module by the moduleEigengenes function. Execute the command line:
   "MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
   MEs = MEList$eigengenes"
   Click Run.
   NOTE: The module eigen represented the overall OTU expression level in the module. It was not a specific OTU, but the first principal component of each cluster obtained by singular network value decomposition.
6. Perform the cluster function based on the correlation coefficient of module eigen. Use the mergeCloseModules function to merge the modules with a value lower than 0.25. Execute the command line:
   "merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)"
   Click Run.
7. Finally, use the plotDendroAndColors function for visualization to obtain the module assignment display diagram of each co-expression network (Figure 2). Use the table function to extract the module attribution corresponding of each OTin the module assignment table. Execute the command line:
   "plotDendroAndColors(geneTree, mergedColors, "Merged dynamic",dendroLabels = FALSE,
      hang = 0.03,addGuide = TRUE, guideHang = 0.05,
      main = "ES_Gene dendrogram and module colors")"
   Click Run.
   ​NOTE: In the module assignment diagram of the co-expressing network, different colors represent different modules, and gray represents OTUs that cannot be classified into any module. A greater number of OTUs in the gray module indicates that the early-stage preprocessing quality of the expression matrix is poor.

4. Module comparison

NOTE: This method can be used to compare the network modules of two ecological microbial communities. In this article, compare the differences of microbial network modules between endosphere and rhizoplane, endosphere and rhizosphere, rhizosphere and rhizoplane.

1. Preservation test
   1. Load the parameters and results of the two data sets saved in the previous steps.
   2. Set the network module assignment result of a group of microbial data as the reference group, whereas the other group as the test group.
   3. Use the modulePreservation function to calculate the values of conservativeness statistical parameters Z_summary and medianRank. Execute the command line:
      "system.time({mp=modulePreservation(multiExpr,
      multiColor,referenceNetworks=1,
      nPermutation=100, randomSeed=1,quickCor=0,verbose=3)})"
      Click Run.
      NOTE: This result can quantified the conservativeness between modules. Z_summary>10 indicates that two modules are highly preserved, whereas Z_summary<2 denotes non- preserved modules. medianRank expresses the relative preservation of the module assessed by ranking. Higher medianRank values denote non- preserved modules. (The key command lines refer to the Supplement S3.)
   4. Use the plot function to visualize the results (Figure 3). Get the parameters Z_summary and medianRank (Table 1).
      NOTE: The network modules that satisfy both the Z_summary value less than 2 and the median Rank value at the top, is the most highly non-preserved module in the two ecological microbial communities.
   5. Based on the results of the aforementioned two statistical parameters to identify the module with most highly non-preserved module of the two networks.
2. Correlation analysis of the module membership
   1. Set the module assignment results of the two networks were set as the reference and the test group, respectively.
      NOTE: The settings need to be the same as Preservation test.
   2. Use the corPvalueStudent function to extract the kME (module membership) value of each OTU in several candidate modules.
      Execute the command line:
      "Pvalue = as.data.frame(corPvalueStudent(as.matrix
      (ModuleMembership), Samples))"
      Click Run.
      NOTE: kME stands for the degree of module membership. ME stands for module eigen, which represents the overall level of OTU expression in the module. kME is the correlation coefficient between each OTU and the ME. Quantify the importance of OTU in the network by the kME value of OTU. (The key command lines refer to the Supplement S4.)
   3. Then, use the verboseScatterplot function to calculate the correlation coefficient of the kME value of the corresponding OTUs in the two networks and draw the correlation analysis diagram (Figure 4).
      Execute the command line:
      "verboseScatterplot(abs(TModuleMembership
      [TmoduleGenes, Tcolumn]),
         abs(NModuleMembership[NmoduleGenes, Ncolumn]),
         xlab = paste("kME in", "ES"),
         ylab = paste("kME in", "RP"),
         main = paste("lightyellow"),
         cex.main = 1.7, cex.lab = 1.6, cex.axis = 1.6, col =          modulecolor)"
      Click Run.
   4. Select the module with the smallest correlation coefficient of the kME value of the OTU of the two networks. Consider this module to have the largest difference of the two networks.

5. Analysis of the microbial differential network module

1. Obtain data of the dominant bacteria phyla through statistical analysis of the OTU sequence set of the module with the largest difference.
   NOTE: The OTU sequence set of the module with the largest difference is summed by the taxonomy of phyla. The dominant bacteria phyla accounted for more than 10%.
2. Then, use the exportNetworkToCytoscape function to obtain the file containing the interaction relationship information of the OTU in the largest differential module.
   Execute the command line:
   "cyt = exportNetworkToCytoscape(modTOM,
   edgeFile = paste("NEW-ES_CytoscapeInput-edges-", modules , ".txt", sep=""),
   nodeFile = paste("NEW-ES_CytoscapeInput-nodes-", modules, ".txt", sep=""),
   weighted = TRUE,threshold = 0.5, nodeNames = modProbes,
   altNodeNames = modGenes, nodeAttr = moduleColors[inModule])"
   Click Run.
3. Import the file into Cytoscape. Set the threshold to 0.5 and adjust other parameters as needed.
4. Construct a co-occurrence network of differential microorganisms (Figure 5).
5. Obtained the information of the core genus that has the most important regulatory role in the network.
   NOTE: According to the the kME value of OUT, the core genus can be defined.
6. Finally, the functions of the core genus were assessed and its influence on the entire difference network was analyzed.

Representative Results

The representative results in this article were downloaded from the 2014 California Abaker rice root microbiome data in the NCBI database (PRJNA386367)⁹. The data includes the rhizosphere, rhizoplane, and endosphere microbiome samples from rice plants grown for 14 weeks in a submerged rice field. We used the WGCNA algorithm to select the power value that satisfied the three networks that were close to the scale-free network (Figure 1) and developed three co-expression networks (Figure 2). In the endosphere, rhizoplane, and rhizosphere soil microbial co-expression network, 23, 22, and 21 modules were identified, respectively. These results indicate that the numbers of microbial interaction networks in the three niches were basically equal.

We further compared the differences in the microbial network modules within the endosphere, rhizoplane, and rhizosphere soil. The following preservation test results of the modules of the three niches groups were obtained. Three extremely non-preserved modules existed between the rhizosphere soil and the rhizoplane (Figure 3a, Table 1). Additionally, nine extremely non-preserved modules were present between the rhizosphere soil and endosphere (Figure 3b, Table 2) and six extremely non-preserved modules between the rhizoplane and endosphere (Figure 3c, Table 3). Furthermore, extremely non-preserved modules among the three niches were found, indicating the presence of large differences in the composition of microorganisms among the three niches. The results of the module membership correlation analysis of the obtained non-conservative modules are illustrated in Figure 4. From the figure, a significantly different module with the least correlation between each two niches is visible among the three niches, which represents the most significant difference between them.

In the rhizosphere-rhizoplane difference network (Figure 5a), the dominant phylum was Proteobacteria (72.97%). In the rhizosphere-endosphere difference network (Figure 5b), the dominant phyla were Proteobacteria (66.36%), Actinobacteria (10.1%), and Bacteroidetes (10.9%). In the rhizoplane-endosphere difference network (Figure 5c), the dominant phyla were Proteobacteria (41.41%), Bacteroidetes (10.10%), Firmicutes (12.12%), and Verrucomicrobia.

Three core genera (Figure 5a), including Rhodobacter and Novosphingobium, six core genera (Figure 5b), including Blvii28 and Dechloromonas, and five core genera (Figure 5c), including Cellvibrio and Geobacter, exerted important regulatory functions in the three differential co-occurrence networks. All core genera, except for Dechloromonas, had influence on only one network, indicating the availability of considerable differences in the relative abundance of microbial populations and species among the three niches of rice roots, which critically affected the abundance and diversity of the existing root microbial communities.

The core genus Azospirillum, present in the rhizosphere-endosphere difference network of rice, participated in nitrogen fixation and promoted plant growth¹⁰. Additionally, the Geobacter genus, which was significantly enriched in the rhizoplane-endosphere difference network, may be the main factor inducing the reduction of insoluble Fe and Mn oxides in many soils and sediments¹¹. These microorganisms interact with a range of microbial communities in the root niches and actively participate in the regulation of microbial networks, which might be critically important for the growth and development of rice roots.

Figure 1. Evaluation of powerβ in the datasets. (a) Evaluation of power β in the ES dataset; (b) Evaluation of power β in the RS dataset; (c) Evaluation of power β in the RP dataset, distribution of the scale-free index R2 (left), distribution of the mean connectivity (right) along different soft power indices. The value of the best power was achieved when R2 tended to saturation and was not lower than 0.8. The three co-expression networks had to be set to the same power value to ensure their comparability. (RS: rhizosphere soil, RP: rhizoplane, and ES: endosphere). Please click here to view a larger version of this figure.

Figure 2. OTU dendrogram obtained by average linkage hierarchical clustering. (a) OTU dendrogram from ES network; (b) OTU dendrogram from RS network; (c) OTU dendrogram from RP network. The color row underneath the dendrogram indicates the module assignment determined by the Dynamic Tree Cut algorithm. (RS: rhizosphere soil, RP: rhizoplane, and ES: endosphere) Please click here to view a larger version of this figure.

Figure 3. Results of the preservation test. (a) The analysis results are based on the RS co-expression network module allocation as the reference group and the RP co-expression network module allocation results as the test group; (b) The analysis results are based on the ES co-expression network module allocation as the reference group and the RS co-expression network module allocation results as the test group; (c) The analysis results are based on the ES co-expression network module allocation as the reference group and the RP co-expression network module allocation results as the test group. Z_summary > 10 indicates that two modules are highly preserved, whereas Z_summary < 2 denotes non- preserved modules. medianRank expresses the relative preservation of the module assessed by ranking. Higher medianRank values denote non- preserved modules. (RS: rhizosphere soil, RP: rhizoplane, and ES: endosphere) Please click here to view a larger version of this figure.

Figure 4. Correlation analysis of the module memberships. (a) Module correlation of the kME value between the RS and RP network; (b) Module correlation of the kME value between RS and RP network; (c) Module correlation of the kME value between RP and ES network (RS: rhizosphere soil, RP: rhizoplane, ES: endosphere, and KME: module membership). Please click here to view a larger version of this figure.

Figure 5. Co-occurrence network of differential microbialpopulation in the root of rice. (a) Co-occurrence network of differential microbial populations in RS-RP; (b) Co-occurrence network of differential microbial populations in ES-RS; (c) Co-occurrence network of differential microbial populations in ES-RP. The analysis was performed using Cytoscape software. Different colors represent different gates in the figure. (RS: Rhizosphere soil, RP: rhizoplane, and ES: endosphere). Please click here to view a larger version of this figure.

module	medianRank	Zsummary	module	medianRank	Zsummary
greenyellow	2	14	yellow	13	5.2
pink	2	17	grey60	14	3.1
midnightblue	3	10	royalblue	15	2.1
brown	4	18	black	16	2.7
lightcyan	6	8.5	tan	17	3.2
purple	7	10	salmon	18	1.3
blue	7	22	magenta	18	2.2
green	8	12	darkred	20	-0.24
cyan	10	4.8	gold	20	14
lightgreen	11	5.1	darkgreen	22	-1.1
lightyellow	12	5.3	grey	22	0.21
red	12	6.1

Table 1. Result of Zsummary and medianRank between the rhizosphere soil and the rhizoplane.

module	medianRank	Zsummary	module	medianRank	Zsummary
salmon	1	19	lightcyan	13	1.1
black	3	5.3	red	15	0.95
yellow	4	5.8	midnightblue	15	-0.0016
lightgreen	5	0.27	royalblue	16	0.83
greenyellow	7	3	magenta	16	0.52
darkturquoise	7	1.2	darkgreen	17	0.16
grey60	9	1.1	tan	18	0.64
blue	10	3.9	lightyellow	19	0.52
purple	10	2.3	darkred	19	-0.18
brown	12	2.3	pink	19	-0.71
cyan	12	0.78	gold	21	11
green	13	1.7

Table 2. Result of Zsummary and medianRank between the rhizosphere soil and endosphere.

module	medianRank	Zsummary	module	medianRank	Zsummary
black	1	15	darkturquoise	13	1.7
salmon	2	27	midnightblue	13	1.6
yellow	3	13	lightgreen	13	0.64
cyan	4	5.4	darkgreen	14	1.5
blue	8	3.9	darkred	16	1.5
lightcyan	9	2.6	purple	17	2.3
pink	10	3.6	greenyellow	18	0.8
royalblue	10	1.5	lightyellow	18	0.42
brown	12	2.9	magenta	19	0.2
green	12	1.9	gold	21	18
red	12	1.9	grey60	21	-0.21
tan	13	2.5

Table 3. Result of Zsummary and medianRank between the rhizoplane and endosphere.

Supplement S1: Please click here to download this File.

Supplement S2: Please click here to download this File.

Supplement S3: Please click here to download this File.

Supplement S4: Please click here to download this File.

Discussion

Correlation networks have been increasingly used in bioinformatics applications. WGCNA is a systems biology method for descriptive analysis of the relationships between various elements of a biological system¹². R software package was used in earlier work on WGCNA^13,14,15. The package includes functions for network construction, module detection, calculations of topological properties, data simulation, visualization, and capability for interfacing with external software. WGCNA has been extensively employed to analyze gene expression data from brain cancer¹⁶, yeast cell cycle¹⁷, mouse genetics^18,19, primate brain tissue^20,21, diabetes²², and plants²³. Use the weighted gene correlation network analysis to construct the network must include at least 8 sample. In this article, we focused on gene co-expression networks that describe the interactions between microbial populations in different environments. We obtained differential networks between microbial populations in divergent environments and identified the key species in each network. The notion that key species are important to community has been widely employed in food-web research²⁴. Some of the species in a complex microbial community can be essential to maintain the stability and functionality of the community, such as Bacteroides in the intestinal flora²⁵. The analysis of microbial communities can be substantially simplified by targeting specific species of potential importance.

Our representative results highlight the differences in microbial communities, which can be identified using the above-described method. Here, microorganisms in different niches of the rice root system were subjected to WGCNA. The differential among the three niches were identified using conservative and module membership analyses. We determined the key species in the difference modules and obtained information on the differences in the composition of the microbial communities in the three niches. Meanwhile, the co-occurrence network revealed the presence of a significant interaction between the changing microorganisms in the roots of rice. Our findings provide direct evidence for the significance and feasibility of WGCNA in the evaluation of microbial community differences in various environments.

Materials

Name	Company	Catalog Number	Comments
R	The University of Auckland	version 4.0.2	R is a free software environment for statistical computing and graphics. It compiles and runs on a wide variety of UNIX platforms, Windows and MacOS.
RStdio	JJ Allaire	version 1.4.1103	The RStudio IDE is a set of integrated tools designed to help you be more productive with R and Python.
Cytoscape		version 3.7.1	Cytoscape is an open source software platform for visualizing complex networks and integrating these with any type of attribute data.
NCBI database			The National Center for Biotechnology Information advances science and health by providing access to biomedical and genomic information.