A549 and MRC-5 cell lines were purchased from Cell Bank of Chinese Academy of Sciences (Shanghai, China). Cells were maintained in Dulbecco's modified Eagle's medium and advanced minimum essential medium, respectively, supplemented with 10% newborn calf serum, 100 U/ml penicillin, and 100 µg/ml streptomycin at 37 °C under a humidified atmosphere of 5% CO₂ and 95% air. A549 or MRC-5 cells were collected after treatment with 10 ng/ml TGF-β1 for 12, 24, and 48 h as our previously described (28–30).

Total RNA was harvested and quantified, and its integrity was verified by denaturing gel electrophoresis. Samples with a 28S:18S ratios of approximately 2:1 were accepted for further analysis, which was performed on an Affymetrix Human Transcriptome Array 2.0 microarray for mRNAs and lncRNAs and Affymetrix Gene Chip miRNA 4.0 microarray for miRNAs. The detection of mRNAs, miRNAs, and lncRNAs was performed according to the manufacturer's instructions. Finally, differentially expressed mRNAs, miRNAs, and lncRNAs were identified through fold-change filtering after the quantile normalization of raw data.

Gene ontology (GO) analysis was applied to predict the main functions of the target genes according to the GO project. Fisher's exact test and χ² were used to classify the GO category. Moreover, false discovery rate (FDR) was calculated to correct the P-value. The smaller the FDR, the smaller the error in judging the P-value. The standard difference screening was P<0.01.

Pathway analysis was performed to determine the significant pathway of the differential genes according to the Kyoto Encyclopedia of Genes and Genomes (KEGG), Biocarta, and Reatome. Fisher's exact test and χ² results were used to select the significant pathway, and the significance threshold was defined by P-value and FDR. The standard difference screening was P<0.05.

MiRNA-target gene network was established on the basis of GO and KEGG predicted data for the illustration of the relationship between miRNAs and their target genes.

Total RNA was isolated with TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA quantity and quality were measured on a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity was assessed by standard denaturing agarose gel electrophoresis. Complementary DNA synthesis was performed with a Moloney Murine Leukemia Virus reverse transcriptase kit (Invitrogen) following the manufacturer's instructions. qRT-PCR was performed with a SYBR-Green-based PCR master mix kit (Takara, Shiga, Japan) on a Rotor Gene 3000 real-time PCR system from Corbett Research (Sydney, Australia). Primers of lncRNAs as the following: linc00941 sense: 5′-GCGGTAGCCTTCTCTGAACTG-3′, antisense: 5′-GTTGCATAACCTGACCTGCC-3′; AF086191 sense: 5′-GCAGAGTGGAGCCTTCTCAT-3′, antisense: 5′-TATGCAAACTCCCATGTGGC-3′. Quantification of miR-500a-5p, miR-628-3p, miR-128-3p, miR-223-3p and miR-30a was performed with a stem-loop real-time PCR miRNA kit bought from Guangzhou RiboBio Co., Ltd. (Guangzhou, China).

Protein concentration was quantified using a bicinchoninic acid protein assay kit and boiled with the sample buffer in a water bath for 5 min. Protein samples were separated with 15% SDS-PAGE gels for 2 h and transferred onto a polyvinylidene difluoride membrane, which was subsequently blocked in 5% non-fat milk for 2 h. Blots were probed using the primary antibodies. The anti-vimentin, a-SMA, E-cadherin and Snail antibody was from ProteinTech Group. After three times washing with tris buffered saline tween, the horseradish peroxidase-conjugated secondary antibodies were added. Antigen-antibody complexes were visualized by enhanced chemiluminescence.

Cells were fixed by treatment with fresh 2.5% glutaraldehyde at 4°C for at least 4 h, post-fixed in 1% osmium tetroxide for 1.5 h, dehydrated in gradient ethanol, infiltrated with Epon812, embedded and cultured at 37°C for 12 h, 45°C for 12 h and 60°C for 24 h. Ultrathin sections prepared with an ultracut E ultramicrotome were stained with uranyl acetate and lead citrate and observed using a JEM-1400 TEM (JEOL Ltd., Tokyo, Japan).

All data are expressed as the means ± SD. Statistical analysis was performed using SPSS13.0 software by one-way ANOVA. P<0.05 was considered to indicate a statistically significant difference.

