Identification of transforming growth factor-beta-regulated microRNAs and the microRNA-targetomes in primary lung fibroblasts

Background Lung fibroblasts are involved in extracellular matrix homeostasis, which is mainly regulated by transforming growth factor-beta (TGF-β), and are therefore crucial in lung tissue repair and remodeling. Abnormal repair and remodeling has been observed in lung diseases like COPD. As miRNA levels can be influenced by TGF-β, we hypothesized that TGF-β influences miRNA expression in lung fibroblasts, thereby affecting their function. Materials and methods We investigated TGF-β1-induced miRNA expression changes in 9 control primary parenchymal lung fibroblasts using miRNA arrays. TGF-β1-induced miRNA expression changes were validated and replicated in an independent set of lung fibroblasts composted of 10 controls and 15 COPD patients using qRT-PCR. Ago2-immunoprecipitation followed by mRNA expression profiling was used to identify the miRNA-targetomes of unstimulated and TGF-β1-stimulated primary lung fibroblasts (n = 2). The genes affected by TGF-β1-modulated miRNAs were identified by comparing the miRNA targetomes of unstimulated and TGF-β1-stimulated fibroblasts. Results Twenty-nine miRNAs were significantly differentially expressed after TGF-β1 stimulation (FDR<0.05). The TGF-β1-induced miR-455-3p and miR-21-3p expression changes were validated and replicated, with in addition, lower miR-455-3p levels in COPD (p<0.05). We identified 964 and 945 genes in the miRNA-targetomes of unstimulated and TGF-β1-stimulated lung fibroblasts, respectively. The TGF-β and Wnt pathways were significantly enriched among the Ago2-IP enriched and predicted targets of miR-455-3p and miR-21-3p. The miR-455-3p target genes HN1, NGF, STRADB, DLD and ANO3 and the miR-21-3p target genes HHEX, CHORDC1 and ZBTB49 were consistently more enriched after TGF-β1 stimulation. Conclusion Two miRNAs, miR-455-3p and miR-21-3p, were induced by TGF-β1 in lung fibroblasts. The significant Ago2-IP enrichment of targets of these miRNAs related to the TGF-β and/or Wnt pathways (NGF, DLD, HHEX) in TGF-β1-stimulated fibroblasts suggest a role for these miRNAs in lung diseases by affecting lung fibroblast function.


Introduction
Lung fibroblasts play a key role in extracellular matrix (ECM) homeostasis and maintenance of the normal lung architecture and are therefore crucial players in lung damage and repair.They have been postulated as important players in lung diseases with disturbed ECM homeostasis and aberrant repair such as chronic obstructive pulmonary disease (COPD).Transforming growth factor beta (TGF-β), an important cytokine in tissue repair and remodeling, stimulates ECM production by fibroblasts mainly through signaling via the downstream SMAD proteins [1][2][3].In addition to stimulation of ECM production, this multifunctional cytokine regulates ECM homeostasis by influencing the expression and activity of matrixmetalloproteinases and their inhibitors [4].
MicroRNAs (miRNAs) have been shown to play a key role in the regulation of cellular activity.MiRNAs are small non-coding RNAs with an average length of 22 nucleotides.They are estimated to regulate over 60% of the protein-coding genes [5].MiRNAs regulate gene expression by partially binding to complementary homologous sequences present on their target mRNA transcripts which then leads to either inhibition of protein translation or mRNA degradation.
The TGF-β signaling pathway is involved in the maturation of a subset of miRNAs, indicating that the levels of several miRNAs might be influenced by TGF-β [6].SMAD proteins are able to mediate the maturation of miRNAs [7] and can indirectly regulate transcription of miRNAs with a SMAD-binding element in their promoter [8,9].Conversely, several miRNAs have been reported to influence TGF-β signaling by regulating the expression of TGF-β1, TGF-β receptors 1 and 2, SMAD1-5 and SMAD7 [6,10,11].
COPD is one of the lung diseases with, mostly smoking-induced, disturbed repair capacity of the lung, resulting in emphysematous lung tissue destruction and (small) airway wall fibrosis [12].TGF-β levels are higher in COPD patients than in control subjects, underscoring the importance of the TGF-β signaling pathway in COPD [12,13].
Several miRNAs have been found to be deregulated in lung tissue, serum or sputum of COPD patients (reviewed by Osei et al. [14]).The expression level of miR-146a was decreased in IL-1β/TNF-α-treated primary lung fibroblasts of COPD patients compared to those of control subjects [15].Moreover, altered levels of circulating miRNAs have been proposed as potential biomarkers for COPD [16].
To our knowledge, only limited information is available regarding the interplay between miRNAs and the TGF-β-induced repair response in primary lung fibroblasts, nor is there information on the actual target gene repertoire of miRNAs (i.e. the miRNA-targetome) in these cells.Therefore, the aim of the current study was to investigate the TGF-β1-induced miRNA expression changes in primary lung fibroblasts and to identify the transcripts that are likely to be affected by these miRNAs using the lung fibroblast miRNA-targetomes in unstimulated and TGF-β1-stimulated primary lung fibroblasts.Furthermore a comparison was made between control and COPD fibroblasts, being one of the lung diseases in which ECM remodeling plays an essential role.

