RNA‑sequencing analysis of aberrantly expressed long non‑coding RNAs and mRNAs in a mouse model of ventilator‑induced lung injury.

Long non-coding RNAs (lncRNAs) are closely associated with the regulation of various biological processes and are involved in the pathogenesis of numerous diseases. However, to the best of our knowledge, the role of lncRNAs in ventilator‑induced lung injury (VILI) has yet to be evaluated. In the present study, high‑throughput sequencing was applied to investigate differentially expressed lncRNAs and mRNAs (fold change >2; false discovery rate <0.05). Bioinformatics analysis was employed to predict the functions of differentially expressed lncRNAs. A total of 104 lncRNAs (74 upregulated and 30 downregulated) and 809 mRNAs (521 upregulated and 288 downregulated) were differentially expressed in lung tissues from the VILI group. Gene ontology analysis demonstrated that the differentially expressed lncRNAs and mRNAs were mainly associated with biological functions, including apoptosis, angiogenesis, neutrophil chemotaxis and skeletal muscle cell differentiation. The top four enriched pathways were the tumor necrosis factor (TNF) signaling pathway, P53 signaling pathway, neuroactive ligand‑receptor interaction and the forkhead box O signaling pathway. Several lncRNAs were predicted to serve a vital role in VILI. Subsequently, three lncRNAs [mitogen‑activated protein kinase kinase 3, opposite strand (Map2k3os), dynamin 3, opposite strand and abhydrolase domain containing 11, opposite strand] and three mRNAs (growth arrest and DNA damage‑inducible α, claudin 4 and thromboxane A2 receptor) were measured by reverse transcription‑quantitative polymerase chain reaction, in order to confirm the veracity of RNA‑sequencing analysis. In addition, Map2k3os small interfering RNA transfection inhibited the expression of stretch‑induced cytokines [TNF‑α, interleukin (IL)‑1β and IL‑6] in MLE12 cells. In conclusion, the results of the present study provided a profile of differentially expressed lncRNAs in VILI. Several important lncRNAs may be involved in the pathological process of VILI, which may be useful to guide further investigation into the pathogenesis for this disease.

Introduction

Mechanical ventilation (MV) is an irreplaceable therapy for patients in intensive care units and emergency departments. Although the techniques for MV have significantly improved in the past 20 years (1), the risk of aggravating lung injury in patients with acute respiratory distress syndrome (ARDS) (2), and damaging healthy lung tissue, remains (3,4). Lung injury caused by MV is also known as ventilator-induced lung injury (VILI), and is characterized by tissue disruption, lung edema and pulmonary inflammation (5–7). A recent study from the ARDS Network indicated that VILI may increase the mortality of critically ill patients by ~10% (8). However, in addition to the protective pulmonary ventilation strategy, there are few specific treatments for VILI, as the potential mechanism remains poorly understood. Therefore, it is imperative to identify the pathogenesis of VILI to provide an experimental foundation for clinical treatment.

In the transcriptome, numerous protein-coding mRNAs, as well as non-protein coding transcripts, are closely associated with the pathogenesis of several diseases. Long non-coding RNAs (lncRNAs) are defined as non-protein coding transcripts >200 nucleotides long (9). In recent years, a growing body of evidence has indicated that lncRNAs are involved in various important biological processes, including proliferation (10), apoptosis (11), brain development (12,13), inflammation (14,15) and immunoregulation (16,17), leading to numerous diseases (18–20); lncRNAs may therefore provide alternative targets for the treatment of these diseases (19,21).

Although the crucial roles of lncRNAs have been investigated in numerous diseases, to the best of the authors' knowledge, no study has been conducted to investigate the profile and function of lncRNAs in VILI. In the present study, differentially expressed lncRNAs and mRNAs were investigated by high-throughput sequencing, and bioinformatics analysis was employed to predict their functions in VILI.

Materials and methods

Animals

