Identification of key genes implicated in the suppressive function of human FOXP3+CD25+CD4+ regulatory T cells through the analysis of time‑series data.

Human forkhead box P3 (FOXP3)+ cluster of differentiation (CD)25+CD4+ regulatory T cells (Tregs) are a type of T cell that express CD4, CD25 and FOXP3, which are critical for maintaining immune homeostasis. The present study aimed to determine the mechanisms underlying Treg function. The GSE11292 dataset was downloaded from the Gene Expression Omnibus, which included data from Treg cells at 19 time points (0‑360 min) with an equal interval of 20 min, and corresponding repeated samples. However, data for Treg cells at time point 120 min were missing. Using the Mfuzz package, the key genes were identified by clustering analysis. Subsequently, regulatory networks and protein‑protein interaction (PPI) networks were constructed and merged into integrated networks using Cytoscape software. Using Database for Annotation, Visualization and Integrated Discover software, enrichment analyses were performed for the genes involved in the PPI networks. Cluster 1 (including 292 genes), cluster 2 (including 111 genes), cluster 3 (including 194 genes) and cluster 4 (including 103 genes) were obtained from the clustering analysis. GAPDH (degree, 40) in cluster 1, Janus kinase 2 (JAK2) (degree, 10) and signal transducer and activator of transcription 5A (STAT5A) (degree, 9) in cluster 3, and tumor necrosis factor (TNF) (degree, 26) and interleukin 2 (IL2) (degree, 22) in cluster 4 had higher degrees in the PPI networks. In addition, it was indicated that several genes may have a role in Treg function by targeting other genes [e.g. microRNA (miR)‑146b‑3p→TNF, miR‑146b‑5p→TNF, miR‑142‑5p→TNF and tripartite motif containing 28 (TRIM28)→GAPDH]. Enrichment analyses indicated that IL2 and TNF were enriched in the immune response and T cell receptor signaling pathway. In conclusion, GAPDH targeted by TRIM28, TNF targeted by miR‑146b‑3p, miR‑146b‑5p and miR‑142‑5p, in addition to JAK2, IL2, and STAT5A may serve important roles in Treg function.

Introduction

Regulatory T cells (Tregs) are a subgroup of T cells that suppress proliferation of effector T cells, sustain tolerance to self-antigens and regulate the immune system (1). As a type of T cell that expresses cluster of differentiation (CD)4, CD25 and forkhead box P3 (FOXP3), human FOXP3+CD25+CD4+ Tregs are critical for maintaining immune homeostasis (2). Tregs are deemed to inhibit tumor immunity and contribute to the growth of cancerous cells, suggesting that high levels of Tregs may indicate poor prognosis for patients with cancer (3). A previous study also demonstrated that regulation of Tregs is conducive to autoimmune disease and organ transplantation (4). Therefore, it is necessary to explore the mechanisms implicated in Treg function.

Numerous genes have been reported to be associated with Treg function. For example, cytotoxic T lymphocyte antigen 4 may be important for the immune suppression of natural Tregs by affecting the activation effects of antigen-presenting cells on other T cells (5–7). Ectopic expression of lymphocyte-activation gene 3 can significantly weaken the proliferative capacity of CD4+ T cells and facilitate their inhibitory effect on effector T cells (8). Prostaglandin E2 (PGE2) promotes the mRNA and protein expression of FOXP3 and increases its promoter activity, thus suggesting that PGE2 in human lymphocytes may regulate FOXP3 expression and the function of Tregs (9,10). In vivo, toll-like receptor 2 (TLR2) modulates Treg function, thus indicating that TLRs may control immune responses via Tregs (11,12). Furthermore, a previous study demonstrated that indoleamine 2,3-dioxygenase 1 predominantly controls the response of Tregs to inflammatory stimuli in the physiological environment (13). Latent transforming growth factor-β is expressed on activated Tregs, and may serve a role in mechanisms underlying infectious tolerance and Treg-mediated suppression (14). In 2012, He et al (15) performed high time-resolution genome-wide gene expression analysis to investigate the genes involved in human Tregs; the results demonstrated that plasminogen activator urokinase was essential for the suppressive function of Tregs. Nevertheless, the potential molecular mechanisms underlying Treg function remain unclear.

The present study, using data deposited by He et al (15), further identified the key genes implicated in Treg function. After searching the microRNA (miRNA)-mRNA pairs, transcription factor (TF)-mRNA pairs and protein-protein interaction (PPI) relationships, regulatory networks and PPI networks were constructed and merged into integrated networks. Finally, enrichment analyses were performed for the genes involved in the PPI networks to predict their possible functions.