Materials and methods Subjects
Parenchymal lung fibroblasts from 9 current or ex-smoking control subjects, undergoing lung surgery for tumor resection were subjected to miRNA profiling using a microarray approach.These primary fibroblasts were isolated from parenchymal lung tissue that was located far away from the tumor and lacked abnormalities as checked by histology on haematoxylin and eosin stained slides [17,18].The replication study group consisted of parenchymal fibroblasts from 10 control subjects and 15 COPD patients and was used to replicate the array findings and to investigate the miRNA expression in parenchymal fibroblasts of COPD patients.The controls had no history of lung disease, other than the lung tumor.COPD patients with a history of another lung disease such as asthma or interstitial lung disease were excluded.
This study was conducted according to national ethical and professional guidelines on the use human body material ("Code of conduct; Dutch federation of biomedical scientific societies"; https://www.federa.org/codes-conduct)and the Research Code of the University Medical Center Groningen (https://www.umcg.nl/EN/Research/Researchers/General/ResearchCode/Paginas/default.aspx).
Lung fibroblasts used in this study are derived from left-over lung material after lung surgery and transplant procedures.Currently, this material is not subject to the act on medical research involving human subjects in the Netherlands and therefore an ethics waiver was provided by the Medical Ethical Committee of the University Medical Center Groningen (METc UMCG).All samples and clinical information were de-identified before experiments were performed.

Primary lung fibroblast culture
Primary parenchymal fibroblasts were isolated, cultured and stored in liquid nitrogen until further use as described previously [19].The fibroblasts were cultured with Ham's F12 medium supplemented with 10% (v/v) fetal calf serum (FCS), 100 U/ml penicillin/streptomycin and 200 mM L-glutamine (all from Lonza, Breda, The Netherlands) at 37˚C in a 90% (v/v) humidified atmosphere with 5% (v/v) CO 2 .The experiments were performed on these fibroblasts at passage 5.For the initial array analysis control fibroblasts in the discovery group were stimulated with 100 U/ml TGF-β1 (equivalent to 3.1 ng/ml) (R&D Systems, Abingdon, UK) in complete Ham's F12 medium containing 0.5% (v/v) FCS for 24 h.In the replication group, fibroblasts were stimulated with 2.5 ng/ml and 7.5 ng/ml TGF-β1 for 24 h.

RNA isolation
Total RNA was isolated from primary parenchymal lung fibroblasts using TRIzol (Invitrogen, Carlsbad, CA, USA), according to the protocol of the manufacturer.RNA in the total (T) fraction of the immunoprecipitation of argonaute-2 (Ago2-IP) was isolated using miRNeasy Mini Kit (Qiagen, Venlo, The Netherlands), whereas RNA from the IP fraction was isolated using miRNeasy Micro Kit (Qiagen) according to manufacturer's protocol.The RNA concentration was measured with a NanoDrop 1000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA).

MiRNA expression profiling
RNA samples were hybridized with the Human GeneChip miRNA 1.0 array (Affymetrix), containing 847 probe sets for human miRNAs.Robust Multi-array Average was performed using GeneSpring GX version 13.1.1software (Agilent Technologies, Santa Clara, CA, USA) for the quantile normalization of the probe-level intensity measurements.MiRNAs were filtered based on the following criteria: expression values in the range of 75 th -100 th percentile in at least 50% of all samples (9 out of 18).After filtering, 205 miRNAs were left for further analyses.Statistical analyses in GeneSpring GX software for miRNA expression profiling were performed using the paired samples t-test and Benjamini-Hochberg false discovery rate (FDR) to correct for multiple testing.An FDR of p<0.05 was considered statistically significant.Unsupervised hierarchical clustering with Pearson's correlation was performed to generate a heatmap of differentially expressed miRNAs using Genesis software version 1.7.6 (Graz University of Technology, Graz, Austria) [20].All differentially expressed miRNAs with a fold change (FC) !1.5 and a normalized signal intensity value of at least 75 in unstimulated or TGF-β1-stimulated fibroblasts were selected for validation.
Next, qRT-PCR was performed in qPCR MasterMix Plus (Eurogentec) and TaqMan micro-RNA assay (Applied Biosystems).Small nucleolar RNA, C/D box 48 (RNU48) was used as reference gene.Samples including the no template control of both cDNA synthesis and qRT-PCR were run in triplicate and the formula 2 -ΔCp was used to calculate the relative miRNA expression levels.The data were analyzed using LightCycler1480 software release 1.5.0 (Roche Diagnostics GmbH).

Immunoprecipitation of Ago2-RISC complex
To identify the miRNA-targetome of primary parenchymal lung fibroblasts of control subjects, Ago2-RIP-Chip was performed as described previously [22,23].We used primary parenchymal lung fibroblasts from two ex-smoking controls (control 1 and 2) with and without TGF-β1 (7.5 ng/ml) stimulation.In brief, 17-20 million cells were harvested and lysed using polysome lysis buffer.Cell lysates were incubated overnight with monoclonal mouse anti-human Ago2 (Clone 2E12-1C9, Abnova, Taipei City, Taiwan) coated Sepharose G beads (Abcam, Cambridge, UK).Lysates generated from the same number of cells were incubated with normal mouse polyclonal IgG 1 isotype control (12-371, Merck Millipore, Amsterdam, The Netherlands) coated beads which served as a negative control for the IP procedure.To monitor the efficiency of the Ago2-IP, protein samples of the total (T), flow through (FT) and IP fractions were loaded on a 7.5% (w/v) polyacrylamide gel.Western blotting was performed as described previously [22].

mRNA expression profiling and identification of the miRNA-targetomes
To determine which transcripts are enriched in the IP fraction after Ago2-IP, mRNA expression profiling was performed using G3 Human Gene Expression 8x60K v3 Microarrays (Agilent Technologies).RNA of the T and IP fraction (50-80 ng) was labeled with cyanine 3 (Cy3) and cyanine 5 (Cy5) using the Two-Color Low Input Quick Amp Labeling Kit (Agilent Technologies) according to manufacturer's protocol Version 6.9.1.NanoDrop 1000 Spectrophotometer (Thermo Scientific) was used to quantify the cRNA yield and to determine the specific activity.
Dye-swap hybridizations were performed overnight using Gene Expression Hybridization Kit (Agilent Technologies) and slides were scanned with the Agilent SureScan Microarray Scanner (Agilent Technologies) using the two-color gene expression protocol GE2_1200_Jun14.Feature Extraction software version 12.0.1.1 was used to extract information from probe features from our microarray scan data.The data was analyzed using GeneSpring GX version 13.1.1software (Agilent Technologies).Only probes with consistent values in the dye-swap experiments (0.5>Cy3/Cy5<2) were used for further analyses.Out of the 58,341 probes on the microarray, 52,843 showed consistent values in the IP fraction and 51,588 in the T fraction of control 1, of which 47,917 were overlapping between the IP and T fractions.In the TGF-β1-stimulated fibroblasts of control 1, 51,991 probes in the IP fraction and 52,210 probes in the T fraction were consistent between the dyes of which 47,655 probes were overlapping between the fractions.For control 2, 50,963 and 49,118 probes were consistent for the IP and T fractions, respectively, and 44,551 probes were in overlap between the two fractions.In the TGF-β1-stimulated fibroblasts of control 2, 52,011 probes in the IP fraction and 50,654 probes in the T fraction were consistent between the dyes of which 46,404 probes were overlapping between the fractions.For each of the four IP-experiments, the list of probes with consistent signal intensities between the two dyes was further filtered using the flagged present status (detectable above background according to the feature extracting software).All probes flagged present in at least one out of two conditions (IP and T) were included in the further analyses (37-43% of probes discarded of the consistent list).Finally, the lists of the two unstimulated samples were overlapped (24,502 remaining) and similarly, the two TGF-β1-stimulated samples were overlapped (25,939 remaining).Both the IP/T ratios of the two unstimulated fibroblast samples and those of the two TGF-β1-stimulated fibroblasts samples were statistically correlated using the Pearson correlation test.The miRNA-targetomes were defined as the top 1,500 of the most IP-enriched transcripts in the fibroblast samples of both control subjects.To identify which of the miR-455-3p and miR-21-3p target genes are influenced by TGF-β1 stimulation, we ranked the Ago2-IP enriched and predicted targets by the IP/T ratio in each of the four IP experiments (Control 1 and 2, both with and without TGF-β1).For miR-455-3p that is broadly conserved, only the predicted targets with conserved binding sites were considered in these analyses (TargetScan version 7.1 [24]).MiR-21-3p is poorly conserved; therefore the predicted miR-21-3p targets with a cumulative weighted context ++ score <-0.2 were included.The target genes that are affected by the TGF-β1-dependent modulation of miR-455-3p or miR-21-3p are those genes that are more prominently IP-enriched in TGF-β1-stimulated fibroblasts compared to unstimulated fibroblasts.This was defined according to the following criteria 1) predicted target gene is present in the top 1,500 most enriched transcripts in the TGF-β1-stimulated fibroblasts; 2) predicted target gene is more enriched in the TGF-β1-stimulated sample with a difference in ranking of !300; 3) the more pronounced enrichment in the TGF-β1-stimulated fibroblasts is consistent in both controls.