Experiments were performed on 16 male ICR mice (age, 7–9 weeks, 25–30 g; Shanghai SLAC Laboratory Animal Co., Shanghai, China). Mice were housed in a temperature-controlled room (22±2°C) under a 12-h light/dark cycle, and were provided ad libitum access to food and water. The humidity was maintained at 50–60%. All experimental protocols and animal handling procedures were approved by the Ethics Committee on Experimental Animals of Shanghai Jiaotong University School of Medicine (Shanghai, China). Mice were divided into two groups: Sham-operated group and VILI group (n=8/group). A total of eight pairs of right lungs were used for RNA-sequencing (RNA-Seq) analysis, whereas eight pairs of left lungs were used for reverse transcription-quantitative polymerase chain reaction (RT-qPCR).

Establishment of a VILI mouse model

Following anesthesia with 100 mg/kg ketamine and 8 mg/kg xylazine (intraperitoneal), mice were fixed in a supine position, endotracheal intubation was performed and the tube was connected to a rodent ventilator (Inspira; Harvard Apparatus Ltd., Holliston, MA, USA). Ventilation was performed at a tidal volume of 30 ml/kg and a respiratory rate of 70 breaths/min for 4 h as previously described (22,23). Mice in the sham group underwent intubation but were allowed to breathe freely. Gas with fraction of inspired oxygen at 21% was used, and the inspiratory/expiratory ratio was 2:1. At the end of the experiment, animals were sacrificed by anesthetic overdose and the lung tissues were harvested for further investigation.

RNA-Seq analysis

Total RNA was extracted from lung tissues in the VILI and sham groups with TRIzol® (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA), according to the manufacturer's protocol. RNA concentration was determined using NanoDrop 2000 (NanoDrop; Thermo Fisher Scientific, Inc., Wilmington, DE, USA), and RNA quality was evaluated with Bioanalyzer 2200 (Agilent Technologies, Inc., Santa Clara, CA, USA) and by 1% agarose gel electrophoresis. RNA with RNA integrity number >8.0 was used to construct a cDNA library. A total of five pair of lungs met the criterion and were used for RNA-Seq analysis.

The cDNA libraries were subsequently constructed with Ion Total RNA-Seq kit version 2.0 (Thermo Fisher Scientific, Inc.). Subsequently, proton sequencing was performed using the cDNA libraries and Ion PI Sequencing 200 kit version 2.0 (Thermo Fisher Scientific, Inc.), according to the manufacturer's protocol. Briefly, samples were mixed and processed on the Ion OneTouch 2 System (Thermo Fisher Scientific, Inc.). Subsequently, they were enriched on the OneTouch 2 ES station for the preparation of template-positive Ion PI Ion Sphere Particles. Finally, the mixed samples were loaded onto 1 P1v2 Proton Chip (Thermo Fisher Scientific, Inc.) for sequencing. Data analysis was performed by NovelBio Bio-Pharm Technology Co., Ltd. (Shanghai, China).

RNA-seq mapping and identification of differentially expressed genes

The Mapsplice v2.2.0 program was used for RNA-seq data mapping (24). Per kilobase per million was used to measure the expression level of each gene (25). The expression of different genes between two groups was evaluated with the EB-Seq (26) algorithm. When the absolute fold change (FC) was >2 and the false discovery rate (FDR) was <0.05, the gene was defined as differentially expressed (26).

Gene ontology (GO) and pathway analysis

