Identification of potential key genes associated with severe pneumonia using mRNA-seq.

This study aimed to identify the potential key genes associated with severe pneumonia using mRNA-seq. Nine peripheral blood samples from patients with severe pneumonia alone (SP group, n=3) and severe pneumonia accompanied with chronic obstructive pulmonary disease (COPD; CSP group, n=3), as well as volunteers without pneumonia (control group, n=3) underwent mRNA-seq. Based on the sequencing data, differentially expressed genes (DEGs) were identified by Limma package. Following the pathway enrichment analysis of DEGs, the genes that were differentially expressed in the SP and CSP groups were selected for pathway enrichment analysis and coexpression analysis. In addition, potential genes related to pneumonia were identified based on the information in the Comparative Toxicogenomics Database. In total, 645 and 528 DEGs were identified in the SP and CSP groups, respectively, compared with the normal controls. Among these DEGs, 88 upregulated genes and 80 downregulated genes were common between the two groups. The functions of the common DEGs were similar to those of the DEGs in the SP group. In the coexpression network, the commonly downregulated genes (including ND1, ND3, ND4L, and ND6) and the commonly upregulated genes (including TSPY6P and CDY10P) exhibited a higher degree. In addition, 131 DEGs (including ND1, ND3, ND6, MIR449A and TAS2R43) were predicted to be potential pneumonia-related genes. In conclusion, the present study demonstrated that the common DEGs may be associated with the progression of severe pneumonia.

Introduction

Pneumonia is defined as the inflammation and consolidation of lung tissue due to an infectious agent. It is a leading cause of mortality worldwide, particularly in children and the elderly (1,2). Typical symptoms of pneumonia include fever, chills, pleuritic chest pain, and productive cough of purulent sputum (1). Severe pneumonia is defined as the admission to an intensive care unit (ICU) and results in extremely high mortality rates (3). Chronic obstructive pulmonary disease (COPD) is a major public health problem and is a leading cause of chronic morbidity and mortality worldwide (4). COPD is particularly prevalent in the elderly and it is a frequent comorbidity in patients with community-acquired pneumonia (CAP) (5). During treatment with corticosteroids in elderly patients with COPD, serious pneumonia typically occurs with a high risk (6,7). The severity of COPD, including the presence of pulmonary emphysema and the airflow limitation grade, has been demonstrated to be independently associated with the development of severe pneumonia (8). Furthermore, CAP patients with COPD have a similar mortality rate compared to those without COPD (9). Therefore, it is necessary to elucidate the molecular mechanism underlying severe pneumonia accompanied with COPD and compare it with the mechanism underlying severe pneumonia.

Several molecular mechanisms associated with severe pneumonia have been described. For instance, four risk single-nucleotide polymorphisms located in chromosomes 1 and 17 have been identified as being significantly correlated with the susceptibility to developing severe pneumonia in A/H1N1 infection (10). Staphylococcal α-hemolysin is an essential virulence factor in severe Staphylococcus aureus pneumonia and α-hemolysin-mediated activation of the NLRP3 inflammasome may induce necrotic pulmonary injury (11). Furthermore, previous reports have indicated that severe pneumonia is associated with methicillin-resistant S. aureus carrying the staphylococcal cassette chromosome mec type IV and Panton-Valentine leukocidin genes (12,13). Furthermore, a recent study reported that high expression of interleukin-10 (IL-10) and interferon-induced protein (IP-10) in human immunodeficiency virus-infected infants were associated with more severe hypoxic pneumonia (14). The IL-6–174 GG genotype is correlated with lower severity and mortality in patients with pneumococcal CAP (15). However, to date, there have been no studies reporting the molecular mechanisms underlying severe pneumonia accompanied with COPD, or comparing the similarities and differences between patients with severe pneumonia alone and severe pneumonia patients with COPD at the genetic level.

In the present study, mRNA-seq and bioinformatics were used to analyze the differential expression of genes in the peripheral blood from patients with severe pneumonia alone and severe pneumonia accompanied with COPD. Genes that were differentially expressed in the two groups of patients were analyzed. These findings may provide novel information for the understanding of the molecular mechanisms underlying severe pneumonia and severe pneumonia accompanied with COPD.