Gene ontology and pathway analysis
Gene set enrichment analysis (GSEA software v2.2.2) was performed using the Molecular Signatures Database (MSigDB v5.1) to check for enrichment of miRNA target genes sets and other gene ontologies [25].The probes were ranked based on the mean IP/T ratio of the two control subjects.
GeneNetwork analysis was performed to identify enrichment of biological processes and pathways (KEGG and reactome) amongst the TargetScan (version 7.1) predicted and Ago2-IP enriched targets of miR-455-3p and miR-21-3p [24].GeneNetwork is based on an independent gene expression dataset and predicts (currently unknown) gene functions based on coexpression.This information is used to assign genes to specific pathways and biological processes [1,26].

Statistical analyses
IBM SPSS Statistics 22 software was used to compare differences in subject characteristics between control subjects and COPD patients, and between the two study groups using Mann-Whitney U test.Differences in ECM gene, α-SMA gene and miRNA expression after TGF-β1 stimulation were analyzed using one-sided paired Wilcoxon signed rank tests.The Mann-Whitney U test was used to compare expression of genes and miRNAs between controls and COPD patients.Chi-square test was used to test differences in the percentage of predicted targets in the Ago2-IP fraction for miR-455-3p and miR-21-3p in the miRNA-targetome compared to the percentage of predicted targets in all expressed genes.A p-value below 0.05 was considered statistically significant.