GO analysis was performed to analyze the biological functions of differentially expressed genes (27). GO annotations were obtained from gene in NCBI (https://www.ncbi.nlm.nih.gov/gene/) and the Gene Ontology (http://www.geneontoloy.org/). Kyoto Encyclopedia of Genes and Genomics (KEGG) database (http://www.genome.jp.kegg/) was employed to analyze the pathways associated with these differentially expressed genes (28). In addition, differentially expressed genes enriched in ≥2 biological pathways were selected to build the ‘path-act network’ a graphical representation of pathways of the interaction of these pathway terms using cytoscape 3.5.1 (29).

Gene co-expression network analysis

Gene co-expression network analysis, analyzed by R software 3.3.3, was employed to predict the key lncRNAs that function in the important biological processes and signal transduction pathways based on the correlation between differentially expressed lncRNAs and mRNAs (30). The significant correlation pairs (Pearson's correlation coefficient >0.95) were selected to build the network. The core regulatory genes were determined by k-core scoring in the network (31). The k-core indicates the hub or nodal status of a gene having interactions with other genes in the network. The higher the k-core score of a gene, the more central the gene location in the gene co-expression network.

RT-qPCR analysis

To confirm the veracity of RNA-Seq analysis, the expression of three lncRNAs and three mRNAs was measured by RT-qPCR using the SYBR-Green method (FastStart Universal SYBR Master; Roche Diagnostics, Basel, Switzerland) and a MiniOpticon real-time PCR detection system (Bio-Rad Laboratories, Inc., Hercules, CA, USA). Total RNA was extracted from lung tissues using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) and cDNA was synthesized (RevertAid First Strand cDNA Synthesis kit; Thermo Fisher Scientific, Inc.) with the reaction conditions: 25°C 10 min, 42°C 1 h, 72°C 10 min and 4°C hold. All of the primers (Table I) used for RT-qPCR were designed and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). PCR was performed under the following conditions: Predenaturation for 2 min at 95°C with cycle, denaturation for 15 sec at 95°C, annealing for 15 sec at 60°C and extension for 1 min at 72°C with 40 cycles from denaturation to extension. The expression of each gene was calculated using the 2−ΔΔCq (32) method and normalized to that of β-actin. The experiments were repeated twice.

Table I.

Primers used for the amplification of various genes (long non-coding RNAs and mRNAs) in mouse lung tissues.

Description	Sequence (5′-3′)	Product size (bp)	Accession no.
Map2k3os, F	AGCAAAGCAACAGCCTCACT	1,267	NR_027800.1
Map2k3os, R	CACGGGCTCTCTGTGCTTAT
Dnm3os, F	GCCTGGCTGGACAGAGTTGT	7,928	NR_002870.2
Dnm3os, R	TCAATGGCTGGTGGTCATTC
Abhd11os, F	CAGTCACCAGGCCTTGACTC	519	NR_026688.1
Abhd11os, R	CGCTTCTTAGCAATGGCTTC
Gadd45a, F	CTGCAGAGCAGAAGACCGAA	1,224	NM_007836.1
Gadd45a, R	GGGTCTACGTTGAGCAGCTT
Cldn4, F	GGCGTCTATGGGACTACAGG	1,827	NM_009903.2
Cldn4, R	GAGCGCACAACTCAGGATG
Tbxa2r, F	GCTCATCTACCTGCGTGTGG	1,809	NM_001277265.1
Tbxa2r, R	CAGCCTGGAGCTGTGAACTG
β-actin, F	CTGTATGCCTCTGGTCGTAC	214	NM_007393.3
β-actin, R	TGATGTCACGCACGATTTCC

Abhd11os, abhydrolase domain containing 11, opposite strand; Cldn4, claudin 4; Dnm3os, dynamin 3, opposite strand; F, forward; Gadd45a, growth arrest and DNA damage-inducible α; Map2k3os, mitogen-activated protein kinase kinase 3, opposite strand; R, reverse; Tbxa2r, thromboxane A2 receptor.

Cell culture and treatment