Materials and methods

Clinical samples

This study was approved by the Medical Ethics Committee of the Chinese People's Liberation Army (PLA) General Hospital (Beijing, China). In total, 18 patients with pneumonia who received therapy in Chinese PLA General Hospital between June 2013 and December 2013 were included in this study, including 9 patients with severe pneumonia alone (the SP group), 9 patients with severe pneumonia accompanied with COPD (the CSP group). Another 9 volunteers without pneumonia were enrolled as normal controls (the C group). For inclusion, patients with severe pneumonia met at least one of the following criteria: i) Disorders of consciousness; ii) respiratory rate ≥30/min; iii) diastolic blood pressure <60 mmHg, PaO2/FiO2 <300, and mechanical ventilation; iv) systolic blood pressure ≤90 mmHg; v) septic shock; vi) bilateral or multilobar pneumonia by chest radiograph, or lesions enlargement within 48 h after admission ≥50%; and vii) oliguria: urine volume <20 ml/h or <80 ml/4 h, or acute renal failure requiring dialysis treatment (16).

Peripheral blood samples were collected from each patient and volunteer. Written informed consent was obtained before sampling. In each group, three samples with an equal volume were randomly pooled into one sample for sequencing. Therefore, three sequencing samples were generated for each group: WLL1–3 for the SP group, WLL4–6 for the CSP group, and WLL7-9 for the C group.

RNA extraction and Illumina sequencing

Firstly, plasma was separated from the nine sequencing samples, respectively. Total RNA was extracted and purified from the plasma samples using an miRNeasy Serum/Plasma Kit (Qiagen GmbH, Hilden, Germany). The quality control of RNA was measured using the Agilent Bioanalyzer 2100 (Agilent Technologies, Inc., Palo Alto, CA, USA). The bioanalyzer software automatically generated the value of RNA Integrity Number (RIN, 1 to 10) based on the ratio of the 18S to 28S ribosomal subunits to determine the level of RNA degradation in gel electrophoresis, which removed individual interpretation in RNA quality control. RNAs with RIN ≥8.0 were used in the study. Subsequently, poly (A) mRNA was enriched by oligo (dT) magnetic beads (Dynabeads® oligo dT, Dynal Inc., Great Neck, NY), and then interrupted into shot fragments by fragmentation buffer (Agilent Technologies, Inc., Santa Clara, CA, USA). RNA fragments were then reverse transcribed into cDNA. The concentration of cDNA in the library was quantified into 1 ng/µl with a Qubit 2.0 fluorometer, and then the cDNA were detected using the Agilent Bioanalyzer 2100 (Agilent Technologies, Inc.). Libraries were pooled according to the data size and effective cDNA concentration. Clusters of the cDNA libraries were generated on an Illumina cBot system (Illumina, Inc., San Diego, CA, USA). Finally, the cDNA libraries were sequenced using the rapid mode 2×150 bp paired-end reads on an Illumina HiSeq™ 4000 system (Illumina, Inc.). Raw sequencing data was uploaded to the public National Center for Biotechnology Information database under the BioProject accession no. PRJNA325383.

Data filtering

Adapter sequences in the raw reads were removed using Cutadapt (version 1.2.1; cutadapt.readthedocs.io/en/stable/) (17), with at least 10-bp overlap (AGATCGGAAG) and allowing for a 20% base error rate. Using a 5-bp window, the mean quality value of each window was controlled to be at least 20. Shortening of the 3′ terminal reads was permitted, although the shortest length was at least 50 bp. Reads with unknown sequences (‘N’) were removed.

Statistics and alignment of reads

The number of raw and clean reads, as well as the percentage of clean reads, was summarized to ensure the validity and reliability of the sequencing data. Furthermore, the index of the reference genome was created using Bowtie 2 (version 2.1.0; sourceforge.net/projects/bowtie-bio/files/bowtie2/2.1.0/), and clean reads were aligned to the human genome (hg19) using TopHat 2.1.1 software (ccb.jhu.edu/software/tophat/index.shtml).

Differential expression analysis of genes