Materials and methods

Microarray data

GSE11292 microarray data were downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database, based on the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array platform (Affymetrix, Inc., Santa Clara, CA, USA). To investigate the key genes important for human Treg function, data from Treg cells at 19 time points (0–360 min) with an equal interval of 20 min, and their corresponding repeated samples, were collected from GSE11292. Notably, data from Treg cells at the 120 min time point were missing. Briefly, the stable human Tregs isolated from peripheral blood were from the same batch of Tregs as used in a previous study by Probst-Kepper et al (16). Tregs were cultured in Iscove's modified Dulbecco's medium (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA) containing 10% fetal calf serum (Gibco; Thermo Fisher Scientific, Inc.), 50 µM µ-mercaptoethanol (Sigma-Aldrich; Merck KGaA, Darmstadt, Germany), 100 µg/ml streptomycin and 100 U/ml penicillin. Human Tregs were added with anti-CD3/anti-CD28-coated Dynabeads (Invitrogen; Thermo Fisher Scientific, Inc.) in a proportion of 1:1 and interleukin-2 (IL2; 100 U/ml; Novartis International AG, Basel, Switzerland) and were assigned into 1.5 ml microtubes (Eppendorf, Hamburg, Germany; 4×106 cells/tube) for each time point. The cells were stored at 4°C to settle the cells and beads, and then were cultured at 37°C (15). GSE11292 data used in the present study were downloaded from the public GEO database; therefore, ethical approval and patient consent were not required.

Data preprocessing and clustering analysis

Using the Affy package (17) in R, the raw data were preprocessed with background correction, quantile normalization, probe summarization and transformation from probe ID to gene symbol. Subsequently, soft clustering analysis was performed for the two group samples (Treg cells and their repeated samples) using the Mfuzz package (18,19). The parameters of minimum standard deviation and acore were set at 0.5 and 0.9, respectively.

Construction of regulatory networks