To determine the function of differentially expressed lncRNAs, the effects of Map2k3os were measured on the release of inflammatory cytokines from alveolar epithelial cells induced by cyclic stretch. The MLE12 [Obio Technology (Shanghai) Corp., Ltd., Shanghai, China] murine lung epithelial cell line was cultured with Dulbecco's modified Eagle's medium (HyClone; GE Healthcare Life Sciences, Logan, UT, USA) supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) at 37°C in an atmosphere containing 5% CO2 (33). Map2k3os small interfering RNA (siRNA) was synthesized by Shanghai GenePharma Co., Ltd. (Shanghai, China), as was the negative control siRNA. The sequence of mouse Map2k3os siRNA was: 5′-CCCAUGGGAAGCAAAGCAATT-3′. The sequence of negative control siRNA was: 5′-UUCUCCGAACGUGUCACGUTT-3′. MLE-12 cells were seeded onto Collagen I-coated Bioflex 6-well culture plates (BFI products, Inc., Hollywood, FL, USA). When the cells reached 80% coverage of one well, they were transfected with siRNA (60 pmol/ml) using Xfect siRNA transfection reagents (Takara Bio, Inc., Otsu, Japan). Following transfection for 24 h, the medium was renewed and the cells were subjected to cyclic stretch (20% linear elongation, sinusoidal wave, 30 cycles/min) for 4 h, as previously described (34), using the Flexcell-FX-5,000 Tension system (Flexcell International Corporation, Burlington, NC, USA). Control cells were placed in the same culture plates and incubated next to the stretched cells at 37°C in an atmosphere containing 5% CO2. Following cyclic stretch, the cells were harvested to analyze the expression of interleukin (IL)-1β, IL-6 and tumor necrosis factor (TNF)-α by RT-qPCR conducted as previously described in the RT-qPCR subsection. The primer sequences (Sangon Biotech Co., Ltd.) were as follows: IL-1β forward, 5′-TCGCAGCAGCACATCAACAAGAG-3′ and reverse, 5′-TGCTCATGTCCTCATCCTGGAAGG-3′; IL-6 forward, 5′-ACTTCCATCCAGTTGCCTTCTTGG-3′ and reverse, 5′-TTAAGCCTCCGACTTGTGAAGTGG-3′; and TNF-α forward, 5′-GCGACGTGGAACTGGCAGAAG-3′ and reverse, 5′-GCCACAAGCAGGAATGAGAAGAGG-3′.

ELISA

Following cyclic stretch, the conditioned medium was collected. IL-1β (cat. no. F10770), IL-6 (cat. no. F10830) and TNF-α (cat. no. F11630) concentrations were determined using ELISA kits (Shanghai Westang Biotechnology Co., Ltd., Shanghai, China) for mice, according to the manufacturer's protocols.

Statistical analysis

All data were analyzed with SPSS version 20.0 (IBM Corp., Armonk, NY, USA). Student's t-test was applied to identify the differently expressed lncRNAs and mRNAs following RNA-Seq and RT-qPCR. According to previous studies (35,36), Fisher's exact test and χ2 test were used to examine the significance of GO category and pathway analysis. FDR was used to correct the multiple testing to minimize the error of P-value (37,38), and in RNA-Seq to identify the differentially expressed genes (FDR<0.05). One-way analysis of variance followed by Dunnett's test for multiple comparisons was use to compare the difference between multiple groups. P<0.05 was considered to indicate a statistically significant difference.

Results

Differentially expressed lncRNAs and mRNAs

LncRNA and mRNA expression between the VILI and sham groups was analyzed using the DE-Seq algorithm, and differentially expressed genes were identified. A total of 104 lncRNAs and 809 mRNAs were identified as differentially expressed between the groups by high-throughput sequencing. Among the 104 lncRNAs, 74 were upregulated and 30 were downregulated. Of the 809 mRNAs, 521 were upregulated and 288 were downregulated. Volcano plot was used to identify differentially expressed genes (lncRNAs and mRNAs) between the two groups; when the expression was similar, the genes were marked blue, whereas the differentially expressed genes (FC>2 and FDR<0.05) were marked red (Fig. 1). The top 20 differentially expressed genes are presented in Table II (10 upregulated and 10 downregulated lncRNAs). To validate the results from RNA-Seq, the expression levels of three lncRNAs (Map2k3os, Dnm3os and Abhd11os) and three mRNAs (Gadd45a, Cldn4 and Tbxa2r) were detected by RT-qPCR. The results demonstrated that the change in the expression of these lncRNAs and mRNAs determined by RT-qPCR was similar to that from RNA-Seq (Fig. 2A and B).