Subject characteristics
The clinical characteristics of both study groups are shown in Table 1.The discovery group consisted of 4 ex-and 5 current smoking controls.The replication group consisted of 4 exand 6 current smoking controls and 15 COPD patients.Both age and number of pack-years did not differ between both study groups, or between controls and COPD patients.Furthermore, the Forced Expiratory Volume in one second (FEV 1 ) and the Forced Expiratory Volume in one second/Forced Vital Capacity (FEV 1 /FVC) ratio in controls did not differ between discovery and replication groups.

Differentially expressed miRNAs after TGF-β1 stimulation
In the discovery group, TGF-β1 stimulation caused a significant induction of the known TGFβ-regulated genes FN1, COL1A1, and α-SMA (p<0.01,Fig 1A) confirming successful TGF-β1 stimulation.We found 29 miRNAs to be differentially expressed after TGF-β1 stimulation in the discovery group (FDR<0.05,Fig 2A).Out of these 29 miRNAs, 8 miRNAs were downregulated and 21 miRNAs were upregulated.Fourteen out of 29 miRNAs had a fold change of at least 1.5, of which 12 showed increased expression after TGF-β1 stimulation (Table 2).Based on the signal intensity, 6 of these 14 miRNAs were selected for qRT-PCR validation and replication (marked in bold in Table 2).