Combined with the validated miRNA-mRNA pairs in miRecords database (http://c1.accurascience.com/miRecords/) (20), and the predicted miRNA-mRNA pairs in miRanda (http://www.microrna.org) (21), MirTarget2 (http://mir db.org/miRDB) (22), PicTar (http://pictar.mdc-berlin.de) (23), PITA (https://genie.weizmann.ac.il/pubs/mir07/mir07_data.html) (24) and TargetScan (http://www.targetscan.org) (25) databases, the miRNA-mRNA relationships involving the genes in each cluster were predicted. P<0.05 and the involvement of at least 2 genes were used as the thresholds for screening significant miRNAs. Subsequently, the screened miRNA-mRNA pairs were visualized in miRNA-mRNA regulatory networks using Cytoscape software (http://www.cyto scape.org) (26).

Based on the transcriptional regulatory data in the ENCyclopedia of DNA Elements (ENCODE) project (http://www.genome.gov/Pages/Research/ENCODE) (27), the transcriptional regulatory relationships between the genes in each cluster were searched and identified. In addition, transcriptional regulatory networks were constructed using Cytoscape software (26).

PPI network construction and network integration

The STRING (28) database (http://string-db.org) was applied to perform a PPI analysis for the genes in each cluster. The PPI pairs with a combined score (required confidence) >0.4 were selected, after which, a PPI network was constructed using Cytoscape software (29). Nodes were considered proteins in the PPI network, whereas their degrees corresponded to the number of interactions associated with them. Nodes with higher degrees were considered hub nodes (30). Finally, the miRNA-mRNA regulatory network, transcriptional regulatory networks and PPI networks were integrated separately for the genes in each cluster.

Functional and pathway enrichment analyses

The Gene Ontology database (GO; http://www.geneontology.org/) classifies functions according to three terms: Molecular function, biological process and cellular component (31). The Kyoto Encyclopedia of Genes and Genomes database (KEGG; http://www.genome.jp/kegg/) contains information regarding biological systems from systemic functional, genomic and chemical aspects (32). Using Database for Annotation, Visualization and Integrated Discovery software (33), GO functional and KEGG pathway enrichment analyses were separately conducted for the genes involved in PPI networks. P<0.05 and the involvement of at least 2 genes were used as the thresholds for screening significant terms.

Results

Clustering analysis

After preprocessing, cluster 1 [including 292 genes; such as tripartite motif containing 28 (TRIM28) and GAPDH], cluster 2 (including 111 genes), cluster 3 (including 194 genes) and cluster 4 [including 103 genes; such as tumor necrosis factor (TNF)] were obtained from soft clustering analysis. Genes in cluster 1 were significantly downregulated after 200 min and were significantly upregulated after 300 min. Genes in cluster 2 were continually downregulated, whereas genes in cluster 3 were continually upregulated. Genes in cluster 4 were significantly upregulated prior to 80 min and expression flattened after that time point (Fig. 1).

Figure 1.

Four clusters (cluster 1, 2, 3 and 4) were obtained by soft clustering analysis. The red, blue and turquoise colors indicate the match degrees between changes of genes and the major changes of the clusters. Red, blue and turquoise represent high, moderate and low match degrees respectively.

Network construction and integration

Using the following thresholds: P<0.05 and targeting at least 2 genes, miRNAs targeting the genes [such as miRNA (miR)-146b-3p→TNF, miR-146b-5p→TNF and miR-142-5p→TNF] in each cluster were enriched (Table I). There were no miRNAs significantly enriched for the genes in cluster 2. Based on the ENCODE project, the transcriptional regulatory relationships between the genes in each cluster were searched and identified. In the transcriptional regulatory network for genes in cluster 1, GAPDH was targeted by TRIM28. However, no transcriptional regulatory relationships were found for the genes in cluster 3. There were 656, 40, 129 and 234 PPIs demonstrated in clusters 1, 2, 3 and 4. The top 10 nodes with the highest degrees in the PPI networks for each cluster are presented in Table II, including GAPDH (degree, 40) in cluster 1, Janus kinase 2 (JAK2; degree, 10) and signal transducer and activator of transcription 5A (STAT5A; degree, 9) in cluster 3, and TNF (degree, 26) and IL2 (degree, 22) in cluster 4. Finally, the regulatory and PPI networks were integrated separately for the genes in clusters 1 (Fig. 2), 2 (Fig. 3), 3 (Fig. 4) and 4 (Fig. 5).

Table I.

miRNAs predicted for the genes in each cluster.

Cluster	miRNA	Gene number	Gene symbol	P-value
1	hsa-miR-204	2	ARPC1B, CDC25B	3.22×10−2
	hsa-miR-338-3p	2	BAX, COX8A	4.29×10−2
	hsa-miR-338-5p	2	BAX, COX8A	4.29×10−2
	hsa-miR-137	2	BAX, COX8A	4.67×10−2
3	hsa-miR-28-5p	3	BCL6, JAK2, STAT5A	1.45×10−3
	hsa-miR-28-3p	3	BCL6, JAK2, STAT5A	2.26×10−3
	hsa-miR-339-3p	2	JAK2, RGS2	3.22×10−3
	hsa-miR-339-5p	2	JAK2, RGS2	1.11×10−2
	hsa-miR-127-3p	2	BCL6, PDGFA	3.40×10−2
4	hsa-miR-146b-3p	6	CXCL8, IL10, MYC, NEDD9, NFKBIA, TNF	2.84×10−4
	hsa-miR-146b-5p	6	CXCL8, IL10, MYC, NEDD9, NFKBIA, TNF	3.61×10−4
	hsa-miR-125a-5p	6	CD40LG, CD80, CD83, IL10, LIF, MYC	2.57×10−3
	hsa-miR-139-5p	2	MYC, PTGS2	8.51×10−3
	hsa-miR-125a-3p	5	CD40LG, GADD45B, IL10, LIF, MYC	9.01×10−3
	hsa-miR-1224-5p	2	GADD45A, TNF	1.12×10−2
	hsa-miR-19a	2	NR4A2, TNF	1.12×10−2
	hsa-miR-455-5p	2	IL10, PTGS2	1.12×10−2
	hsa-miR-671-5p	2	GADD45A, RELB	1.12×10−2
	hsa-miR-142-5p	3	AHR, IL10, TNF	2.16×10−2
	hsa-miR-630	2	GADD45A, YES1	2.50×10−2
	hsa-miR-34a	2	MYC, SIRT1	3.83×10−2

miRNA, microRNA.

Table II.

Top 10 nodes with higher degrees in the protein-protein interaction networks for each cluster.

Gene	Degree
Cluster 1
GAPDH	40
TSPO	37
ISG15	32
HSP90AB1	31
OAS1	24
OASL	22
MX1	21
PSMC3	19
RPL8	19
RPS3	18
Cluster 2
IFIT1	5
TUBA1A	5
IRF7	4
MAVS	4
TUBA4A	4
IFIT3	4
HDAC4	3
RGS19	3
RPS6KA1	3
FOS	3
Cluster 3
JAK2	10
STAT5A	9
WDR43	8
CTNNB1	8
BCL6	8
WDR36	8
TSLP	6
HSPD1	6
MAK16	6
CCR8	5
Cluster 4
MYC	27
TNF	26
IL2	22
ICAM1	21
IL4	18
CD40LG	17
IL8	17
RELB	16
IL10	15
NFKBIA	14

Figure 2.

Integrated network of the miRNA-mRNA regulatory network, transcriptional regulatory network and PPI network for the genes in cluster 1. Green nodes represent downregulated genes; triangles indicate transcription factors; blue diamonds indicate miRNAs; arrows represent transcriptional regulatory relationships; T-shaped lines indicate miRNA regulatory relationships; lines represent PPI relationships. miRNA, microRNA; PPI, protein-protein interaction.

Figure 3.

Integrated network of the transcriptional regulatory network and PPI network for the genes in cluster 2. Green nodes represent downregulated genes; triangles indicate transcription factors; arrows represent transcriptional regulatory relationships; lines represent PPI relationships. PPI, protein-protein interaction.

Figure 4.

Integrated network of the miRNA-mRNA regulatory network and PPI network for the genes in cluster 3. Red nodes represent upregulated genes; blue diamonds indicate miRNAs; T-shaped lines indicate miRNA regulatory relationships; lines represent PPI relationships. miRNA, microRNA; PPI, protein-protein interaction.

Figure 5.

Integrated network of the miRNA-mRNA regulatory network, transcriptional regulatory network and PPI network for the genes in cluster 4. Red nodes represent upregulated genes; triangles indicate transcription factors; blue diamonds indicate miRNAs; arrows represent transcriptional regulatory relationships; T-shaped lines indicate miRNA regulatory relationships; lines represent PPI relationships. miRNA, microRNA; PPI, protein-protein interaction.

GO functional and KEGG pathway enrichment analyses were conducted for the genes involved in the PPI networks. The top 10 functions enriched for the genes involved in the PPI networks are listed in Table III. Genes in the PPI networks were enriched in functions including negative regulation of protein metabolic process (cluster 1; P=2.39×10−9), defense response (cluster 2; P=2.84×10−3), response to organic substance (cluster 3; P=2.20×10−5), and immune response (cluster 4; P=3.43×10−8; which involved IL2 and TNF). The top 10 pathways enriched for the genes involved in the PPI networks are presented in Table IV, including proteasome (cluster 1; P=8.10×10−4), toll-like receptor signaling pathway (cluster 2; P=2.26×10−2), cytokine-cytokine receptor interaction (cluster 3; P=6.80×10−4), and T cell receptor signaling pathway (cluster 4; P=1.41×10−5; which involved IL2 and TNF).

Table III.

Top 10 functions enriched for the genes involved in protein-protein interaction networks.

Cluster	Description	Gene number	Gene symbol	P-value
1	GO:0051248~negative regulation of protein metabolic process	18	HSP90AB1, PSMB10, CLN3, PPP2R1A, NDUFA13, CDK5, PRKCD, FLNA, PSMB8, TGFB1, TIMP1, PSMB9, PSMC3, PSMB3, BAX, PSMD3, VPS28, PSMD8	2.39×10−9
	GO:0032268~regulation of cellular protein metabolic process	27	HSP90AB1, PSMB10, EIF5A, ITGB2, STUB1, TGFB1, TIMP1, NR1H2, PSMB3, PSMD3, PSMD8, CD28, CLN3, PPP2R1A, CD3E, NDUFA13, CDK5, PRKCD, RPS5, PSMB8, PSMB9, PSMC3, BAX, CD81, HSPB1, VPS28, PPP1R15A	7.78×10−9
	GO:0032269~negative regulation of cellular protein metabolic process	17	HSP90AB1, PSMB10, CLN3, PPP2R1A, NDUFA13, CDK5, PRKCD, PSMB8, TGFB1, TIMP1, PSMB9, PSMC3, PSMB3, BAX, PSMD3, VPS28, PSMD8	9.96×10−9
	GO:0031400~negative regulation of protein modification process	13	PSMB10, PPP2R1A, CDK5, PRKCD, PSMB8, TGFB1, PSMB9, PSMC3, BAX, PSMB3, PSMD3, VPS28, PSMD8	1.86×10−7
	GO:0006955~immune response	30	PSMB10, IFITM2, IFITM3, ACP5, CD70, IL32, OAS1, TGFB1, IFI35, MIF, TUBB, TMEM173, ZAP70, IL2RG, CD27, CD28, CD7, BST2, NCF4, HLA-B, PRKCD, PSMB8, HLA-G, BCAP31, PSMB9, GPI, CYBA, CORO1A, OASL, LIME1	3.38×10−7
	GO:0070271~protein complex biogenesis	23	ARL2, PPP2R1A, OXA1L, POLR2E, CD3E, AP2S1, ALDOC, POLR2J, POLR2I, ARPC4, CDK5, TGFB1, FLNA, MIF, CYBA, TUBB, BAX, ALOX5AP, GTF2F1, CAPG, VAMP3, TUBA1C, SCO2	5.68×10−6
	GO:0006461~protein complex assembly	23	ARL2, PPP2R1A, OXA1L, POLR2E, CD3E, AP2S1, ALDOC, POLR2J, POLR2I, ARPC4, CDK5, TGFB1, FLNA, MIF, CYBA, TUBB, BAX, ALOX5AP, GTF2F1, CAPG, VAMP3, TUBA1C, SCO2	5.68×10−6
	GO:0031399~regulation of protein modification process	17	PSMB10, PPP2R1A, CD3E, ITGB2, CDK5, PRKCD, STUB1, PSMB8, TGFB1, PSMB9, PSMC3, PSMB3, BAX, CD81, PSMD3, VPS28, PSMD8	8.00×10−6
	GO:0031397~negative regulation of protein ubiquitination	9	PSMB10, PSMC3, PSMB3, PSMD3, VPS28, CDK5, PSMD8, PSMB8, PSMB9	1.30×10−5
	GO:0065003~macromolecular complex assembly	26	OXA1L, POLR2E, ALDOC, POLR2J, AP2S1, POLR2I, ARPC4, TGFB1, MIF, TUBB, ALOX5AP, SCO2, TUBA1C, ARL2, PPP2R1A, CD3E, CDK5, FLNA, CYBA, PIH1D1, DGAT1, GTF2F1, BAX, CAPG, SNRPB, VAMP3	1.65×10−5
2	GO:0006952~defense response	7	MAVS, TNFAIP8L2, HDAC4, FOS, CCR5, IRF7, TLR1	2.84×10−3
	GO:0006954~inflammatory response	5	HDAC4, FOS, CCR5, IRF7, TLR1	6.92×10−3
	GO:0006853~carnitine shuttle	2	SLC25A20, CPT2	7.08×10−3
	GO:0045892~negative regulation of transcription, DNA-dependent	5	HDAC4, IRF7, FOSB, KLF4, DNAJB6	9.48×10−3
	GO:0051253~negative regulation of RNA metabolic process	5	HDAC4, IRF7, FOSB, KLF4, DNAJB6	1.00×10−2
	GO:0032365~intracellular lipid transport	2	SLC25A20, CPT2	1.64×10−2
	GO:0015838~betaine transport	2	SLC25A20, CPT2	1.64×10−2
	GO:0015879~carnitine transport	2	SLC25A20, CPT2	1.64×10−2
	GO:0006955~immune response	6	MAVS, TNFAIP8L2, GPR183, CCR5, TLR1, GPR65	2.18×10−2
	GO:0016481~negative regulation of transcription	5	HDAC4, IRF7, FOSB, KLF4, DNAJB6	2.22×10−2
3	GO:0010033~response to organic substance	15	BCL10, EIF2C2, STAT5A, TAF9B, HSPA1A, HSPA1B, CTNNB1, CYP7B1, ID2, TFRC, GNG10, PRKRA, HSPA6, JAK2, HSPD1, DNAJB4	2.20×10−5
	GO:0050867~positive regulation of cell activation	6	BCL10, IL5, STAT5A, BCL6, JAK2, HSPD1	3.07×10−4
	GO:0001817~regulation of cytokine production	7	BCL10, REL, STAT5A, BCL6, JAK2, HSPD1, IL1A	3.97×10−4
	GO:0006325~chromatin organization	9	HIST1H2AC, KDM2B, HIST2H2BE, HIST1H2BG, HIST1H2AE, MORF4L2, EED, HIST1H3D, BCOR, HIST1H3H	8.67×10−4
	GO:0006334~nucleosome assembly	5	HIST1H2AC, HIST2H2BE, HIST1H2BG, HIST1H2AE, HIST1H3D, HIST1H3H	1.03×10−3
	GO:0031497~chromatin assembly	5	HIST1H2AC, HIST2H2BE, HIST1H2BG, HIST1H2AE, HIST1H3D, HIST1H3H	1.18×10−3
	GO:0002761~regulation of myeloid leukocyte differentiation	4	IL5, ID2, STAT5A, CTNNB1	1.33×10−3
	GO:0010629~negative regulation of gene expression	10	EIF2C2, ID2, JARID2, TAF9B, PRKRA, EED, BCL6, BCOR, RBPJ, CTNNB1	1.36×10−3
	GO:0065004~protein-DNA complex assembly	5	HIST1H2AC, HIST2H2BE, HIST1H2BG, HIST1H2AE, HIST1H3D, HIST1H3H	1.39×10−3
	GO:0034728~nucleosome organization	5	HIST1H2AC, HIST2H2BE, HIST1H2BG, HIST1H2AE, HIST1H3D, HIST1H3H	1.51×10−3
4	GO:0042127~regulation of cell proliferation	20	IL4, NAMPT, TNF, FOSL2, IL8, PTGS2, KLF10, CTLA4, NFKBIA, IL13, SIRT1, SLAMF1, IL10, IL12RB2, LIF, CD80, MYC, PLAU, LTA, IL2	6.12×10−9
	GO:0016265~death	19	TRAF1, IER3, FOSL2, TNF, BCL2A1, NR4A2, NFKBIA, FASLG, BIRC3, SIRT1, AHR, GADD45G, ZC3H12A, SIAH2, TNFAIP3, GADD45B, GADD45A, MYC, LTA	1.07×10−8
	GO:0001775~cell activation	13	IL4, ZBTB32, ICAM1, TNF, IL8, RELB, IL21R, SLAMF1, IL10, CD80, CD40LG, LTA, IL2	1.67×10−8
	GO:0006915~apoptosis	17	TRAF1, IER3, TNF, BCL2A1, NFKBIA, FASLG, BIRC3, SIRT1, AHR, GADD45G, ZC3H12A, SIAH2, TNFAIP3, GADD45B, GADD45A, MYC, LTA	3.30×10−8
	GO:0006955~immune response	18	IL4, ICAM1, CCL3, TNF, IL18RAP, IL8, RELB, CTLA4, FASLG, IL13, GEM, CCL4, IL10, LIF, CD83, CD40LG, LTA, IL2	3.43×10−8
	GO:0012501~programmed cell death	17	TRAF1, IER3, TNF, BCL2A1, NFKBIA, FASLG, BIRC3, SIRT1, AHR, GADD45G, ZC3H12A, SIAH2, TNFAIP3, GADD45B, GADD45A, MYC, LTA	4.06×10−8
	GO:0042981~regulation of apoptosis	19	TRAF1, IL4, IER3, TNF, PTGS2, KLF10, BCL2A1, NR4A2, NFKBIA, FASLG, PIM3, BIRC3, SIRT1, IL10, CD40LG, TNFAIP3, MYC, LTA, IL2	5.39×10−8
	GO:0008219~cell death	18	TRAF1, IER3, FOSL2, TNF, BCL2A1, NFKBIA, FASLG, BIRC3, SIRT1, AHR, GADD45G, ZC3H12A, SIAH2, TNFAIP3, GADD45B, GADD45A, MYC, LTA	6.26×10−8
	GO:0043067~regulation of programmed cell death	19	TRAF1, IL4, IER3, TNF, PTGS2, KLF10, BCL2A1, NR4A2, NFKBIA, FASLG, PIM3, BIRC3, SIRT1, IL10, CD40LG, TNFAIP3, MYC, LTA, IL2	6.26×10−8
	GO:0010941~regulation of cell death	19	TRAF1, IL4, IER3, TNF, PTGS2, KLF10, BCL2A1, NR4A2, NFKBIA, FASLG, PIM3, BIRC3, SIRT1, IL10, CD40LG, TNFAIP3, MYC, LTA, IL2	6.62×10−8

Table IV.

The pathways enriched for the genes involved in protein-protein interaction networks (only the top 10 pathways are presented for cluster 4).

Cluster	Description	Gene number	Gene symbol	P-value
1	hsa03050:Proteasome	7	PSMB10, PSMC3, PSMB3, PSMD3, PSMD8, PSMB8, PSMB9	8.10×10−4
	hsa05016:Huntington's disease	13	UQCRC1, POLR2E, NDUFB7, CREB3, SLC25A6, AP2S1, POLR2J, COX8A, POLR2I, GPX1, BAX, NDUFS3, AP2M1	1.16×10−3
	hsa04142:Lysosome	9	CLN3, AP1M1, PSAP, ACP5, NEU1, CTSA, ATP6V0D1, ATP6V0B, GBA	6.71×10−3
	hsa04514:Cell adhesion molecules (CAMs)	8	ITGB7, ICAM2, ICAM3, ITGB2, HLA-B, SELPLG, HLA-G, CD28	3.85×10−2
	hsa04650:Natural killer cell mediated cytotoxicity	8	RAC2, ICAM2, LCK, CD247, ZAP70, ITGB2, HLA-B, HLA-G	3.98×10−2
	hsa05130:Pathogenic Escherichia coli infection	5	ARPC1A, ARPC1B, TUBB, ARPC4, TUBA1C	4.70×10−2
2	hsa04620:Toll-like receptor signaling pathway	3	FOS, IRF7, TLR1	2.26×10−2
3	hsa04060:Cytokine-cytokine receptor interaction	9	CCR8, TSLP, IL1R2, IL5, CCL20, CXCL3, CXCL2, IL1RAP, IL1A	6.80×10−4
	hsa04062:Chemokine signaling pathway	7	CCR8, BRAF, CCL20, CXCL3, GNG10, CXCL2, JAK2	2.71×10−3
	hsa05322:Systemic lupus erythematosus	5	HIST1H2AC, HIST2H2BE, HIST1H2BG, HIST1H2AE, HIST1H3D, HIST1H3H	6.60×10−3
	hsa04640:Hematopoietic cell lineage	4	IL1R2, IL5, TFRC, IL1A	2.76×10−2
	hsa03010:Ribosome	4	RPL23, RPL5, RPL37A, RPS24	2.84×10−2
4	hsa04060:Cytokine-cytokine receptor interaction	15	IL4, CCL3, TNF, IL18RAP, IL8, IL21R, FASLG, IL13, CCL4, IL10, IL12RB2, LIF, CD40LG, LTA, IL2	3.12×10−9
	hsa05330:Allograft rejection	7	IL4, TNF, CD80, CD40LG, FASLG, IL10, IL2	2.24×10−7
	hsa05320:Autoimmune thyroid disease	7	IL4, CD80, CD40LG, CTLA4, FASLG, IL10, IL2	1.91×10−6
	hsa04620:Toll-like receptor signaling pathway	8	CCL3, TNF, CD80, IL8, MAP2K3, MAP3K8, NFKBIA, CCL4	9.04×10−6
	hsa04660:T cell receptor signaling pathway	8	IL4, TNF, CD40LG, MAP3K8, CTLA4, NFKBIA, IL10, IL2	1.41×10−5
	hsa04630:Jak-STAT signaling pathway	9	LIF, IL12RB2, IL4, SPRY1, IL21R, IL13, MYC, IL10, IL2	1.71×10−5
	hsa04010:MAPK signaling pathway	11	DUSP5, TNF, DUSP2, MAP2K3, RELB, MAP3K8, GADD45G, FASLG, GADD45B, GADD45A, MYC	2.21×10−5
	hsa05310:Asthma	5	IL4, TNF, CD40LG, IL13, IL10	6.12×10−5
	hsa04940:Type I diabetes mellitus	5	TNF, CD80, FASLG, LTA, IL2	2.68×10−4
	hsa04672:Intestinal immune network for IgA production	5	IL4, CD80, CD40LG, IL10, IL2	4.89×10−4

Discussion

In the present study, cluster 1 (including 292 genes), cluster 2 (including 111 genes), cluster 3 (including 194 genes) and cluster 4 (including 103 genes) were obtained from soft clustering analysis. Genes in cluster 1 were significantly downregulated after 200 min and were significantly upregulated after 300 min. Genes in cluster 2 and cluster 3 had evidently opposite tendencies. Genes in cluster 4 were significantly upregulated prior to 80 min and expression plateaued thereafter. The miRNA-mRNA pairs, TF-mRNA pairs and PPI relationships were searched, respectively. There were no miRNAs significantly enriched for the genes in cluster 2, and no transcriptional regulatory relationships were determined for the genes in cluster 3. There were 656, 40, 129 and 234 PPIs for genes in clusters 1, 2, 3 and 4, respectively. In particular, GAPDH (degree, 40) in cluster 1, JAK2 (degree, 10) and STAT5A (degree, 9) in cluster 3, and TNF (degree, 26) and IL2 (degree, 22) in cluster 4 exhibited high degrees in the PPI networks.

The gapA gene that encodes GAPDH is conserved in numerous serotypes of Haemophilus parasuis, and the GAPDH (pCgap) DNA vaccine may contribute to the immune response and inhibit infection with H. parasuis (34). A previous study demonstrated that GAPDH in Streptococcus agalactiae can function as a virulence-associated immunomodulatory protein (35). As a component of heterochromatin complexes, TRIM28 is phosphorylated following stimulation by the T cell antigen receptor, and is implicated in T cell activation and tolerance (36). In the transcriptional regulatory network for genes in cluster 1, GAPDH was targeted by TRIM28, indicating that TRIM28 may have a role in Treg function through targeting GAPDH. JAK2 propagates receptor-binding signals through inflammatory cytokines, and can serve as a relevant biological target in the control of allograft rejection or grading acute graft-versus-host disease without broader immune impairment (37). Genome-wide association studies have reported that the JAK-STAT signaling pathway is highly correlated with human autoimmunity, and targeting various JAKs has been applied in immune-mediated disease (38). These findings indicated that JAK2 may also be associated with Treg function.

In the miRNA-mRNA regulatory network for genes in cluster 4, TNF was targeted by miR-146b-3p, miR-146b-5p and miR-142-5p. Suppression or stimulation of the costimulators of TNF receptor family (TNFR) members may be used to treat cancer, autoimmunity, infectious disease and transplantation (39). Through quantitative polymerase chain reaction and flow cytometry, previous studies have indicated that anti-TNF antibody (infliximab) can increase FOXP3 expression in CD4+CD25high Tregs and restore the suppressive function of Tregs, thus suggesting that TNF may have a role in controlling autoimmunity through suppressing CD4+CD25+ Treg activity (40,41). In addition, the stimulation of glucocorticoid-induced TNFR-related protein conquers self-tolerance/ignorance and promotes T cell-mediated antitumor activity with minimal autoimmunity (42,43). miR-146a is ubiquitously expressed in Tregs and has an important role in congenital and acquired immune responses (44,45). The results of a western blot analysis and enzyme-linked immunosorbent assay indicated that miR-142-3p controls the levels of cyclic adenosine monophosphate via regulating adenylyl cyclase 9 mRNA in CD4+CD25+ Treg cells and CD4+CD25− T cells (46,47). Therefore, miR-146b-3p, miR-146b-5p and miR-142-5p may serve roles in Treg function via regulating TNF.

STAT5 binds to the FOXP3 promoter, indicating that activation of IL-2 receptor β-dependent STAT5 contributes to Treg differentiation by mediating FOXP3 expression (48). STAT5a/b have been demonstrated to serve a nonredundant, essential role in regulating Tregs, and have an opposing role compared with STAT3 in regulating FOXP3 (49). Flow cytometric analysis indicated that STAT5B transfers a crucial IL-2-mediated signal that promotes the accumulation of functional Tregs in vivo (50). STAT5 activation maintains the expression of FOXP3 in CD4+CD25− effector T cells and Tregs, thus indicating the influential role of cytokines on FOXP3 expression (51). These findings indicated that IL2 and STAT5A may be implicated in Treg function. Enrichment analyses demonstrated that IL2 and TNF were enriched in immune response and T cell receptor signaling pathway, suggesting that IL2 and TNF may affect Treg function via the immune response and T cell receptor signaling pathway.

In conclusion, cluster 1 (including 292 genes), cluster 2 (including 111 genes), cluster 3 (including 194 genes) and cluster 4 (including 103 genes) were obtained from a soft clustering analysis. Subsequently, GAPDH was revealed to be targeted by TRIM28, and TNF was targeted by miR-146b-3p, miR-146b-5p and miR-142-5p; these interactions in addition to JAK2, IL2 and STAT5A may have important roles in Treg function. However, these findings, which were obtained by bioinformatics analysis, require further experimental verification.