Figure 1.

Differentially expressed genes between sham and ventilator-induced lung injury groups. The volcano plot includes all of the genes in the two groups. Differentially expressed genes are marked as red dots (FDR<0.05), other genes are marked as blue or green dots. Blue dots meant FDR≥0.05, green dots meant −2≤FC≤2. FC, fold change; FDR, false discovery rate.

Table II.

Top 20 differentially expressed lncRNAs (10 upregulated and 10 downregulated lncRNAs) between ventilator-induced lung injury group and sham group.

AccID	FDR	Regulation	Location	Description
Gm4610	9.71×10−27	Up	chr3	Predicted gene 4610
LOC102636324	1.55×10−16	Up	chr13	Uncharacterized LOC102636324
Gm11827	2.74×10−15	Up	chr4	Predicted gene 11827
LOC101056163	1.01×10−12	Up	chr8	14-3-3 protein sigma-like
Gm5582	3.60×10−12	Up	chr6	Predicted gene 5582
Serpina3h	5.07×10−12	Up	chr12	Serpina3h protein
Krt8-ps	9.79×10−10	Up	chr7	Keratin 8, pseudogene
Gm9523	4.19×10−9	Up	chr5	Methionine adenosyltransferase II, alpha pseudogene
LOC102633421	9.00×10−9	Up	chr11	14-3-3 protein sigma-like
Mir22hg	3.83×10−8	Up	chr11	Mir22 host gene (non-protein coding)
Gm16119	3.98×10−8	Down	chr2	Uridine-cytidine kinase 1-like 1, opposite strand
LOC102632993	6.82×10−8	Down	chr7	Uncharacterized LOC102632993
LOC102635638	4.02×10−7	Down	chr9	Uncharacterized LOC102635638
9630028I04Rik	8.16×10−7	Down	chr17	RIKEN cDNA 9630028I04 gene
A330009N23Rik	4.05×10−6	Down	chr15	RIKEN cDNA A330009N23 gene
AI197445	3.76×10−5	Down	chr13	Expressed sequence AI197445
D930048N14Rik	1.05×10−4	Down	chr11	Protein D930048N14Rik
Panct2	1.57×10−4	Down	chr1	Pluripotency-associated noncoding transcript 2
5930430L01Rik	3.56×10−4	Down	chr5	Katanin p60 ATPase-containing subunit A-like 1
2610507I01Rik	5.23×10−4	Down	chr11	RIKEN cDNA 2610507I01 gene

FDR, false discovery rate.

Figure 2.

(A) Three lncRNAs (Map2k3os, Dnm3os and Abhd11os) identified by RNA-Seq were verified by RT-qPCR. **P<0.01 vs. the sham group. (B) Three mRNAs (Gadd45a, Cldn4 and Tbxa2r) identified by RNA-Seq were verified by RT-qPCR. **P<0.01 vs. the sham group. Abhd11os, abhydrolase domain containing 11, opposite strand; Cldn4, claudin 4; Dnm3os, dynamin 3, opposite strand; Gadd45a, growth arrest and DNA damage-inducible α; lncRNA, long non-coding RNA; Map2k3os, mitogen-activated protein kinase kinase 3, opposite strand; RNA-Seq, RNA-sequencing; RT-qPCR, reverse transcription-quantitative polymerase chain reaction; Tbxa2r, thromboxane A2 receptor.

GO analysis

GO analysis indicated that the differentially expressed mRNAs were associated with numerous important biological processes, cellular components and molecular functions. The present study indicated that 189 GO terms (147 upregulated and 42 downregulated) associated with biological processes were enriched (P<0.01, FDR<0.05). The primary upregulated GO categories were mainly involved in cell apoptosis, neutrophil chemotaxis, inflammatory response, cell proliferation and response to glucocorticoid (Fig. 3A), whereas the downregulated categories were mainly associated with the positive regulation of transcription (DNA-templated), cell-cell signaling, cellular response to interferon-γ, cell cycle and positive regulation of cytosolic calcium concentration (Fig. 3B).

