This article is mentioned in:

Introduction

With improvements in living standards, increasing attention is focused on intestinal health. Intestinal health is defined as a state, in which there are no signs of gut disease, no risk of intestinal disease and no intestinal diseases (1). The adhesion of intestinal flora to the intestinal wall can have a beneficial or deleterious role in hosts, and determines the health status of the gut (2). Therefore, investigations of intestinal flora are essential for intestinal health. A previous study showed that the diversity of the colonic mucosa was closely associated with bacterial microflora in patients with active inflammatory bowel disease (3). Disturbances of microbiota may also induce immune dysregulation and can develop into inflammatory bowel disease (4). In addition, gut bacteria can degrade various xenobiotics, including halogenated aromatic compounds, which may be beneficial for human health (5).

Various genes associated with intestinal flora have been identified for investigating intestinal health. In the microbial community, the diversity of genes can provide several types of enzymes and pathways associated with intestinal health (6). Jia et al (7) suggested that cytochrome P450, family 4, subfamily F, polypeptide 14 and tachykinin, precursor 1 interact with each other and are involved in pathways, including the neuropeptide signaling pathway, oxidation reduction and metabolism by regulating intestinal microbiota depletion. In addition, FBJ murine osteosarcoma viral oncogene homolog has been confirmed to enhance gut health by altering the ecology of the gut microbiota and improving the proteolysis of feces (8). Another important gene, interleukin 6 has been confirmed as a key molecule in gut barrier dysfunction (9). Although several genes associated with intestinal health have been investigated, further information is required. Therefore, molecular mechanisms, including critical genes, pathways and their interactions, require investigation.

In the present study, microarray analysis was performed for the screening of differentially expressed genes (DEGs) between microbiota-depleted mice and control mice. Subsequently, functional and pathway enrichment analyses of the DEGs were processed. Finally, literature associated with the DEGs was mined and their associations were analyzed. The results of the present study may identify potential important genes for further investigations on intestinal flora and human health.

Materials and methods

Data acquisition

The gene expression profiles of GSE22648 were downloaded from the Gene Expression Omnibus database (www.ncbi.nlm.nih.gov/geo) with the platform GPL6887 (Illumina MouseWG-6 v2.0 expression beadchip) (10). A total of 11 colon intestinal epithelial cell samples were obtained from mice, including five microbiota-depleted mice and six control mice, and analyzed.

Data preprocessing and DEG screening

The Series Matrix File (ftp.ncbi.nlm.nih.gov/geo/series/GSE22nnn/GSE22648/matrix) was downloaded to transfer probe names into gene symbols based on the platform annotation information. Using the aggregate function of the R statistical package, version 3.4.0 (stat.ethz.ch/R-manual/R-devel/library/stats/html/aggregate.html), the mean value was considered to be the expression value of a gene when multiple probes mapped to a single gene. For probes with missing values, nearest neighbor averaging (11) in the impure package (12) of R was used for offset of the margin value (k value was defaulted to 10). Quantile normalization was applied for standardization based on the preprocessCore package of R. The standard matrix was then obtained.

The DEGs between two groups were screened using the Limma package (13) (bioconductor.org/packages/release/bioc/html/limma.html; version 3.5) and the differences of mean expression values were assessed using the unpaired t-test method. In addition, the Benjamini-Hochberg algorithm (14) was used for p-value correction. Finally, the thresholds of the DEGs were log2 fold change |log2FC| >0.585 and adjusted p-value of P<0.05.

Functional and pathway enrichment of DEGs

The online tool of the Database of Annotation, Visualization and Integrated Discover (DAVID) (15) was used for Gene Ontology (GO; www.geneontology.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.jp/kegg) pathway enrichment analysis of DEGs. The cut-off criterion was P<0.05.

Mining of DEG-associated literature

The DEG-associated literature was mined using the GenCLip 2.0 online tool (ci.smu.edu.cn) (16). The Gene Cluster with Literature Profiles module generated statistically overrepresented key words grouped by a fuzzy cluster algorithm to annotate the input genes. The key words were generated based on the occurrence frequencies of free terms in the gene-associated literature or were provided by the user. The associations among the genes and keywords were linked to the relevant MEDLINE abstracts in which the co-occurrence of genes and keywords are highlighted. The Literature Mining Gene Networks module was used to construct a gene-network of the input genes and generate sub-networks based on the user defined query terms. It also calculated the probability of the random occurrence of the networks through random simulation, and provided GO and pathway enrichment analyses of genes. Co-citation networks of DEGs were constructed for the mining of ‘hot genes’ which were observed to be repeatedly associated with the other genes. At the same time, the DEGs were clustered to mine the functional information of key genes in previous studies with thresholds of P<0.0001 and hit ≥4.

