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).

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.

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).

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).

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, 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.

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.

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% CO₂ (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% CO₂. 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′.

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.

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 χ² 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.

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).

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).

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.

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.

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.