Replication of TGF-β1 effects on miRNA expression in COPD and control fibroblasts
In the replication group, we included primary parenchymal lung fibroblasts from control subjects as well as from COPD patients.Again, the TGF-β1 stimulation of 2.5 ng/ml (p<0.01,Fig 1B) and 7.5 ng/ml was successful as it significantly induced the mRNA expression of FN1, COL1A1, and α-SMA in the primary parenchymal lung fibroblasts of both controls and COPD patients (p<0.01,Fig 1C).
In the primary lung fibroblasts of control subjects, we were able to replicate the expression changes of 2 out of 5 miRNAs (miR-455-3p and miR-21-3p) after 2.5 ng/ml TGF-β1 (p<0.01,Fig 3A) as well as after 7.5 ng/ml TGF-β1 stimulation (p<0.01,Fig 3B).In the fibroblasts of COPD patients, we found significant upregulation of miR-455-3p and miR-21-3p after 2.5 ng/ ml as well as after 7.5 ng/ml TGF-β1 (p<0.0001).In addition, we found a significant downregulation of miR-490-3p in COPD fibroblasts after stimulation with 2.5 ng/ml (p<0.01,Fig 3A) and a significant induction of miR-143-3p in COPD fibroblasts after stimulation with 7.5 ng/ ml TGF-β1 (p<0.01,Fig 3B), which both were not observed in the replication of the control subjects.
As a secondary aim we also compared the expression levels of the 5 selected miRNAs in primary parenchymal lung fibroblasts between control subjects and COPD patients.MiR-455-3p levels were significantly lower in COPD patients (p<0.05) both at basal levels and after TGF-β1 stimulation (p<0.05,Fig 3A and 3B).In addition, a lower miR-143-3p expression was