Figure 3.

GO analysis based on differential genes between sham and ventilator-induced lung injury groups. (A) Top 15 GO categories for upregulated genes. (B) Top 15 GO categories for downregulated genes (FDR<0.05). FDR, false discovery rate; GO, gene ontology.

Pathway enrichment analysis

The pathway enrichment analysis was employed to investigate the pathways associated with important differentially expressed genes. -Lg P-value was used to describe the significance level of pathway enrichment. There were 60 upregulated pathways and 21 downregulated pathways. The top 20 differentially expressed pathways with upregulated expression were mainly associated with the TNF, hypoxia-inducible factor (HIF)-1, p53 and Forkhead box O signaling pathways, and those with downregulated expression were associated with neuroactive ligand-receptor interaction, Hippo signaling pathway, phenylalanine metabolism and calcium signaling pathway (Fig. 4A). The differentially expressed genes in these pathways provide evidence for further investigations of the potential mechanisms underlying the pathogenesis of VILI. The interactions between these pathways are illustrated in Fig. 4B.

Figure 4.

Pathway analysis of differentially expressed genes based on the KEGG database. (A) Top 20 signaling pathways, as determined by KEGG. (B) Interactions between the differential signaling pathways (path-act network, a graphical representation of pathways of the interaction of these pathway terms using cytoscape 3.5.1). Red dots represent upregulated pathways, green dots represent downregulated pathways and yellow dots represent bidirectional pathways. KEGG, Kyoto Encyclopedia of Genes and Genomics.

LncRNA-mRNA co-expression network

LncRNA-mRNA co-expression network was constructed based on the differentially expressed genes detected between the VILI and sham group lung tissues. There were 182 network nodes and 509 connections between genes in the VILI group. In this co-expression network, 272 pairs had positive association and 237 pairs exhibited negative association (Fig. 5). The properties of the constructed gene networks were determined by k-cores and the degree of differences. In the co-expression network, lncRNAs with differential degree ≥10 and k-core ≥8 included RIKEN cDNA 4933407K13 gene (4933407K13Rik), RIKEN cDNA E230013L22 gene (E230013L22Rik), RIKEN cDNA 3010003L21 gene (3010003L21Rik), RIKEN cDNA A330009N23 gene (A330009N23Rik), protein D930048N14Rik (D930048N14Rik), predicted gene 11827 (Gm11827) and predicted gene 10872 (Gm10872). The top 20 differentially expressed lncRNAs, as determined by the gene co-expression analysis, are presented in Table III. These differentially expressed lncRNAs may be responsible for the potential mechanism underlying VILI.

Figure 5.

Co-expression network of differential genes in the ventilator-induced lung injury group.

Table III.

Top 20 differentially expressed long non-coding RNAs in the gene co-expression network.

AccID	DifKcore	DifDegree	Description
D930048N14Rik	10	19	Protein D930048N14Rik
Gm10872	9	14	predicted gene 10872
E230013L22Rik	8	19	RIKEN cDNA E230013L22 gene
Gm11827	8	13	predicted gene 11827
4933407K13Rik	8	19	RIKEN cDNA 4933407K13 gene
3010003L21Rik	8	17	RIKEN cDNA 3010003L21 gene
A330009N23Rik	8	17	RIKEN cDNA A330009N23 gene
B130024G19Rik	7	14	Putative uncharacterized protein
Panct2	7	10	Pluripotency-associated noncoding transcript 2
Chd3os	4	4	Putative uncharacterized protein RNF21
AI197445	4	4	Expressed sequence AI197445
Ywhaq-ps	3	6	Tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, theta polypeptide, pseudogene 1
9630028I04Rik	3	6	RIKEN cDNA 9630028I04 gene
LOC102635838	2	6	Uncharacterized LOC102635838
Snhg12	2	3	Small nucleolar RNA host gene 12
Map2k3os	2	6	Mitogen-activated protein kinase kinase 3, opposite strand
4930481A15Rik	2	2	MCG1045479, isoform CRA_b
1810006J02Rik	1	3	RIKEN cDNA 1810006J02 gene
LOC102632993	1	1	Uncharacterized LOC102632993
LOC102637368	1	1	Uncharacterized LOC102637368