EMT is the phenomenon in which epithelial cells obtain mesenchymal characteristics, including change in shape, expression of mesenchymal markers, such as collagen and a smooth muscle actin (a-SMA); loss of epithelial cell marker E-cadherin; and increase of transcription repressor Snail. Under an optical microscope, the shapes of the cells changed from irregular to spindle. In the TEM images, the number of muscle fibers in the cytoplasm of TGF-β1-activated cells are larger than those of normal control cells (Fig. 1A). Meanwhile, the expression levels of collagen I, a-SMA, E-cadherin, and Snail were significantly increased, whereas that of E-cadherin expression significantly decreased (Fig. 1B-F). These findings indicated the successful establishment of the cell model.

The expressed profiles were determined by using Affymetrix Human Transcriptome Array 2.0 microarray for mRNAs and lncRNAs and Affymetrix Gene Chip miRNA4.0 microarray for miRNAs in A549 and TGF-β1-activated A549 at 12, 24, and 48 h, respectively. Hierarchical clustering analysis results showed that the mRNAs, miRNAs, and lncRNAs of the cell models are differentially expressed (Fig. 2A-C). Basing on our previous studies and other published works, we analyzed the data of TGF-β1-activated A549 only at 24 and 48 h. A total of 242 mRNAs were differentially coexpressed in TGF-β1-activated 24 and 48 h cell samples compared with the normal control group. Of these mRNAs, 139 were upregulated, 103 were downregulated, and 143 miRNAs were differentially coexpressed in TGF-β1-activated 24 and 48 h cell samples compared with the normal control group. Among these miRNAs, 83 were upregulated, 60 were downregulated, and 127 lncRNAs were differentially coexpressed in TGF-β1-activated 24 and 48 h cell samples compared with the normal control group. Among which, 71 were upregulated and 56 were downregulated.

To validate the microarray data, we randomly selected two miRNAs (miR-500a-5p and miR-628-3p) and two lncRNAs (linc00941 and AF086191) and determined their expression levels through qRT-PCR. These RNAs were selected because they have relatively higher differences than other RNAs. As shown in Fig. 2, the two miRNAs and two lncRNAs were differentially expressed in A549 cell models compared with those in controls (Fig. 2D and E). To further validate this result, we tested these ncRNAs on MRC-5 cells derived from human embryonic lung fibroblast cells (Fig. 2F and G). All results are in agreement with those in A549 cell models (12). As expected, we also identified previously discovered fibrotic miR-30a (28,29), thus validating our microarray data.

The functional classification of 242 coexpressed mRNAs in TGF-β1-activated 24 and 48 h cell samples was analyzed using GO analysis. According to the GO annotation tool, differentially upregulated coexpressed mRNAs were principally enriched for GO terms related to collagen catabolism, cell-matrix adhesion, protein amino acid dephosphorylation, cytoskeleton organization and biogenesis, cell-substrate junction assembly, cell differentiation, response to hypoxia, DNA replication, substrate-bound cell migration, and protein amino acid O-linked glycosylation (Fig. 3A). At the same time, we compared the downregulated coexpression mRNAs in TGF-β1-activated 24 and 48 h cell samples with that in the normal control group, which were principally enriched for GO terms related to ion transport, signal transduction, electron transport chain, proteolysis, oxidation reduction, induction of apoptosis by extracellular signals, inflammatory response, Jun-N-terminal kinase cascade, protein kinase cascade, and apoptosis (Fig. 3B). The TGF-β1 induction in EMT-related genes in pulmonary epithelial cells contributes to the upregulation or downregulation of the coexpression of these genes by endowing cells with a specific adhesive, motile, and invasive potential.

Biological pathways associated with coexpressed mRNAs in TGF-β1-activated 24 and 48 h cell samples were analyzed using the KEGG database (Fig. 3C and D). The 10 upregulated and 10 downregulated pathways listed were identified and associated with focal adhesion, extracellular cell matrix (ECM)-receptor nitration, regulation of actin cytoskeleton, hematopoietic cell lineage, gap junction, galactose metabolism, autophagy regulation, p53 signaling pathway, DNA polymerase, tight junction, and cell adhesion molecules.

Furthermore, we further analyzed the interaction network of the 242 coexpressed mRNAs on the basis of protein-protein interaction. In total, 21 upregulated mRNAs and nine downregulated mRNAs were included in the interaction network (Fig. 3E), which showed four types of interaction, namely, activated, inhibited, complementary, and binding interaction, among the coexpressed 30 mRNAs. CTNNA, ITGA11, and ITGB3 are the network nodes or connections in the interaction network. In addition, MMP1, MMP2, and COLLA1, which are the physiological indicators of fibrogenesis during EMT progression, were also highly and positively expressed.