The read count of each gene was calculated by HTSeq 0.6.1p2 (huber.embl.de/users/anders/HTSeq) and subsequently normalized using Trimmed Mean of M values normalization (18) in the Limma package (version 3.24.15; bioconductor.org/packages/release/bioc/html/limma.html). Normalized data were transformed to gene expression matrix using the voom method (19).

Subsequently, differentially expressed genes (DEGs) in the comparison groups of SP vs. C and CSP vs. C were identified using the edgeR package of Bioconductor (bioconductor.org/packages/release/bioc/html/edgeR.html). Only the genes meeting the criteria of |log2 FC (fold change) |>1 and P<0.05 were identified as DEGs.

Functional analysis of the DEGs

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed for the identified DEGs using the clusterProfiler (version 2.2.7; bioconductor.org/packages/release/bioc/html/clusterProfiler.html) package of R (20). The P-value of each pathway term was calculated by Fisher's exact test (21), and only the pathway terms with P<0.05 were considered significant.

Venn analysis of the DEGs

Gene symbols and their log2 FC values of the identified DEGs were input into the VennPlex software (version 1.0.02; irp.nia.nih.gov/bioinformatics/vennplex.html) (22) to compare the DEGs in the SP and CSP groups. Common upregulated and downregulated genes between the two groups were identified and defined as common DEGs.

Coexpression analysis of the common DEGs

Pearson correlation coefficient (PCC) was used to screen the coexpression pairs among all of DEGs. P-value was calculated based on the Z-score of the coexpression pair, and |PCC| >0.9 and P<0.05 were set as the cut-off criteria.

Furthermore, coexpression pairs of the common DEGs were selected for the construction of the coexpression network. Here, each pair must contain at least one common DEG. The network was visualized by Cytoscape 3.3.0 (cytoscape.org/). In the network, a ‘node’ represents a protein (or gene); a ‘line’ represents a coexpression pair of the proteins. The ‘degree’ of each node is equal to the number of nodes that coexpress with this node. The higher the degree is, the closer the connections with other nodes are, indicating a more important role of the node in the network.

Analysis of the potential pneumonia-related genes

Comparative Toxicogenomics Database (CTD; ctdbase.org/) provides manually curated information about chemical-gene/protein interactions, chemical-disease and gene-disease relationships (23). In this study, the keyword ‘pneumonia’ was used to search the genes associated with pneumonia (April 26, 2016). The DEGs in the SP and CSP groups were then compared with the pneumonia-related genes in CTD using a venn tool (bioinformatics.lu/venn.php) in order to identify the potential pneumonia-related genes from the DEGs.

Results

Data summary of quality control and sequence alignment

Following quality control of the raw sequencing data, the percentage of clean data was >97% in the nine samples (Table I), indicating the sequencing data as of high quality. The sequence alignment of the clean reads to the hg19 indicated that overall mapping rates in the majority of samples were between 60 and 75%. The concordant pair alignment rate was 50–65% in the majority of samples (Table II).

Table I.

Summary of the sequencing data after quality control.

Sample	Raw data (Mb)	Clean data (Mb)	Percent (%)
WLL1	15857.06	15457.13	97.48
WLL2	30910.76	30140.40	97.51
WLL3	14380.84	13989.76	97.28
WLL7	16242.65	15811.34	97.34
WLL8	17571.73	17114.66	97.40
WLL9	12614.38	12283.53	97.38
WLL4	12177.13	11909.79	97.80
WLL5	21674.30	21182.59	97.73
WLL6	19091.95	18678.84	97.84

WLL1-3, samples from the patients with severe pneumonia alone

WLL4-6, samples from the severe pneumonia patients accompanied with chronic obstructive pulmonary disease; WLL7-9, control samples. Percent=clean data/raw data.

Table II.

Data summary of the sequence alignment.