Results

Data preprocessing and DEG screening

Following data preprocessing, a total of 30,865 expression data of transcripts were obtained from 11 samples. As shown in Fig. 1, the median of each sample was at the same level.

Figure 1. Box-plot of raw data standardization. Data are shown (A) prior to and (B) following normalization. Sample names are on the horizontal axis and the logarithm transformations of the expression values are shown on the vertical axis.

A total of 115 DEGs were screened with the thresholds of |log2FC| >0.585 and adjusted p-value of P<0.05, which included 58 upregulated genes and 57 downregulated genes. The cluster heatmap of the DEGs is shown in Fig. 2.

Figure 2. Cluster heatmap of differently expressed genes and samples.

GO functional and KEGG pathway enrichment analysis

As shown in Table I, the screened upregulated genes were significantly enriched into various GO terms, including microsome, oxidation reduction and heme binding. They were also enriched in the pathway of drug metabolism. The 57 downregulated DEGs were enriched in different functions, including DNA packaging, nucleosome assembly and nucleosome. These downregulated DEGs were also enriched in pathways, including arachidonic acid metabolism, systemic lupus erythematosus and linoleic acid metabolism.

Table I. Top five GO and pathway enrichment terms of differentially expressed genes.

Table I.

Top five GO and pathway enrichment terms of differentially expressed genes.

Ontology	Term	Count	P-value
Upregulated
GO-CC	GO:0005792~microsome	6	5.59E-05
GO-CC	GO:0042598~vesicular fraction	6	6.56E-05
GO-CC	GO:0005624~membrane fraction	8	1.62E-04
GO-CC	GO:0005626~insoluble fraction	8	2.01E-04
GO-CC	GO:0000267~cell fraction	8	4.21E-04
GO-BP	GO:0055114~oxidation reduction	8	7.98E-04
GO-MF	GO:0020037~heme binding	7	8.59E-07
GO-MF	GO:0046906~tetrapyrrole binding	7	1.14E-06
GO-MF	GO:0009055~electron carrier activity	7	6.14E-06
GO-MF	GO:0070330~aromatase activity	4	6.55E-05
GO-MF	GO:0005506~iron ion binding	7	1.21E-04
KEGG_PATHWAY	mmu00982:Drug metabolism	4	0.00128
Downregulated
GO-BP	GO:0006323~DNA packaging	6	7.41E-06
GO-BP	GO:0006334~nucleosome assembly	5	4.43E-05
GO-BP	GO:0031497~chromatin assembly	5	4.93E-05
GO-BP	GO:0034728~nucleosome organization	5	5.19E-05
GO-BP	GO:0065004~protein-DNA complex assembly	5	5.19E-05
GO-CC	GO:0000786~nucleosome	5	3.96E-05
GO-CC	GO:0044427~chromosomal part	8	4.92E-05
GO-CC	GO:0032993~protein-DNA complex	5	8.38E-05
GO-CC	GO:0005694~chromosome	8	1.46E-04
GO-CC	GO:0000785~chromatin	6	1.51E-04
KEGG_PATHWAY	mmu00590:Arachidonic acid metabolism	4	0.00171
KEGG_PATHWAY	mmu05322:Systemic lupus erythematosus	4	0.00318
KEGG_PATHWAY	mmu00591:Linoleic acid metabolism	3	0.00792

[i] GO, Gene Ontology; CC, cellular component; BP, biological process; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

DEG-associated literature mining

By mining the citation associations between the DEGs in PubMed, at least eight genes were found to have pairwise co-citation associations. In addition, four upregulated DEGs and four downregulated DEGs were identified with high co-citation levels (Fig. 3 and Table II). Therefore, four downregulated DEGs, including baculoviral IAP repeat containing 5 (BIRC5), aurora kinase A (AURKA), aurora kinase B (AURKB) and inner centromere protein antigens 135/155 kDa (INCENP) were co-cited in several studies. The four upregulated genes were laminin α3, bone morphogenetic protein 1, solute carrier family 22 (organic cation/zwitterion transporter), member 4 (SLC22A4) and PDZ domain containing 3.