Identification of the miRNA-targetomes in primary lung fibroblasts
To define the mRNAs that are actively targeted by miRNAs in primary parenchymal lung fibroblasts we profiled the mRNAs recovered from Ago2-IP from unstimulated and TGF-β1-stimulated primary parenchymal lung fibroblasts of two control subjects.The efficiency of the IP procedure was confirmed by Western blotting of Ago2, showing a clear enrichment of the Ago2 protein in the IP fraction of the Ago2-IP sample and not in the IgG1 control (Fig 4A).Additional qRT-PCR of three miRNAs also showed a clear increase in miRNA levels in the Ago2-IP fractions (Fig 4B).The Ago2-IP results were consistent between the two controls, as the IP/T ratios of the 24,502 probes corresponding to 17,129 unique genes in the unstimulated lung fibroblasts and the 25,939 probes corresponding to 18,037 unique genes in the TGF-β1-stimulated lung fibroblasts were highly correlated between the two subjects (R 2 = 0.8164, p<0.0001 and R 2 = 0.7257, p<0.0001, respectively).GSEA analysis using the ranking of the IP/ T ratio of all genes that are consistently expressed in unstimulated lung fibroblasts of the two subjects demonstrated an enrichment of 17 miRNA-target gene sets in the top 30 enriched gene sets (Table 3).The top 5 most significant miRNA gene sets that are identified in primary lung fibroblasts were the miRNA (seed families of) Let-7, miR-29, miR-26, miR-181 and miR-101.Similar GSEA analysis was also performed on all consistently expressed genes in the TGF-β1-stimulated lung fibroblasts of the two subjects and showed 15 miRNA-target gene sets in the top 30 enriched gene sets (Table 3).In these TGF-β1-stimulated lung fibroblasts, the most significant active miRNA gene sets were the miRNA (seed families of) Let-7, miR-26, miR-202, miR-27 and miR-29.Together these analyses support the robustness of our IP approach.
The top 1,500 most IP-enriched probes that were consistently enriched in the fibroblasts of the two control subjects are defined as the parenchymal lung fibroblast "miRNA-targetome".Of the top 1,500 most IP-enriched probes in the unstimulated primary parenchymal lung fibroblasts of the two control subjects, we identified 1,110 probes overlapping between the two control subjects, which corresponded to 964 unique genes, and these were defined as the miRNA-targetome of unstimulated cells (Fig 5A).Of the top 1,500 most IP-enriched probes in the TGF-β1-stimulated lung fibroblasts of the two control subjects, we identified 1,083 probes overlapping between the control subjects, which corresponded to 945 unique genes, and these were identified as the miRNA-targetome of the TGF-β1-stimulated cells (Fig 5B).The IP/T ratios of the top 1,500 most IP-enriched probes in the unstimulated and TGF-β1-stimulated lung fibroblasts of the two control subjects and the identified genes in the miRNA-targetomes are shown in S1 Tables.
MiR-455-3p and miR-21-3p predicted targets are enriched in the lung fibroblast targetomes As the TGF-β1 effect on miR-455-3p and miR-21-3p was consistent throughout the experiments, we decided to focus on these two miRNAs.We determined which transcripts might be affected by the TGF-β1-induced changes in miR-455-3p and miR-21-3p expression by identifying TargetScan predicted targets of these two miRNAs within the miRNA-targetome lists.
In the unstimulated as well as in the TGF-β1-stimulated fibroblasts of both subjects, a significant increase of genes with at least one conserved miR-455-3p binding site was observed in the top 1,500 most IP-enriched genes, compared to the predicted miR-455-3p targets amongst all expressed genes (p<0.0001,Fig 6).In total, 57 and 50 predicted miR-455-3p target genes were enriched in the miRNA-targetome of unstimulated and TGF-β1-stimulated lung fibroblasts, respectively.Also for miR-21-3p a significant increase in target genes was observed in the top 1,500 most enriched genes (p<0.0001,Fig 6).For miR-21-3p, 55 and 54 predicted target genes were enriched in the miRNA-targetome of unstimulated and TGF-β1-stimulated lung fibroblasts, respectively.Identification of the processes and pathways of the miR-455-3p and miR-21-3p target genes To investigate the possible functions of miR-455-3p and miR-21-3p in the unstimulated and TGF-β1-stimulated lung fibroblasts, we subjected the predicted and IP-enriched targets of these miRNAs to GeneNetwork analysis, focusing on Biological processes, Kegg pathways and Reactome gene sets.Per collection, the top 10 most significant biological processes and pathways of the 57 and 50 predicted miR-455-3p targets of unstimulated and TGF-β1-stimulated fibroblasts are depicted in Table 4. Several biological processes and pathways are related to very general processes such as gene expression, differentiation, proliferation, various types of cancer and signaling pathways.Interestingly, the TGF-β signaling pathway (p 5x10 -6 ) and the Wnt signaling pathway (p 7x10 -6 ) were significantly enriched among these genes both in the unstimulated and in the TGF-β1-stimulated lung fibroblasts.Furthermore, the TGF-β receptor signaling pathway was significantly enriched in the TGF-β1-stimulated fibroblasts (p = 2x10 -8 ).Several of the predicted miR-455-3p targets were allocated to the TGF-β-related processes and pathway and to the Wnt signaling pathway by GeneNetwork as annotated or unannotated genes (Table 5).Amongst the 55 and 54 predicted miR-21-3p targets present in the miRNA-targetome of unstimulated and TGF-β1-stimulated fibroblasts, respectively, the top 10 most significant biological processes and pathways per collection, are depicted in Table 6.These processes and pathways are related to kinase activity, cancer, cell cycle, metabolism and various signaling pathways.The TGF-β signaling pathway was significantly enriched among the genes with a predicted miR-21-3p binding site both in unstimulated (p 2x10 -3 ) and TGF-β1-stimulated lung fibroblasts (p 2x10 -3 ).Moreover, among the predicted miR-21-3p target genes in the unstimulated lung fibroblasts, the Wnt signaling pathway was also significantly enriched (p = 4x10 -3 ).The predicted miR-21-3p targets involved in TGF-β and Wnt signaling pathways were allocated as annotated genes by GeneNetwork (Table 5).

Identification of target genes modulated by TGF-β1-induced miR-455-3p and miR-21-3p
Next, we investigated which target genes are affected by stimulation with TGF-β1.These target genes are defined as being more prominently enriched in the IP of TGF-β1-stimulated fibroblasts compared to unstimulated fibroblasts.Five miR-455-3p predicted target genes HN1, NGF, STRADB, DLD and ANO3 showed an increased enrichment in the IP upon TGF-β1 stimulation.Of these genes, NGF and DLD were predicted by GeneNetwork to be involved in the TGF-β1-related biological processes and signaling pathway (marked bold in Table 5).Furthermore, DLD is predicted to play a role in the Wnt signaling pathway.
Of the predicted miR-21-3p target genes, increased Ago2-IP enrichment upon TGF-β1 stimulation was observed for HHEX, CHORDC1 and ZBTB49.HHEX was one of the genes predicted to be involved in the TGF-β signaling pathway (marked bold in Table 5).