Sample	Left mapped reads	Right mapped reads	Overall mapping rate (%)	Concordant alignment pair rate (%)
WLL1	40655723	31737852	68.50	53.20
WLL2	68307024	48549459	56.70	39.30
WLL3	36746057	28219712	67.80	52.80
WLL4	31492239	27325889	72.50	62.10
WLL5	57866494	49099958	74.00	62.30
WLL6	50706071	42750383	73.40	61.40
WLL7	38850721	28955703	62.60	46.50
WLL8	45900593	36937197	70.70	57.80
WLL9	32531013	25214426	68.70	55.00

WLL1-3, samples from the patients with severe pneumonia alone

WLL4-6, samples from the severe pneumonia patients accompanied with chronic obstructive pulmonary disease; WLL7-9, control samples. Left mapped reads and right mapped reads represent the number of reads in both terminals that are mapped to the reference genome.

Identification of DEGs

Based on the criteria, a total of 645 DEGs (323 upregulated and 322 downregulated genes) were identified in the SP group, along with 528 DEGs (292 upregulated and 236 downregulated genes) in the CSP group, compared with the C group. These DEGs were able to distinguish the two group samples (Fig. 1A and B).

Figure 1.

Heat maps of the DEGs. (A) Heat map of DEGs between the severe pneumonia alone group and the control group. (B) Heat map of DEGs between the severe pneumonia with COPD group and the control group. Green indicates downregulated and red indicates upregulated. DEGs, differentially expressed genes; COPD, chronic obstructive pulmonary disease. WLL1-3, samples taken from patients with severe pneumonia alone; WLL4-6, samples taken from severe pneumonia patients accompanied with COPD; WLL7-9, samples taken from the normal controls.

Pathway enrichment analysis of DEGs

To further investigate the potential biological functions of the identified DEGs, KEGG pathway enrichment analysis was performed. For the DEGs in the SP group, the upregulated gene HLA-DRB5 was significantly enriched in the pathways associated with asthma, allograft rejection and intestinal immune network for IgA production, which were related to immunity. The downregulated genes were mainly enriched in the pathways associated with aminoacyl-tRNA biosynthesis (TRNS1 and TRNM) and oxidative phosphorylation (COX3 and ND4L; Table III).

Table III.

Enriched pathway terms of differentially expressed genes in the patients with severe pneumonia alone compared with the controls.

Pathway ID	Description	P-value	Genes
Upregulated
hsa05310	Asthma	3.12E-02	HLA-DRB5
hsa05330	Allograft rejection	3.91E-02	HLA-DRB5
hsa05332	Graft-versus-host disease	4.30E-02	HLA-DRB5
hsa04940	Type I diabetes mellitus	4.50E-02	HLA-DRB5
hsa04672	Intestinal immune network for IgA production	4.89E-02	HLA-DRB5
Downregulated
hsa00970	Aminoacyl-tRNA biosynthesis	1.18E-30	TRNS1, TRNM, TRNW, TRNG, TRNH, TRNL1, TRNN,
			TRNC, TRNA, TRNP, TRNE, TRNT, TRNV, TRNS2, TRND,
			TRNI, TRNQ, TRNK, TRNY, TRNL2, TRNF, TRNR
hsa05012	Parkinson's disease	1.35E-10	COX3, ND4L, COX2, ND1, COX1, ND4, ATP8,
			CYTB, ND3, ND5, ND2, ATP6, ND6
hsa00190	Oxidative phosphorylation	1.64E-10	COX3, ND4L, COX2, ND1, COX1, ND4, ATP8,
			CYTB, ND3, ND5, ND2, ATP6, ND6
hsa05016	Huntington's disease	1.59E-03	GPX1, COX3, COX2, COX1, ATP8, CYTB, ATP6
hsa05010	Alzheimer's disease	4.89E-03	COX3, COX2, COX1, ATP8, CYTB, ATP6
hsa04260	Cardiac muscle contraction	5.92E-03	COX3, COX2, COX1, CYTB

Furthermore, for the DEGs in the CSP group, the upregulated genes (such as UGT2B17) were markedly enriched in the pathways of metabolism, including starch and sucrose metabolism, and ascorbate and aldarate metabolism. The downregulated genes were significantly enriched in the aminoacyl-tRNA biosynthesis (TRNS1 and TRNM) and oxidative phosphorylation (COX3 and ND4L) pathways (Table IV), which is similar to those of the downregulated genes in the SP group.