DifDegree, differential degree; DifKcore, differential key core.

Effects of Map2k3os siRNA on inflammatory cytokine release from MLE12 cells

Compared with in the control cells, it was identified that Map2k3os siRNA decreased the expression of Map2k3os by 89%, as determined by RT-qPCR (Table IV). Cyclic stretch significantly increased the mRNA expression levels of TNF-α, IL-1β and IL-6 in MLE12 cells, and protein levels in the supernatants. Conversely, Map2k3os siRNA transfection inhibited the expression of these stretch-induced cytokines (Fig. 6). These results indicated the effects of Map2k3os on epithelial secretion of cytokines in cyclic stretch-induced VILI.

Table IV.

Effects of siRNA on long non-coding RNA Map2k3os expression in MLE12 cells (n=4).

Group	Relative expression	Relative reduction (%)
Control siRNA	1±0.184	0
Map2k3os siRNA	0.106±0.042[a]	89

Map2k3os, mitogen-activated protein kinase kinase 3, opposite strand; siRNA, small interfering RNA. Data is presented as the mean ± standard error.

P<0.01 compared with control siRNA group.

Figure 6.

mRNA expression levels of (A) IL-1β, (B) IL-6 and (C) TNF-α in MLE12 cells, as determined by RT-qPCR. Protein level of (D) IL-1β, (E) IL-6 and (F) TNF-α in the culture medium, as determined by ELISA. (n=6). **P<0.01 compared with control group without cyclic stretch or siRNA. ##P<0.01 compared with control siRNA group. CS, cyclic stretch; IL, interleukin; Map2k3os, mitogen-activated protein kinase kinase 3, opposite strand; RT-qPCR, reverse transcription-quantitative polymerase chain reaction; siRNA, small interfering RNA; TNF, tumor necrosis factor.

Discussion

In the present study, RNA-Seq analysis identified 104 lncRNAs and 809 mRNAs, which were differentially expressed in lung tissues between the VILI group and the sham group. Bioinformatics analysis predicted several lncRNAs that were likely to alter the transcription of related mRNAs, and thus participate in the development of VILI. To the best of the authors' knowledge, the present study is the first to analyze the function of lncRNAs in VILI using high-throughput sequencing and bioinformatics analysis.

At present, there have been few studies regarding the transcriptomics of VILI (39–41). These previous studies used microarray analysis, and to the best of the authors' knowledge, a previous study tested the lncRNAs. Ma et al (39) focused on early stress response genes in rodents (mouse and rat) following VILI, and identified 41 upregulated genes and 7 downregulated genes in the injured lungs. GO analysis revealed that these genes may be involved in cell cycle arrest, immune response, inflammatory response, blood coagulation, chemotaxis, apoptosis, cell-cell signaling and negative regulation of cell proliferation (39). Another study (40) only reported inflammasome-related gene expression, including IL-1α, caspase-activator domain-10, and IL-1 receptor-1 and −2. Furthermore, a previous study (41) used microarray analysis to identify the effects of non-muscle myosin light chain kinase isoform on VILI. This study also identified alterations in the expression of several genes in VILI, including thioredoxin domain-containing 9, C-X-C motif chemokine ligand 2 and MYB proto-oncogene, transcription factor. As in previous studies, all of the biological processes reported and the majority of genes identified were also observed in the present study. The difference in the differentially expressed genes among available studies is possibly ascribed to the difference in species (mouse and rat vs. mouse) and methodology (microarray vs. high-throughput sequencing). In the present study, not only mRNAs, but also lncRNAs, were analyzed in the lung tissues, and GO analysis, pathway enrichment analysis, and lncRNA-mRNA co-expression networks were employed to predict the function of lncRNAs. Therefore, the findings may provide evidence regarding the mechanism underlying VILI and aid the development of novel targets to treat VILI.