Figure 3. Co-citation network of differently expressed genes. BIRC5, baculoviral IAP repeat containing 5; AURK, aurora kinase; ICENP, inner centromere protein antigens; LAMA3, laminin α3; BMP1, bone morphogenetic protein 1; SLC22A4, solute carrier family 22 (organic cation/zwitterion transporter), member 4; PDZ3, PDZ domain containing 3. The lines represent the co-citation relationship of two genes.

Table II. Results of literature-associated mining of differentially expressed genes.

Table II.

Results of literature-associated mining of differentially expressed genes.

Gene	Co-genes (n)	Co-cites (n)	Total (n)	log2 fold change
BIRC5	3	143	5,645	−0.971
AURKA	3	114	990	−1.054
AURKB	3	169	537	−0.700
INCENP	3	127	247	−0.710
LAMA3	1	1	301	0.896
BMP1	1	1	445	0.663
SLC22A4	1	2	271	0.790
PDZD3	1	2	27	0.715

[i] BIRC5, baculoviral IAP repeat containing 5; AURK, aurora kinase; ICENP, inner centromere protein antigens; LAMA3, laminin α3; BMP1, bone morphogenetic protein 1; SLC22A4, solute carrier family 22 (organic cation/zwitterion transporter) member 4; PDZ3, PDZ domain containing 3.

Gene cluster analysis with the literature profiles of all DEGs was processed (Fig. 4 and Table III). The results showed that a total of 19 genes, including BIRC5, AURKB, angiotensin I converting enzyme 2 and free fatty acid receptor 2, were enriched in four modules, including cell division, chromosome segregation, inflammatory bowel disease and inflammatory response. In addition, these modules were associated with various functional terms, including inflammatory response, chromosome condensation, chromosomal passenger complex and phospholipase A2.

Figure 4. Cluster heatmap of differently expressed genes from literature mining.

Table III. Modules of literature-associated differentially expressed genes.

Table III.

Modules of literature-associated differentially expressed genes.

Key word	Hit	Q-value	Genes
#single1 Enrichment Score: 6.21
INFLAMMATORY RESPONSE	17	0.000223	ACE2; BIRC5; CD177; CD276; CHAC1; CYP4B1; F2R; FFAR2; GPX2; IRF8; MYLK; NFKBIZ; PLA2G2A; REG3G; RETNLB; SAA1; SAA2
#cluster1 Enrichment Score: 6.17
CHROMOSOME CONDENSATION	6	0.000023	AURKB; BIRC5; H1F0; HIST2H2BE; INCENP; NCAPH
CHROMOSOME SEGREGATION	7	0.000622	AURKA; AURKB; BIRC5; GSG2; HIST2H2BE; INCENP; NCAPH
CELL DIVISION	14	0.000568	ACE2; AURKA; AURKB; BIRC5; CD177; FBXO6; FNDC3B; GSG2; H1F0; HIST2H2BE; INCENP; MYLK; NCAPH; SIDT1
#cluster2 Enrichment Score: 5.66
CHROMOSOMAL PASSENGER COMPLEX	5	0.000016	AURKB; BIRC5; GSG2; INCENP; NBL1
SPINDLE ASSEMBLY CHECKPOINT	5	0.002422	AURKA; AURKB; BIRC5; GSG2; INCENP
SPINDLE CHECKPOINT	5	0.003230	AURKA; AURKB; BIRC5; GSG2; INCENP
#single2 Enrichment Score: 5.42
INFLAMMATORY BOWEL DISEASE	9	0.000590	AGR2; FFAR2; GPX2; MYLK; NFKBIZ; REG3G; RETNLB; SAA1; SLC22A4
#single3 Enrichment Score: 4.73
PHOSPHOLIPASE A2	6	0.002282	CD177; MYLK; PLA2G2A; PLA2G3; SAA1; SMPDL3B

Discussion

Gut bacteria and their reactions are important in gut health and general human health (17). Understanding the process and potential mechanism of the effect of microbial intestinal metabolism on human health is important for the treatment of intestinal diseases (18). In the present study, four upregulated DEGs and four downregulated DEGs were identified with high co-citation levels. Among these genes, the AURKA, INCENP, BIRC5 and AURKB downregulated genes and SLC22A4 upregulated gene were important DEGs, which were enriched in four modules, including cell division, chromosome segregation, inflammatory bowel disease and inflammatory response.