Table IV.

Enriched pathway terms of differentially expressed genes in the severe pneumonia patients accompanied with chronic obstructive pulmonary disease compared with the controls.

Pathway ID	Description	P-value	Genes
Upregulated
hsa00500	Starch and sucrose metabolism	2.23E-03	UGT2B17, AMY1C
hsa00980	Metabolism of xenobiotics by cytochrome P450	3.82E-03	UGT2B17, GSTM1
hsa00982	Drug metabolism-cytochrome P450	4.04E-03	UGT2B17, GSTM1
hsa00053	Ascorbate and aldarate metabolism	3.48E-02	UGT2B17
hsa05310	Asthma	4.13E-02	HLA-DRB5
hsa00040	Pentose and glucuronate interconversions	4.26E-02	UGT2B17
Downregulated
hsa00970	Aminoacyl-tRNA biosynthesis	5.87E-36	TRNL1, TRNS2, TRNR, TRNV, TRNH, TRNC,
			TRNN, TRNP, TRNY, TRNT, TRNG, TRNE, TRNQ,
			TRNI, TRNF, TRNS1, TRNW, TRNA, TRND,
			TRNK, TRNL2, TRNM
hsa05012	Parkinson's disease	3.63E-13	COX1, COX2, COX3, ND1, ATP6, ND4L, ND4,
			ATP8, ND3, ND2, CYTB, ND5, ND6
hsa00190	Oxidative phosphorylation	4.43E-13	COX1, COX2, COX3, ND1, ATP6, ND4L, ND4,
			ATP8, ND3, ND2, CYTB, ND5, ND6
hsa05010	Alzheimer's disease	5.47E-04	COX1, COX2, COX3, ATP6, ATP8, CYTB
hsa05016	Huntington's disease	8.60E-04	COX1, COX2, COX3, ATP6, ATP8, CYTB
hsa04260	Cardiac muscle contraction	1.28E-03	COX1, COX2, COX3, CYTB
hsa01100	Metabolic pathways	1.63E-02	COX1, COX2, COX3, ND1, ATP6, ND4L, ND4,
			ATP8, ND3, ND2, CYTB, ND5, ND6

Venn analysis of the DEGs and the functions of the common DEGs

To further investigate whether there are genes that are differentially expressed in the two pneumonia groups, the DEGs in the SP and CSP groups were compared. According to the venn diagram, 88 DEGs were upregulated in the two groups (named common-upregulated genes) and 80 were downregulated in the two groups (named common-downregulated genes; Fig. 2A).

Figure 2.

Venn diagram (A) and results of pathway enrichment analysis (B) for the differentially expressed genes. SP, severe pneumonia; CSP, severe pneumonia with chronic obstructive pulmonary disease; C, control; Up, upregulated genes; Down, downregulated genes.

Furthermore, pathway enrichment analysis of the common DEGs revealed that the commonly upregulated gene, HLA-DRB5, was significantly enriched in the pathways associated with asthma, allograft rejection and intestinal immune network for IgA production. The commonly downregulated genes were enriched in the pathways aminoacyl-tRNA biosynthesis (TRNS1 and TRNM) and oxidative phosphorylation (COX3 and ND4L) pathways (Fig. 2B). The functions of the common DEGs were similar to those of the DEGs in the SP group.

Analysis of the coexpression network

In total, 805 coexpression pairs of common genes were chosen to construct a network. The coexpression network consisted of 245 nodes and 805 coexpression pairs. A series of 69 commonly upregulated genes and 65 commonly downregulated genes were included in the network (Fig. 3).

Figure 3.

Coexpression network of the differentially expressed genes. White nodes represent the upregulated genes and grey nodes represent the downregulated genes. Round nodes represent the genes that are differentially expressed only in the patients with severe pneumonia alone. Quadrate nodes represent the genes that are differentially expressed only in the severe pneumonia patients with chronic obstructive pulmonary disease. Triangular nodes represent the genes that are differentially expressed in both patients with severe pneumonia alone and severe pneumonia patients with chronic obstructive pulmonary disease. The solid line represents the positive correlation, namely, the Pearson correlation coefficient >0; the dotted line represents the negative correlation, namely, the Pearson correlation coefficient <0.