The 30 coexpressed mRNAs were further analyzed. Among these mRNAs, 24 mRNAs, which can coexpress and interact with 11 miRNAs and 33 lncRNAs, exist. The coexpressed interaction network indicated that one mRNA might interact with several miRNAs or lncRNAs. Similarly one miRNA or lncRNA can interact with several mRNAs or lncRNAs (Fig. 4). The degree of mRNA-miRNA-lncRNA coexpression were analyzed. An increased degree indicated an increased coexpression in the mRNA-miRNA-lncRNA interaction network. The coexpressed degree of the mRNAs, miRNAs, and lncRNAs in the network nodes and connections were higher than that of the others (Fig. 5A-C). We further analyzed the lncRNA category. The data showed that most of differentially coexpressed lncRNAs are transcribed from the intergenic gene. Some are from sense-overlapping gene. Meanwhile, a few other unknown sources of lncRNAs exist (Fig. 5D).

In the interaction network, miR-223-3p displays the highest coexpression degree compared with other miRNAs. miR-30a has relatively high coexpression degree, whereas miR-128-3p has relatively low coexpression degree (Fig. 5A). Thus, we selected the three miRNAs and determined their expression levels by using qRT-PCR to illuminate the interaction network. As shown in Fig. 5E, the three miRNAs were differentially decreased in the A549 cell models compared with those in controls. To further validate the data, we tested these miRNA expression levels on MRC-5 cells. The expression trend of these miRNAs is similar pattern to those of miRNAs in A549 cell models (Fig. 5F). Then, analysis of the Pearson's correlation coefficient demonstrated that miR-30a is inversely correlated with collagen I, and miR-128-3p and miR-223-3p are inversely correlated with α-SMA (Fig. 5G-I).

Some researchers named the interaction among mRNA, miRNA, and lncRNA functions ceRNA mechanism. A detailed study should be carried out for determining whether the ceRNA mechanism is involved in the interactions among coexpressed mRNAs, miRNAs, and lncRNAs.

The ceRNA mechanism proposed that all types of RNA transcripts communicate through a new ‘language’ mediated by miRNA-binding sites. The miRNAs bind to sequences with partial complementarity on target RNA transcripts resulting in the repression of target gene expression (17). Thus, we further analyzed the miRNAs and their targeted genes in the interaction network of mRNAs, miRNAs, and lncRNAs. First, the target genes of coexpressed miRNA were predicted based on TargetScan, miRanda data, and miRbase. Subsequently, we analyzed the functional classification of these target genes by GO analysis. The results showed that the upregulated target genes were principally enriched for GO terms related to organic substance response, skeletal system development, transferase activity regulation, mesenchyme development, cell proliferation, and cell morphogenesis (Fig. 6A). The downregulated coexpressed miRNA target genes were principally enriched for GO terms related to the positive regulation of development, negative regulation of cell differentiation, positive regulation of molecular function, and cell fate commitment (Fig. 6B). Most of these upregulated or downregulated miRNA target genes are associated with EMT.

In the interaction network of coexpressed miRNAs and their targeted genes, miR-223-3p, miR-376C-3p, miR-491-5p, miR-505-3p, miR-140-5p, miR-551A, miR-21, miR-1306, miR-128-3p, and miR-382-5p from the mRNA-miRNA-lncRNA interaction network were considered as examples (Fig. 6C). One miRNA had many target genes, whereas different miRNAs can have the same target gene. One gene also can have different binding sites with miRNAs. In total, 57 mRNAs might be the 10 miRNAs' targeted genes. However, among these targeted genes, not a single one is the same as the mRNA in mRNA-miRNA-lncRNA interaction network. This finding means that the ceRNA mechanism does not exist in the above coexpressed mRNA-miRNA-lncRNA interaction network. Analysis with DIANA-miRPath determined 10 candidate miRNAs, which were involved in possible disease-related pathways (P-value cutoff at 0.05). miR-223-3p, miR-491-5p, and miR-140-5p were significantly implicated in disease-related pathways and were found to be associated in 12, 10, and eight pathways, respectively. The Venn diagram revealed the intersection (Fig. 6D). Thus, cancer and pancreatic cancer pathways are crossing pathways related to miR-223-3p, miR-491-5p, and miR-140-5p.

In addition, our previous study demonstrated that miR-30a may function as a novel therapeutic target for lung fibrosis by blocking mitochondrial fission dependent on dynamin-related protein1 (Drp-1). miR-30a exerts its regulatory function through regulating Drp-1 promoter hydroxymethylation by directly targeting TET1. To date, our further detailed study has showed no existence of ceRNA mechanism in miR-30a regulatory mode (28,29).