Discussion
In the present study, we investigated the effect of TGF-β1 on miRNA expression in primary lung fibroblasts and used the miRNA-targetomes of these cells to identify transcripts that are likely to be affected by the TGF-β1-dependent modulation of miRNA expression levels.Twenty-nine miRNAs were differentially expressed after TGF-β1 stimulation.TGF-β1-induced miR-455-3p and miR-21-3p expression was validated with qRT-PCR and also replicated in an independent set of control and COPD fibroblasts.In the replication set, lower levels of miR-455-3p were observed in COPD as compared to control fibroblasts.The TGF-β signaling pathway was significantly enriched amongst the Ago2-IP enriched predicted targets of miR-455-3p and miR-21-3p and we have identified several target genes, which were more Ago2-IP-enriched upon TGF-β1 stimulation.
We showed that TGF-β1 regulates expression of multiple miRNAs in primary lung fibroblasts of control subjects and COPD patients with consistent and significant induction of miR-455-3p and miR-21-3p.In a previous study of Milosevic et al., 84 miRNAs were shown to be differentially expressed after TGF-β1 stimulation in commercially available human lung fibroblasts [27].Seven out of the 29 differentially expressed miRNAs in our initial miRNA array screening were overlapping with the 84 miRNAs found by Milosevic et al.In accordance with their results, we showed upregulation of miR-21-3p, miR-23a-5p, miR-503-5p, miR-424-3p, and miR-214-5p upon TGF-β1 stimulation.However, they found downregulation of miR-455-3p and miR-23a-3p [27], whereas we found upregulation of these two miRNAs.These inconsistencies might be due to differences in using a single versus multiple donors, culture methods and TGF-β1 stimulation time (2 h versus 24 h).Notably, upregulation of miR-455-3p and miR-23a-3p after TGF-β stimulation has been observed in other cell types [28][29][30][31].Other studies observed increased expression levels of miR-27a-3p and miR-199a-5p upon TGF-β stimulation in MRC-5 lung fibroblasts in accordance with our initial miRNA array results [8,32].
Using Ago2-RIP-Chip, we defined the miRNA-targetomes of both unstimulated and TGF-β1-stimulated primary lung fibroblasts.By assessing which of the miR-455-3p and miR-21-3p predicted targets are present in these targetomes we gained insight in the potential role of these TGF-β1-induced miRNAs in lung fibroblasts.To our knowledge this is the first time that the miRNA-targetomes have been identified in primary lung fibroblasts.The lists of experimentally proven targets of miRNAs in primary lung fibroblasts with and without TGF-β1 stimulation represent a valuable source for future studies.The predicted miR-455-3p and miR-21-3p targets present in the miRNA-targetomes of unstimulated and TGF-β1-stimulated lung fibroblasts were significantly enriched for genes related to the TGF-β processes and signaling pathway.These findings suggest that there is significant cross-talk between miRNAs and the TGF-β pathway.Moreover, the enrichment of the TGF-β signaling pathway suggests a role for these miRNAs in regulating ECM homeostasis by lung fibroblasts [6].
Previous data suggested that the TGF-β1-induced miR-455-3p expression is regulated through a SMAD2/3 binding element and that miR-455-3p directly targets SMAD2, ACVR2B and CHRDL1 [31].This suggests a feedback loop resulting in decreased SMAD2/3 signaling.From these reported targets of miR-455-3p, ACVR2B was found in the miRNA-targetome of unstimulated primary fibroblasts.SMAD2 and CHRDL1 were expressed in primary lung fibroblasts but were not enriched in the IP fractions.Of the IP-enriched and predicted miR-455-3p target genes, the levels of enrichment of DLD and NGF transcripts in the Ago2-IP fraction were increased upon TGF-β1 stimulation.Both genes were predicted by GeneNetwork to play a role in the TGF-β-related processes and signaling pathway.DLD is a sub-protein of the pyruvate dehydrogenase protein complex, which is involved in the cellular metabolism [33].Interestingly, pyruvate dehydrogenase activity was found to be decreased in the muscles of COPD patients [33].NGF is a neurotrophin which can be produced and released by lung fibroblasts during inflammation.A previous study showed that NGF accelerates the migration of human fetal lung fibroblasts to fibronectin indicating that NGF may stimulate repair [34].Furthermore, NGF itself and TGF-β enhance the NGF release by MRC-5 cells and NGF, similarly to TGF-β, induces α-SMA expression in these cells [35].The same study showed that expression of NGF and its receptors, trkA and p75, was increased in fibrotic lung tissue compared to normal tissue [35].Taken together, NGF may have pro-fibrotic effects in the lung and TGF-β1-induced miR-455-3p may serve to control NGF production inhibiting its pro-fibrotic effects on lung fibroblasts.
A second gene set enriched within the Ago2-IP-enriched miR-455-3p targets is the Wnt signaling pathway, with TCF7L1 being a predicted target and present in the miRNA-targetomes of both unstimulated and TGF-β1-stimulated lung fibroblasts.TCF7L1 mediates the transcription of the canonical Wnt signaling target genes [36].DLD was a second gene predicted to be involved in the Wnt signaling pathway.The Wnt pathway regulates the differentiation of lung fibroblasts into myofibroblasts after activation via TGF-β and this enhances the ECM production [37].Therefore, activation of the Wnt signaling pathway has been suggested to be an important factor in tissue repair.Deregulation of this pathway via miRNAs may thus contribute to the disturbed ECM homeostasis in lung diseases.
The lower expression level of miR-455-3p in COPD patients compared to control subjects indicates its potential role in aberrant tissue repair and remodeling of the lungs via effects on the target genes involved in the Wnt and TGF-β signaling pathways.
As mentioned above, the TGF-β signaling pathway was also enriched amongst the IPenriched and predicted miR-21-3p target genes.HHEX is one of the genes whose Ago2-IP enrichment increases upon TGF-β1 stimulation and predicted to be involved in the TGF-β signaling pathway.The expression of this homeobox containing gene was found to be regulated by the BMP signaling pathway and is implied to be involved in lung organogenesis [38].
Other IP-enriched and predicted miR-21-3p target genes that were more enriched after TGF-β1 stimulation are CHORDC1 and ZBTB49.CHORDC1 has been shown to inhibit kinase activity of ROCK2 in mouse embryonic fibroblasts [39].Rho kinases have been implicated in asthma and COPD and have been suggested as a promising therapeutic target for these lung diseases [40].ZBTB49 is a zinc finger that can inhibit cell proliferation of HEK293 and HCT116 cells and thus may affect cell proliferation in lung fibroblasts as well [41].