The onset of VILI is associated with several biological processes and signaling pathways. MV can damage the lungs by direct mechanical injury and indirectly via the induction of biotrauma. Several signaling pathways and mRNAs identified in the present study have been identified in previous studies. Hyperinflation of mouse lungs in vivo and in vitro can induce the expression of nuclear factor (NF)-κB and IL-6 in alveolar macrophages and alveolar epithelial type II cells (42). In addition, chemokines [chemokine (C-X-C motif) ligand 1, cytokine-induced neutrophil chemoattractant-2α and macrophage inflammatory protein 2 (MIP-2)] and their receptors are activated during lung injury (43). The neutrophil elastase inhibitor sivelestat is able to inhibit neutrophil accumulation, and decrease the release of MIP-2, IL-6 and TNF-α in VILI mice (44). The present study also confirmed the important role of inflammatory response, neutrophil chemotaxis and the NF-κB signaling pathway. In addition, apoptosis was identified as the first key biological process in VILI by high-throughput sequencing. Previous studies also identified that pathological cyclic stretch could induce pulmonary epithelial apoptosis and barrier dysfunction in alveolar epithelial cells (45,46). Conversely, a previous study reported that ventilation at low tidal volumes results in a mild inflammatory response without increasing apoptosis, which is accompanied by accelerated alveolar epithelial cell proliferation (47). This differs from the results of the present study, which may be ascribed to the low tidal volume used in the previous study.

KEGG pathway analysis in the present study also demonstrated that several signaling pathways were associated with the pathogenesis of VILI. TNF is able to regulate pulmonary alveolar permeability, alveolar fluid clearance, adhesion molecule expression and leucocyte recruitment (48). TNF has been identified as a marker of early activation of inflammation and can be rapidly released in response to cyclic stretch (49,50). In addition, HIF-1 can promote transcription of A2B adenosine receptor (A2BAR), whereas deletion of the A2BAR gene could reduce survival time and increase pulmonary albumin leakage following VILI (51,52). Furthermore, although the specific role of the P53 signaling pathway in VILI has not been reported, microarray assay has revealed it is involved in the pathogenesis of VILI (41). In the present study, a number of genes were differentially expressed between the sham group and the VILI group. Some genes have been verified by previous studies (49,53), whereas others may be novel targets with which to explore the mechanisms of VILI.

In the lncRNA-mRNA co-expression network, 4933407K13Rik, E230013L22Rik, 3010003L21Rik, A330009N23Rik, D930048N14Rik, Gm11827 and Gm10872 were most likely to be responsible for the potential mechanism underlying VILI. These lncRNAs may directly and indirectly regulate the transcription of related mRNAs, thus participating in the development of VILI. The top 20 differentially expressed lncRNAs are main target genes for future study. The present study confirmed that transfection with siRNA targeting Map2k3os (one of the top 20 lncRNAs) inhibited stretch-induced cytokine (TNF-α, IL-1β and IL-6) expression in MLE12 cells. However, further experiments are required to study the detailed function and mechanism of these lncRNAs in VILI. In addition, next-generation sequencing technology itself has certain limitations, including flux, which means the number of RNA samples that can be measured by a sequenator at the same time, is not high enough and the length of RNA reads remains relatively short (54), which may be improved with the development of improved technology.

In conclusion, the present study provided novel information regarding the lncRNAs, signaling pathways and the co-expression network involved in VILI. Numerous significant lncRNAs are likely to participate in the pathogenesis of VILI, which may be useful to guide further research into the mechanisms and targeted therapy for this disease.