A previous study indicated that intestinal cell proliferation is associated with epithelial homeostasis, regeneration (19), and various diseases, including obstructive jaundice (20) and cancer (4). In the present study, AURKA was identified with a high number of citations on literature mining. This gene encodes a protein, which is a cell cycle regulated kinase, and is also closely associated with chromosome segregation (21). Ratushny et al (22) suggested that by combining with AURKA, SRC kinase was activated and further affected various biological processes, including tumor invasion. In addition, AURKA locates at the spindle poles and regulates cytokinesis by regulating centrosome maturation (23). In the colon of HNF4α mutant mice, the expression of AURKA was significantly upregulated, and HNF4α was a key regulator for ion transport in the gut mucosa. Therefore, the findings of the present study suggested that AURKA may be involved in the functional modules of chromosome segregation and cell division, and may be important in affecting intestinal microbiome homeostasis in intestinal health.

AURKA is concentrated primarily at the spindle poles in the period of meiosis (MI) and meiosis II, whereas AURKB locates at kinetochores in the period of MI (24). As with AURKA, AURKB was identified in the present study as being enriched in functions of chromosome segregation and cell division. Fenton et al (25) processed a microarray analysis, and showed that leptin may upregulate the expression of AURKB and further induce colon epithelial cell proliferation. In addition, the expression of AURKB can be regulated by microRNA (miR)-24 by recognizing complementary sequences (26). In addition, in irritable bowel syndrome, miR-24 has been confirmed to increase gut permeability (27). Therefore, AURKB may also be a target in the treatment of intestinal disease.

INCENP is a complex partner of AURKB for kinase activation, and the regulation of chromosome segregation and cytokinesis (28). In the present study, INCENP was screened to have a high number of citations, and was enriched in cell division and chromosome segregation. INCENP has also been confirmed to be an important member of the chromosomal passenger complex, which is an essential regulator for the completion of cell division (29). In addition, INCENP/Sli15 phosphorylation can activate Aurora B/IpI1, which is important for the activity of the chromosomal passenger complex (30). The expression of INCENP has also been shown to be increased in several types of cancer, including colorectal cancer (31). Therefore, the present study hypothesized that INCENP may affect intestinal health through interacting with AURKB.

In the present study, BIRC5 was enriched in three modules, including cell division, chromosome segregation and inflammatory response. Neufert et al (32) confirmed that BIRC5 can limit bacterial growth and is involved in mucosal wound healing, which suggests a promising therapeutic approach for disturbed intestinal homeostasis. In addition, activated signal transducer and activator of transcription 3 can induce the transcription of BIRC5 and regulate tissue defense during the process of gastrointestinal infection. Previous studies have also shown that microbiota and commensal bacteria can induce inflammatory response in immunocompromised hosts (33), and when the normal microbial community is disrupted, the risk of inflammation is increased (34), This offers novel insights into the pathogenesis of disease for treating intestinal disorders. The gut microbiota can also affect host metabolism and produce various pro-inflammatory molecules, including lipopolysaccharide, assisting in the development of treatments for metabolic diseases (35). Therefore, BIRC5 may be a critical gene in the effects of microbiota on intestinal health.

Of the upregulated DEGs identified in the present study, SLC22A4 was found to have a high number of citations and was enriched in the module of inflammatory bowel disease. A mutation of SLC22A4, which encodes organic carnitine transporter 1, has been confirmed to have a clinical role in gut health (36). In addition, Tokuhiro et al (37) showed that the organic cation transporter SLC22A4 functioned as a transporter in the inflammatory milieu. A previous study also confirmed that polymorphisms of SLC22A4 were closely associated with early inflammatory bowel disease (38). The development of the gut microbiota is essential for promoting the gut immune system and maintaining system balance (39). Therefore, SLC22A4 may be an essential gene for intestinal health via modulating the immune system.

In conclusion, the present study performed bioinformatics analysis of the effects of gut flora on molecular changes in intestinal epithelial cells. AURKA, INCENP, BIRC5, AURKB and SLC22A4 were identified as potential target genes for intestinal disease treatment through their involvement in various functions, including cell division, chromosome segregation, inflammatory bowel disease and inflammatory response. These results indicate novel therapeutic targets for diseases caused by disturbed intestinal homeostasis. However, there were limitations to the present study; gut microbiota are associated with the delivery mode and the type of feeding in infants, and changes throughout life between infancy and elderly. Therefore, further experiments are required to confirm the results obtained in the present study.

References