Gene co-expression network analysis of the human gut commensal bacterium Faecalibacterium prausnitzii in R-Shiny

Abstract

Faecalibacterium prausnitzii is abundant in the healthy human intestinal microbiota, and the absence or scarcity of this bacterium has been linked with inflammatory diseases and metabolic disorders. F. prausnitzii thus shows promise as a next-generation probiotic for use in restoring the balance of the gut microbial flora and, due to its strong anti-inflammatory properties, for the treatment of certain pathological conditions. However, very little information is available about gene function and regulation in this species. Here, we utilized a systems biology approach-weighted gene co-expression network analysis (WGCNA)-to analyze gene expression in three publicly available RNAseq datasets from F. prausnitzii strain A2-165, all obtained in different laboratory conditions. The co-expression network was then subdivided into 24 co-expression gene modules. A subsequent enrichment analysis revealed that these modules are associated with different kinds of biological processes, such as arginine, histidine, cobalamin, or fatty acid metabolism as well as bacteriophage function, molecular chaperones, stress response, or SOS response. Some genes appeared to be associated with mechanisms of protection against oxidative stress and could be essential for F. prausnitzii's adaptation and survival under anaerobic laboratory conditions. Hub and bottleneck genes were identified by analyses of intramodular connectivity and betweenness, respectively; this highlighted the high connectivity of genes located on mobile genetic elements, which could promote the genetic evolution of F. prausnitzii within its ecological niche. This study provides the first exploration of the complex regulatory networks in F. prausnitzii, and all of the "omics" data are available online for exploration through a graphical interface at https://shiny.migale.inrae.fr/app/faeprau.

Fig 1. Partial least-square discriminant analysis of RNA-seq data.

(A) PLS-DA plot of RNA-seq data showing clear transcriptome-based discrimination among the three publication datasets. Each point represents the transcriptome signature of one sample, with ellipses representing 95% confidence level. (B) Genes with high discriminatory ability were identified from PLS-DA. (C) Histogram showing the distribution of COG categories (Clusters of Orthologous Genes) associated with the genes with the highest discriminatory power in the two components generated by PLS-DA. See S2 Table for the full list of discriminatory genes.

Fig 2. Illustration of transcriptomic data mining.

(A) and (C) Visualization of RNAseq results with volcano plots. (B, D, E, F, G, H) Boxplots showing the read counts for genes of interest among different growth conditions.

Fig 3. Network construction and module detection with WGCNA.

(A) Network topology analysis for various values of soft-thresholding powers, with the scale-free index and the mean connectivity as a function of the soft-thresholding power. (B) Dendrogram of all genes divided into 24 modules, with dissimilarity based on topological overlap, presented with assigned module colors. (C) Dendogram representing the 24 modules identified by WGCNA. The heightcut (red line, heightcut = 0.0) was used to unmerge modules. The grey module represents genes that are not included in any of the other modules. (D) A heatmap depicting the topological overlap matrix (TOM) among all genes in the analysis. The intensity of the red color indicates the strength of the correlation between all pairs of genes.

Fig 4. Identification of the modules associated with conditions in the three original datasets.

(A) Heatmap depicting the correlation between module eigengenes and the original datasets. Pearson coefficient correlations are indicated. The p-value is indicated in parentheses. (B) Heatmap of gene expression levels in the modules across the samples in the three original datasets.

Fig 5. Illustration of enrichment analysis with STRING (version 9.0).

Results are presented for the modules “Lightgreen” (red color: Cobalamin biosynthetic pathway), “Darkred” (blue color: Bacteriophage functions), “Royalblue” (red color: Fatty acid biosynthesis), “Purple” (green color: Transposon-encoded proteins; red color: Molybdenum cofactor biosynthesis; yellow color: Endonuclease/relaxase).

Fig 6. Network plot depicting the top connections in the “turquoise” module.

Nodes represent genes, and node size is correlated with the degree of connectivity of the gene.