Conclusions
We demonstrated that TGF-β1 affects the expression of several miRNAs in primary parenchymal lung fibroblasts, including miR-455-3p and miR-21-3p.Furthermore, we showed that the targets of miR-455-3p and miR-21-3p play a role in the TGF-β and Wnt signaling pathways.In addition, our study showed that TGF-β1-dependent miRNA modulation affects several genes that may influence the function of lung fibroblasts via several processes including tissue repair, cellular metabolism, cell migration and cell proliferation.Interestingly, we found a difference in miR-455-3p expression between COPD patients and control subjects, which suggests a role for this miRNA in COPD pathogenesis.These findings add to our understanding of the role of miRNAs in the functionality of lung fibroblasts in tissue repair and remodeling and support a role for aberrant miRNA regulation in lung diseases such as COPD in which high TGF-β levels have been reported.

Fig 5 .
Fig 5. Defining the miRNA-targetomes of unstimulated and TGF-β1-stimulated lung fibroblasts.The overlap of the top 1,500 most IP-enriched probes in the (A) unstimulated and (B) TGF-β1-stimulated fibroblasts of the two control subjects are defined as the miRNA-targetome.The identified genes in the miRNA-targetomes are listed in S1 Tables.https://doi.org/10.1371/journal.pone.0183815.g005

Fig 6 .Table 5 .a
Fig 6.Significant Ago2-IP-enrichment of miR-455-3p and miR-21-3p predicted targets.For each miRNA, the percentages of predicted targets were calculated in all expressed genes and in the top 1,500 most IP-enriched genes in all four IP experiments.Chisquare test was used to determine whether the number of predicted targets in the Ago2-IP fraction for miR-455-3p and miR-21-3p in the top 1,500 most enriched genes was significant different from the expected based on the number of predicted targets in all expressed genes.https://doi.org/10.1371/journal.pone.0183815.g006

Table 1 . Subject characteristics. Subjects Discovery group Replication group Control subjects d Control subjects d COPD patients
FEV 1 , % predicted, percentage of Forced Expiratory Volume in one second of the predicted normal value for an individual of the same sex, age and height.c FEV 1 /FVC, Forced Expiratory Volume in one second/Forced Vital Capacity ratio expressed in percentage, a measurement for obstruction/restriction in the lungs.