Some genes had a higher degree in the network, including ND3 (degree, 36), ND4L (degree, 36), TRNT (degree, 35), ND6 (degree, 35), TRNP (degree, 35), ATP8 (degree, 35), ND1 (degree, 35), MTND2P28 (degree, 35) and MTATP6P1 (degree, 34). All of these genes were the commonly downregulated genes. Among the commonly upregulated genes, several genes exhibited high degrees, including TSPY6P (degree, 15), CDY10P (degree, 12), TTTY15 (degree, 11), RNA5SP233 (degree, 11), OFD1P3Y (degree, 9) and RN7SL255P (degree, 8; Fig. 3).

Identification of potential pneumonia-related genes

Based on the information in CTD, 131 DEGs were identified to be associated with pneumonia, including 77 DEGs in the SP group and 54 DEGs in the CSP group (Table V). Among these DEGs, 4 commonly upregulated genes (MIR449A, TAS2R43, MIR656, and EIF5AP4) and 16 commonly downregulated genes (NOXO1, COX1, COX2, COX3, ND4, ND1, ND2, ND6, CYTB, ND3, SH3BGRL, ND5, ATP6, CPXCR1, VCX3B, and CT45A3) were included.

Table V.

Potential pneumonia-related genes identified from the differentially expressed genes based on the Comparative Toxicogenomics Database.

Category	Potential pneumonia-related genes
Severe pneumonia alone group	PF4, TIMP1, GPX1, B2M, FTH1, NOXO1, PPBP, SDPR, TMSB4X, COX1, PCBP1,
	SH3BGRL3, NRGN, COX2, LRRC26, RGS18, FLNA, ITGB3, TCEAL7, SAA4, ACOT1,
	HAPLN4, COX3, ND4, HBA1, HLA-E, HIST1H2AH, ND1, PGRMC1, TUBB, ND2,
	RPL15, MIR148B, ND6, MIR205, CYTB, APOBEC3B, ND3, PCDH20, MIR449A,
	SH3BGRL, HLA-DRB5, USP17L9P, ND5, MIR193A, HBA2, SNORA12, MIR559,
	MIR572, TAS2R43, CSNK2B, SPANXB1, SLX1A, FTHL17, LY6G6E, SLC10A3,
	MIR339, RNR1, MIR513B, HLA-DRB6, MIR4300, MIR656, ATP6, POTEF, CPXCR1,
	RN7SK, FAM138A, RPS18P9, SNORD115-30, MIR873, MIR190A, RBMY2EP, RNU1-1,
	EIF5AP4, VCX3B, MIR4499, CT45A3
Severe pneumonia accompanied with chronic obstructive pulmonary disease group	GSTM1, JUND, NOXO1, COX1, COX2, COX3, MRPL12, ND4, ND1, SPANXD, ND2,
	RPS20, ND6, CYTB, MDP1, ND3, EIF4EBP3, MIR93, HNRNPU, MIR449A, MIR100,
	SH3BGRL, RPS4Y2, HLA-DRB5, IGKV1-17, ND5, UGT2B17, TAS2R43, MT1P3,
	POM121L8P, MIR188, LCE3C, GLOD5, CYB5D1, HBG1, MED14OS, INSL5, MIR656,
	SNORD113-3, ATP6, DTX2P1-UPK3BP1-PMS2P11, CPXCR1, DUX4L4, SNORA28,
	MIR190A, F8A3, MIR19B2, OR7E125P, EIF5AP4, MIR216A, VCX3B, TRBV21-1,
	RNU4-2, CT45A3

Discussion

In the present study, based on the mRNA-seq and bioinformatics findings, a total of 645 and 528 DEGs were identified in the patients with severe pneumonia alone and severe pneumonia patients with COPD, respectively, compared with the normal controls. Among these DEGs, 88 upregulated genes and 80 downregulated genes were common between the two groups. In the coexpression network, the commonly downregulated genes (including ND1, ND3, ND4L, and ND6) and the commonly upregulated genes (including TSPY6P and CDY10P) had a higher degree. In addition, 131 DEGs (including COX1, COX2, COX3, ND1, ND3, and ND6) were predicted to be potential pneumonia-related genes.

ND1, ND3, ND4L and ND6 encode the subunits of NADH dehydrogenase, which participate in the mitochondrial oxidative phosphorylation, functioning in the transfer of electrons from NADH to the respiratory chain (24). A previous study has demonstrated that the type I NADH dehydrogenase of the human pathogen Mycobacterium tuberculosis (Mtb) is able to neutralize NOX2-derived reactive oxygen species to inhibit tumor necrosis factor-α-mediated host macrophage apoptosis in innate immunity (25). Differentially expressed NADH dehydrogenase subunits were detected during the immune response of Ostrea edulis against bonamiosis (26). These results indicate the association of NADH dehydrogenase with immunity, which is an essential factor in the progression of pneumonia. In the present study, ND1, ND3 and ND6 were identified as potential pneumonia-related genes. Therefore, these genes encoding NADH dehydrogenase subunits may have crucial roles in the progression of severe pneumonia.

Commonly downregulated genes, COX1, COX2 and COX3, were also predicted as potential pneumonia-related genes. These three genes all encode isoforms of cyclooxygenase, which has a critical role in various normal physiological functions and regulates a variety of pathological conditions (27). Previous studies have demonstrated the associations of COX2 with inflammation and immunity (27,28). In the immune response of pneumonia, COX2 expression is induced in lung lesions (27,29,30), which is inconsistent with the result of the current study. This disparity may be due to the difference of sample types used for COX2 detection. In our future studies, COX2 expression will be determined in lung tissue and peripheral blood of patients with severe pneumonia and normal control subjects.

With the exception of the downregulated genes, several commonly upregulated genes also had a higher degree in the coexpression network, including TSPY6P and CDY10P, which are pseudogenes. In recent years, the biological functions of pseudogenes have been increasingly investigated, and pseudogenes have been reported to serve a role in formation of tumors (31). However, to the best of our knowledge, there have been no reports on the functions of TSPY6P and CDY10P. Therefore, the present study was the first to demonstrate TSPY6P and CDY10P were deregulated in patients with severe pneumonia, and were thus worthy of further research.

Furthermore, in this study, 4 commonly upregulated genes (MIR449A, TAS2R43, MIR656 and EIF5AP4) were identified as potential pneumonia-related genes. MIR449A and MIR656 encode miR-449a and miR-656, respectively. A previous study has demonstrated that miRNA-449a modulates the expression of the inflammatory marker, YKL40, by targeting components of the NOTCH signaling pathway, indicating the important role of miR-449a in the inflammatory response (32). TAS2R43 encodes taste 2 receptor member 43, which belongs to the large TAS2R receptor family. An agonist for TAS2R1 and other TAS2Rs, amarogentin, is able to decrease IL-8 and MMP-1 expression in mast cells (33), suggesting the potential association of TAS2R43 with immunity. Therefore, MIR449A and TAS2R43 may have a critical role in the progression of severe pneumonia, possibly via participating in inflammation and immunity. However, MIR656 and EIF5AP4 have not yet been reported to be related to pneumonia, thus, they are novel genes that are potentially associated with severe pneumonia.

Although a set of upregulated and downregulated genes were common in the samples from patients with severe pneumonia alone and patients with severe pneumonia accompanied with COPD, there were a series of genes that were differentially expressed between the two patient groups according to the venn diagram analysis performed in the present study. The differences in molecular mechanisms between severe pneumonia and severe pneumonia accompanied with COPD will be further investigated in our future studies.

In conclusion, the downregulated genes encoding NADH dehydrogenase subunits (ND1, ND3, ND4L and ND6) and the upregulated genes (TSPY6P, CDY10P, MIR449A, TAS2R43, MIR656, and EIF5AP4) are newly-identified genes that may be associated with the progression of severe pneumonia. These genes will be further researched in our future studies. These results provide novel information for further experimental studies to build upon and contribute to the continued understanding of the molecular mechanisms underlying severe pneumonia.