In Silico Detection and Characterization of microRNAs and Their Target Genes in microRNA Microarray Datasets from Patients with Systemic Sclerosis

Abstract

Interstitial lung disease (ILD) is the main reason of death in patients with systemic sclerosis (SSc). The potential microRNA (miRNA)–messenger RNA (mRNA) interaction networks of SSc-ILD from a systematic biological perspective are unclear. To characterize differentially expressed miRNAs (DE-miRNAs) and differentially expressed genes (DEGs) likely related to SSc-ILD, we downloaded the miRNA microarray dataset (GSE81923) and mRNA datasets (GSE76808 and GSE81292) from the Gene Expression Omnibus database. Comprehensive bioinformatic analyses were conducted to predict target genes for DE-miRNAs and generate an miRNA-hub gene network with SSc-ILD. In total, 26 DE-miRNAs were detected in SSc-ILD, among which 2 were upregulated and 24 were downregulated. Additionally, 178 common DEGs (55 upregulated and 123 downregulated) were identified. miRNAs were primarily enriched in pathways involving inflammation and regulation of fibroblasts. The hub genes identified were MMP7, IER2, HBEGF, CCL4, NFKBIA, JUNB, LIF, SERPINE1, FOSL1, and NAMPT. We discovered the miRNA-mediated regulatory network in SSc-ILD using an integrated bioinformatic analysis. The findings provide novel insight and expand our comprehension of the molecular mechanisms participating in the pathogenesis of SSc-ILD, along with identification of new potential diagnostic biomarkers.

Introduction

Systemic sclerosis (SSc) is an uncommon disease of connective tissues that is specified by excessive collagen deposition, microvascular involvement, and autoimmunity (Servettaz et al., 2006). Pulmonary involvement, particularly interstitial lung disease (ILD), results in high mortality rates in patients with SSc (Elhai et al., 2012). Approximately 79% of SSc patients are found to have ILD during autopsy (D'Angelo et al., 1969). One multicenter observational study of over 5000 patients with SSc showed that ILD accounted for more than one-third of SSc-related deaths (Tyndall et al., 2010). The prevalence of SSc-ILD is high, with more than 15% of patients progressing to severe lung fibrosis. The mechanism and etiology in SSc-ILD formation and progression have been investigated extensively, yet causes of SSc-ILD remain unclear. Distinguishing between individuals with severe disease and patients with slow or stable disease makes it challenging for physicians to identify appropriate therapies (Roth et al., 2011).

microRNAs (miRNAs), known as a class of short noncoding RNAs, have the capability of negatively regulating the expression of genes through interactions with the 3′-untranslated region of the targeted messenger RNA (mRNA) (Wang et al., 2018). miRNAs are recognized for their important roles in numerous cellular processes, such as cellular apoptosis, proliferation, and differentiation (Baumjohann and Ansel, 2013; Wang et al., 2018). To date, miRNAs and their association with pulmonary fibrosis have become a central topic for many biomedical researchers around the globe. miRNA and mRNA show excellent stability and reliability in the blood, which increases their potential use as novel markers for detection and monitoring of SSc-ILD in the clinic (Young et al., 2017). In the past decade, technologies classified as “-omics” have been widely utilized to characterize alternations of gene expression and identify potential biomarkers at the molecular level with detailed insights into the disease (Kim et al., 2006; Milano et al., 2008). SSc-ILD is thought to be associated with multiple putative miRNAs, including miR-29a (Khalil et al., 2015), miR-155 (Christmann et al., 2016), miR-21, miR-182 (Kouri et al., 2015), and miR-92a (Bagnato et al., 2017). The different reports on modified miRNA expression were primarily caused by sample heterogeneity, different sequencing platforms, and distinct analytical approaches. However, the relationship between SSc-ILD and miRNA-mRNA interaction networks remains unclear from a systematic biological perspective.

Microarray datasets were retrieved from the Gene Expression Omnibus (GEO) database to investigate correlations between differentially expressed genes (DEGs) and differentially expressed miRNAs (DE-miRNAs) by bioinformatic analyses. This allowed for characterization of DE-miRNAs in the miRNA expression profiles, along with DEGs in the mRNA expression profiles of patients with SSc-ILD. The DE-miRNA and DEG data were combined for the gene ontology (GO) and pathway enrichment analysis along with functional module analyses of protein–protein interaction (PPI) networks. The interaction network between miRNA and mRNA associated with SSc-ILD provides the fundamentals for expanding our knowledge of biological mechanisms underlying the pathogenesis of SSc-ILD and for identification of diagnostic and therapeutic targets of SSc-ILD.

Materials and Methods

Data collection and inclusion criteria

The National Center for Biotechnology Information (NCBI) GEO database* is a public genomics data repository that contains gene expression datasets and original series and platform records. Microarray data from Idiopathic pulmonary fibrosis (IPF)-relevant miRNA and mRNA expression profiles were obtained and extracted from the GEO database. Keywords, including “Systemic sclerosis,” “interstitial lung disease,” “Homo sapiens,” and “expression profiling by array,” were used to search for publicly available studies, from which 11 GSE studies were retrieved after a systematic review and screened with the inclusion criteria as follows: (1) samples diagnosed with SSc-ILD tissues and normal tissue samples and (2) gene expression profiling of mRNA or miRNA. Finally, the miRNA microarray dataset, GSE81293 (Christmann et al., 2016), and the mRNA microarray datasets, GSE76808 and GSE81292 (Christmann et al., 2014), were used for further analysis. Baseline characteristics and demographic information of patients with SSc-related ILD and healthy controls were previously described by Christmann et al. (2014, 2016).

Microarray data information

The microarray dataset, GSE81293, was from the GPL16384 platform ([miRNA-3] Affymetrix Multispecies miRNA-3 Array) with lung biopsy tissues from 20 matched SSc-ILD and normal lung samples. The dataset, GSE76808, was from the GPL571 platform ([HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array) with samples from seven SSc-ILD patients and four healthy persons. The GSE81292 dataset was from the GPL18991 platform ([HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array), containing 15 lung biopsy samples of SSc-ILD and 5 healthy persons.

Data processing and characterization of DE-miRNAs and DEGs

Data were preprocessed for background adjustment, normalization, and summarization in R language (version 3.5.0). Probe levels were calculated and converted into gene expression levels according to the annotation files in the GEO database. DE-miRNAs and DEGs between the SSc-ILD and normal controls were analyzed with the limma package (version 3.34.8) in the R language (Smyth, 2004). p-Values were adjusted by the Benjamini–Hochberg false discovery rate method, and adjusted p-value (adj. p) <0.05 and |log fold change (FC)| >1 were used as the cutoff values. The heat map of DEGs was generated using the pheatmap function of the R/Bioconductor package (Reimers and Carey, 2006). A volcano plot of DE-miRNAs and DEGs was drawn using the ggplot2 package (version 3.1.0).

Functional and pathway enrichment analysis

GO for DEGs was performed to assign the genes and their products into biological processes (BPs), molecular functions (MFs), and cellular components (CCs). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was utilized with the Database for Annotation, Visualization, and Integrated Discovery (DAVID)^**; to identify functional and metabolic pathways (Huang Da et al., 2009). A bubble chart was generated using the ggplot2 package (version 3.1.0). GO and pathway enrichment analyses of DE-miRNAs were applied by DIANA miRPath (Vlachos et al., 2012), and any p-value <0.05 was rated to be statistically significant.

PPI network and detection of hub genes

The Search Tool for the Retrieval of Interacting Genes (STRING)^** is an online database for discovering gene interactions, which includes both physical and functional associations (Szklarczyk et al., 2011). STRING was used to evaluate interactive associations among DEGs with a confidence score of >0.4 as the cutoff criterion. Cytoscape software (version 3.5.1) (Shannon et al., 2003) was applied to construct the PPI network. The hub genes in the network were detected to characterize key elements in SSc-ILD, using cytoHubba (version 0.1) in Cytoscape (Chin et al., 2014).

Prediction of miRNA targets

To understand the associations between miRNAs and gene expression, the genes targeted by DEGs were predicted from GSE81293 using TargetScan^† (Lewis et al., 2005), miRanda (John et al., 2004), and miRDB^‡ tools (Wong and Wang, 2015). Genes from three prediction methods were considered as the candidate targets. The mRNA-miRNA interactions were generated by Cytoscape to further characterize correlations between miRNAs and potential target genes (Kohl et al., 2011).

Results

Baseline characteristics of patients with SSc-related ILD

All patients with SSc-related ILD from the mRNA microarray datasets, GSE76808 and GSE81292 (Christmann et al., 2014), were female. More than half (52%) of the patients had Diffuse cutaneous SSc (dcSSc) and 48% were positive for the anti-Scl-70 (antitopoisomerase I) antibody. Most patients with SSc-ILD developed progressive disease in the miRNA microarray dataset GSE81293 (Christmann et al., 2016) and showed an average baseline score of 6.88 ± 4.47 compared with 12.76 ± 6.59 after treatment. The detailed characteristics and demographic information of patients with SSc-related ILD and healthy controls were previously described by Christmann et al. (2014, 2016).

Detection of DEGs in two mRNA microarray datasets

As shown in Table 1, a total of 29 SSc-ILD and 9 matched normal lung tissues were analyzed, from which 356 and 413 DEGs were identified from datasets GSE76808 and GSE81292, respectively, using adj. p < 0.05 and |log FC| >1 as the cutoff. Among them, 139 DEGs were upregulated and 217 were downregulated in GSE76808 (Fig. 1a, b), while 165 were upregulated and 248 were downregulated in GSE81292 (Fig. 1c, d). As shown in Table 2, a total of 178 DEGs were identified by the integrated analysis, including 55 upregulated DEGs and 123 downregulated DEGs between the SSc-ILD and healthy controls.

FIG. 1.

Identification of DEGs between SSc-ILD and healthy controls. (a) Volcano plots of DEGs in the GSE76808 dataset. Light gray points in the plot represent upregulated DEGs. Medium gray points in the plot represent downregulated DEGs. (b) Heat map of all DEGs in the GSE76808 dataset. (c) Volcano plots of DEGs in the GSE81292 dataset. Light gray points in the plot represent upregulated DEGs. Medium gray points in the plot represent downregulated DEGs. (d) Heat map of all DEGs in the GSE81292 dataset. DEG, differentially expressed gene; ILD, interstitial lung disease; SSc, systemic sclerosis.

Table 1.

Dataset Details of the microRNA and Gene Expression Profiles Related to Patients with Systemic Sclerosis–Interstitial Lung Disease

	GEO ID	Platform	Organism	Experiment type	Samples (case vs. control)	Country	Year
mRNA	GSE81292	GPL18991	Homo sapiens	Expression profiling by array	15 vs. 5	United States of America	2017
	GSE76808	GPL571	Homo sapiens	Expression profiling by array	14 vs. 4	United States of America	2016
microRNA	GSE81293	GPL16384	Homo sapiens	Noncoding RNA profiling by array	15 vs. 5	United States of America	2016

GEO, Gene Expression Omnibus.

Table 2.

One Hundred Seventy-Eight Differentially Expressed Genes Identified in Patients with Systemic Sclerosis–Interstitial Lung Disease

Regulation	DEGs (gene symbol)
Upregulated	BMP4, CDH3, AQP5, TYMS, BLNK, COL3A1, GGH, CXCL14, CEP55, IGF1, IL13RA2, MEOX2, KIAA0125, IGHM, MMP7, CTSE, GDE1, CCNB1, PDGFD, PTPLAD1, FMO1, ZKSCAN7, ELF5, SERPINF1, ABCC3, RARRES1, PCNA, IRX5, PRDX4, ZBED8, CD27, TNFRSF17, DSG2, AADAC, COL1A1, POSTN, AMPD1, NMU, MEOX1, OLFM1, SEL1L3, CFI, PPP1R3C, IGJ, ENC1, PERP, CD38, EAF2, SLC22A3, GALNT12, IGFBP2, SGPP1, COL15A1, STEAP1, YIPF5
Downregulated	IL7R, MT1G, RGS9, SLCO4A1, RGS2, EDN1, RRAD, SERPINE1, IER2, ARC, CXCL3, IL4R, CYP3A5, PLA1A, S100A3, IDO1, IER3, CEBPD, PPP1R15A, MAP3K8, FUT1, STC1, LRRC32, IRF1, EFCC1, NR4A2, SGK1, NEDD9, SPSB1, ANGPTL4, CLDN5, NAMPT, CSF3, DUSP6, RGS16, EDNRB, LIMS2, FSTL3, CA4, ADM, LINC00312, SPRY4, FAM107A, ADAMTS1, PLK3, ATF3, PTGS2, DUSP1, ITPKC, PIGA, AGER, GFPT2, CXCL2, FOSB, HAS2, RGCC, IL1B, CD83, IL1R2, STC2, CHI3L2, CCL4, NR4A1, KDM6B, BMP2, LAMC3, SERPINB1, HBEGF, NFKBIA, VIPR1, FOS, JUNB, TMEM2, TIPARP, RHOB, RALA, RAMP3, CAMTA1, TRIB1, SLC2A3, C8orf4, PTX3, DUSP5, CD69, ELF3, APOLD1, PER3, C10orf10, CX3CL1, ZNF331, MMP19, IL18R1, ICAM1, NR4A3, JUN, PELI1, MT1X, MT2A, RASIP1, LIF, KLF9, RND1, CYR61, MAFF, TNFAIP3, SIK1, CCL2, FOSL1, MT1M, BHLHE40, IL6, NXF3, IL1RL1, IL18RAP, CXCL8, EGR1, RGS1, SLC19A2, GADD45B, ARID5A, KLF4, ESM1, ATHL1

Regulation

DEGs (gene symbol)

Upregulated

BMP4, CDH3, AQP5, TYMS, BLNK, COL3A1, GGH, CXCL14, CEP55, IGF1, IL13RA2, MEOX2, KIAA0125, IGHM, MMP7, CTSE, GDE1, CCNB1, PDGFD, PTPLAD1, FMO1, ZKSCAN7, ELF5, SERPINF1, ABCC3, RARRES1, PCNA, IRX5, PRDX4, ZBED8, CD27, TNFRSF17, DSG2, AADAC, COL1A1, POSTN, AMPD1, NMU, MEOX1, OLFM1, SEL1L3, CFI, PPP1R3C, IGJ, ENC1, PERP, CD38, EAF2, SLC22A3, GALNT12, IGFBP2, SGPP1, COL15A1, STEAP1, YIPF5

Downregulated

IL7R, MT1G, RGS9, SLCO4A1, RGS2, EDN1, RRAD, SERPINE1, IER2, ARC, CXCL3, IL4R, CYP3A5, PLA1A, S100A3, IDO1, IER3, CEBPD, PPP1R15A, MAP3K8, FUT1, STC1, LRRC32, IRF1, EFCC1, NR4A2, SGK1, NEDD9, SPSB1, ANGPTL4, CLDN5, NAMPT, CSF3, DUSP6, RGS16, EDNRB, LIMS2, FSTL3, CA4, ADM, LINC00312, SPRY4, FAM107A, ADAMTS1, PLK3, ATF3, PTGS2, DUSP1, ITPKC, PIGA, AGER, GFPT2, CXCL2, FOSB, HAS2, RGCC, IL1B, CD83, IL1R2, STC2, CHI3L2, CCL4, NR4A1, KDM6B, BMP2, LAMC3, SERPINB1, HBEGF, NFKBIA, VIPR1, FOS, JUNB, TMEM2, TIPARP, RHOB, RALA, RAMP3, CAMTA1, TRIB1, SLC2A3, C8orf4, PTX3, DUSP5, CD69, ELF3, APOLD1, PER3, C10orf10, CX3CL1, ZNF331, MMP19, IL18R1, ICAM1, NR4A3, JUN, PELI1, MT1X, MT2A, RASIP1, LIF, KLF9, RND1, CYR61, MAFF, TNFAIP3, SIK1, CCL2, FOSL1, MT1M, BHLHE40, IL6, NXF3, IL1RL1, IL18RAP, CXCL8, EGR1, RGS1, SLC19A2, GADD45B, ARID5A, KLF4, ESM1, ATHL1

DEG, differentially expressed gene.

GO and pathway enrichment analysis of DEGs

As shown in Figure 2a, 178 DEGs were annotated and mainly classified into BPs, MFs, and CCs. In the BP category, DEGs were significantly enriched in the organic cyclic compound, regulation of BPs, regulation of cell processes, and cell proliferation. In the CC group, GO terms were higher in the extracellular space, extracellular region, extracellular matrix (ECM), cell surface, and cytoplasm. In the MF category, DEGs were enriched in cytokine activity, transcriptional activator activity, receptor binding, cytokine receptor binding, chemokine activity, core promoter proximal region sequence-specific DNA binding, and RNA polymerase II core promoter proximal region sequence-specific DNA binding. Next, 178 DEGs were enriched in 34 signal pathways, including the tumor necrosis factor (TNF) signaling pathway, cytokine–cytokine receptor interaction, NF-kappa B signaling pathway, NOD-like receptor signaling pathway, Toll-like receptor (TLR) signaling pathway, rheumatoid arthritis, and MAPK signaling pathway. The top 20 pathways are shown in Figure 2b.

FIG. 2.

GO term and KEGG pathway enrichment analysis of 178 DEGs. (a) GO terms of significantly enriched genes for biological processes, molecular functions, and cellular components. (b) Top 20 significantly enriched KEGG pathways. The gradual color represents the adjusted p-value, and the size of the black spots represents the gene numbers. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Construction of the PPI network and hub gene analysis

The PPI network of 178 DEGs, which contained 177 nodes and 522 edges, was built using the Cytoscape software based on the STRING database (Fig. 3a). The hub genes were characterized by the Cytoscape software (Fig. 3b). The 10 top hub nodes with the highest degrees of connectivity were matrix metallopeptidase 7 (MMP7), immediate early response 2 (IER2), heparin-binding EGF-like growth factor (HBEGF), CC motif chemokine ligand 4 (CCL4), NFKB inhibitor α (NFKBIA), JunB proto-oncogene (JUNB), leukemia inhibitory factor (LIF), interleukin 6 family cytokine, serpin family E member 1 (SERPINE1), FOS-like 1, AP-1 transcription factor subunit (FOSL1), and nicotinamide phosphoribosyltransferase (NAMPT).

FIG. 3.

PPI network and hub genes of potential DEGs. (a) PPI network of DEGs. (b) The top 20 hub genes in the PPI network of predicted target genes. PPI, protein–protein interaction.

DE-miRNA pathway enrichment analysis

A total of 26 DE-miRNAs were identified from GSE81293 expression profiling datasets by comparing the SSc-ILD samples with healthy controls. Among them, 2 miRNAs (hsa-let7e and hsa-miR-455) were significantly upregulated, while the other 24 miRNAs were downregulated in the SSc-ILD samples (Table 3). The top 10 dysregulated miRNAs were hsa-miR-4459, hsa-miR-4507, hsa-miR-4417, hsa-miR-3687, hsa-miR-1246, hsa-miR-4668, hsa-miR-4530, hsa-miR-4299, hsa-miR-4516, and hsa-miR-663b. The KEGG analysis showed that 26 DE-miRNAs were enriched in 26 signaling pathways (p < 0.05), of which 20 most significantly enriched pathways, as determined by p-values, are shown in Figure 4a. GO term analysis demonstrated that DE-miRNAs were enriched in 175 functional GO terms, with the 20 most common shown in Figure 4b. A heat map of selected pathways and GO categories of the union genes are shown in Figure 4.

FIG. 4.

Heat map of DE-miRNAs in the GSE81293 dataset. (a) Significant pathways for DE-miRNAs. (b) Significant GO terms for DE-miRNAs. DE, differentially expressed.

Table 3.

Twenty-Six microRNAs Significantly Expressed in Systemic Sclerosis–Interstitial Lung Disease from the GSE81293 Dataset

Regulation	miRNA	log FC	adj. p-value
Upregulated	hsa-let-7e	1.02	0.010
Upregulated	hsa-miR-455	1.25	0.022
Downregulated	hsa-miR-1246	−4.42	0.016
	hsa-miR-3687	−4.35	0.016
	hsa-miR-4417	−3.55	0.011
	hsa-miR-4505	−2.8	0.035
	hsa-miR-4459	−2.75	0.004
	hsa-miR-3663	−2.61	0.030
	hsa-miR-663b	−2.44	0.023
	hsa-miR-4656	−2.41	0.030
	hsa-miR-27a	−2.29	0.049
	hsa-miR-4530	−2.15	0.016
	hsa-miR-4507	−2.15	0.005
	hsa-miR-4440	−2.12	0.035
	hsa-miR-4299	−1.97	0.017
	hsa-miR-4793	−1.90	0.049
	hsa-miR-665	−1.84	0.030
	hsa-miR-4668	−1.54	0.016
	hsa-miR-149	−1.49	0.049
	hsa-miR-762	−1.45	0.025
	hsa-miR-4281	−1.39	0.049
	hsa-miR-4516	−1.31	0.022
	hsa-miR-4463	−1.24	0.030
	hsa-miR-4685	−1.21	0.042
	hsa-miR-1469	−1.13	0.049
	hsa-miR-4508	−1.12	0.049

FC, fold change.

Correlation between DE-miRNAs and DEGs associated with SSc-ILD

The candidate target genes were predicted from the 26 dysregulated miRNAs using three miRNA-target interaction tools, including miRanda/mirSVR, TargetScan, and miRDB. Intersections between the candidate target genes and 178 common DEGs from three datasets were determined. As shown in Table 4, 14 miRNAs and their target genes commonly existed in DEGs of the two mRNA datasets. In total, six DEGs of the hub genes that were commonly predicted from two mRNA datasets were regulated by dysregulated miRNA (Table 5). The C-X3-C motif chemokine ligand 1 (CX3CL1) was predicted as the target of three miRNAs, and Fos-related antigen 1 (FOSL1) and FOSB (known as G0/G1 switch regulatory protein 3, G0S3) as targets of two miRNAs. Cytoscape was used to generate the miRNA-gene network for further elucidation of correlations between miRNAs and potential target genes (Fig. 5).

FIG. 5.

The regulation network of DE-miRNAs and DEGs in SSc-ILD.

Table 4.

Fourteen microRNAs and Their Target Genes Found in Differentially Expressed Genes of the Two Common Messenger RNA Datasets

miRNA	log FC	adj. p-value	Target genes
hsa-miR-4505	−2.8	0.035	KDM6B, LRRC32, NR4A1, FOSB, STC1, FUT1, CX3CL1
hsa-miR-4459	−2.75	0.004	APOLD1, SERPINB1, FAM107A, FOSL1, BMP2, CDH3, SPRY4, ITPKC
hsa-miR-3663	−2.61	0.030	KDM6B, DSG2
hsa-miR-663b	−2.44	0.023	SPRY4, SLC2A3
hsa-miR-4656	−2.41	0.030	EFCC1, ARC
hsa-miR-27a	−2.29	0.049	C8orf4
hsa-miR-4530	−2.15	0.016	CX3CL1, HAS2, TRIB1, JUNB, AQP5
hsa-miR-4507	−2.15	0.005	COL1A1, CTSE, STC1
hsa-miR-4299	−1.97	0.017	POSTN, TIPARP, RND1
hsa-miR-4668	−1.54	0.016	AGER, TNFAIP3, APOLD1, C8orf4, PERP, SGK1, HAS2, CDH3, SPRY4, COL15A1, TRIB1
hsa-miR-762	−1.45	0.025	COL1A1, KDM6B, STC1, CX3CL1, FAM107A, SPSB1, SIK1, ARC, FSTL3
hsa-miR-4516	−1.31	0.022	SEL1L3, GFPT2, CD83, TNFAIP3, LRRC32, FOSB, EDNRB, FAM107A, FOSL1, LIMS2, SPRY4
hsa-miR-4463	−1.24	0.030	ELF3, FAM107A
hsa-miR-4508	−1.12	0.049	KDM6B, CLDN5

Table 5.

Dysfunctional microRNAs Regulated Six Differentially Expressed Genes in the GSE81293 Dataset

Gene name	miRNA	Log FC
Gene name	miRNA	GSE81292	GSE76808
CX3CL1	3	−1.807	−2.026
FOSL1	2	−1.74	−1.815
ELF3	1	−1.04	−1.919
NR4A1	1	−2.367	−1.078
FOSB	2	−4.667	−5.269
JUNB	1	−1.987	−2.028

Discussion

ILD is a common presentation of SSc with high mortality rates, yet the pathogenesis and optimal treatment of SSc-ILD remain unclear. Recently, several miRNAs, including the miR-29 family, miR-30 family, and let-7d (Tanaka et al., 2013; Bagnato et al., 2017), were found to be involved in the pathogenesis of many human diseases, such as cancer and inflammatory diseases (Davidson-Moncada et al., 2010). In addition, these miRNAs were shown to modulate several fibrotic-related genes with multifunctional roles in SSc fibrosis (Angiolilli et al., 2018). The overexpression of miR-29a in SSc fibroblasts reduces the expression of antiapoptotic proteins, including BCL-2, increasing the susceptibility to cell death (Jafarinejad-Farsangi et al., 2015). The inhibition of profibrotic miRNA miR-21, which is commonly overexpressed in SSc, results in increased apoptosis in SSc fibroblasts (Jafarinejad-Farsangi et al., 2016). In addition, miR-30-induced therapies were demonstrated by Tanaka et al. (2013) in patients with IPF and SSc. In the present study, 26 DE-miRNAs (2 upregulated and 24 downregulated miRNAs) were identified from the miRNA expression microarray dataset with SSc-ILD and control lung tissue samples (Table 3) and were most significantly enriched in the signaling pathways associated with inflammation and fibroblast regulation, such as ECM-receptor interaction, transforming growth factor-beta (TGF-β) signaling pathway, Wnt signaling pathway, and cell adherence junctions (Fig. 4).

Any upregulated miRNA typically leads to declined production of the target protein, while downregulation often leads to an excess in production of target proteins (Makino and Jinnin, 2016). Previous studies revealed the effects of miRNAs on fibrogenic activity by targeting TGF-β signaling events in lung epithelial cells and fibroblasts and regulating target genes in lung inflammation and lung fibrosis processes (Lepparanta et al., 2012; Christmann et al., 2014). ECM-receptor interactions are closely related to IPF (Hoffmann et al., 2014). In addition, let-7 regulated epithelial–mesenchymal transition and TGF-β to participate in pulmonary fibrosis (Pandit et al., 2010). MiR-92a decreased TGF-β1-induced Wnt1 inducible signaling pathway protein 1 (WISP1) expression in lung fibroblasts ex vivo (Lappi-Blanco et al., 2013; Berschneider et al., 2014). As cell adhesions and biological adhesions are critical processes for the formation of pulmonary fibrosis (Lappi-Blanco et al., 2013), we hypothesized that DE-miRNAs identified in this study might be associated with progression of SSc-ILD. Future studies are required to functionally and clinically validate DE-miRNAs in progression of SSc-ILD.

Microarray and high-throughput screening technologies have been widely applied to predict new target genes. However, previous studies have focused on a single genetic event or individual cohorts with limited sample sizes. In the present study, miRNA expression microarray data from two cohorts were collected, from which a total of 178 DEGs (55 upregulated and 123 downregulated) were detected. Signaling pathway analysis revealed associations of the cytokine–cytokine receptor interaction and NF-kappa B and TLR signaling pathways with progression of ILD in SSc patients. Innate immune signaling by TLRs has been generally accepted as a crucial factor that drives the persistent fibrotic response in SSc (Bhattacharyya and Varga, 2015). The KEGG analysis and preliminary bioinformatic analysis confirmed the cytokine–cytokine receptor interaction signaling pathway as a critical inflammatory regulator associated with chronically inflamed fibrosis. In addition, many hub genes were characterized by the PPI network construction for the DEGs, with the top 10 being MMP7, IER2, HBEGF, CCL4, NFKBIA, JUNB, LIF, SERPINE1, FOSL1, and NAMPT. Future studies are needed to validate these important hub genes as potential prognostic markers and therapeutic targets for SSc-ILD.

We reported possible joint action of DE-miRNAs and DEGs in SSc-ILD. In the present study, miRNA-mRNA interactions were analyzed using miRNA and mRNA expression profiles to gain additional insight into the essential roles of miRNAs in regulating progression of lung fibrosis in SSc. In total, 14 miRNAs showed regulatory effects on the identified DEGs that were overlapped between the 2 mRNA datasets and 10 DE-miRNAs and 6 DEGs involved in this miRNA-target relationship (Tables 4 and 5; Fig. 5), suggesting possible coregulation of DE-miRNAs and DEGs in SSc-ILD. Among them, CX3CL1 was predicted as the target of three miRNAs, and FOSL1 and FOSB as the target of two miRNAs. CX3CL1 is generated by type 2 pneumocytes and airway epithelium (Marasini et al., 2005) and forms the biological axis (CX3CR1/CX3CL1) with CX3CR1 to recruit plasma cells to the interstitium for generation of autoantibodies that can damage the parenchyma of SSc-ILD lungs. Recently, a large cohort study demonstrated associations of CX3XL1 with antitopoisomerase I antibody and lung fibrosis, and eventually with ILD progression, using multivariable regression analysis (Hoffmann-Vold et al., 2018). Both FOSL1 and FOSB genes encode leucine zipper proteins that are members of the FOS gene family and enable dimerization of proteins in the JUN family, including JUNB, to regulate cell proliferation, differentiation, and transformation. JUN and FOS were positively expressed in fibrous dysplasia (Sakamoto et al., 1999).

Overexpression of FOSL1 increases ILD responding to TLR signaling (Takada et al., 2011). FOSL1 expression in the lungs mediates antifibrotic effects by modulating proinflammatory, profibrotic, and fibrotic gene expression, providing a potential biomarker for pulmonary fibrosis with poor prognosis and treatment options (Rajasekaran et al., 2012). E74-like ETS transcription factor 3 (ELF3) was upregulated by the inflammatory cytokines, interleukin-1β and TNF-α, in bronchial epithelial cell lines (Wu et al., 2008), playing essential functions in regulating epithelial cell differentiation during lung regeneration (Oliver et al., 2012). Nuclear receptor subfamily 4 group A member 1 (NR4A1), a transcription factor expressed in various tissues, was found to promote fibrogenesis by persistent stimulation with TGF-β1 and AKT-dependent phosphorylation (Zeng et al., 2018).

There are a few limitations in the present study. First, the sample size in the microarray dataset was limited. Second, findings of the in silico identification of DEGs and DE-miRNAs and their potential target genes lack experimental validation in vitro and in vivo. Finally, the associations between miRNA and mRNA expression profiles require further experimental validation. Future investigations with larger clinical sample sizes, along with optimized validation procedures, are needed.

Data Availability

Data supporting the results of the present investigation are available from the corresponding author upon reasonable request.

Footnotes

Disclosure Statement

No competing financial interests exist.

Funding Information

This study was financially supported by the National Natural Science Foundation of China (No. 81602004 and No. 81600259).

References

Angiolilli

, Marut

, van der Kroef

, Chouri

, Reedquist

K.A.

, and Radstake

T.R.D.J.

(2018). New insights into the genetics and epigenetics of systemic sclerosis. Nat Rev Rheumatol, 14, 657–673.

Bagnato

, Roberts

W.N.

, Roman

, and Gangemi

(2017). A systematic review of overlapping microRNA patterns in systemic sclerosis and idiopathic pulmonary fibrosis. Eur Respir Rev, 26, pii:160125.

Baumjohann

, and Ansel

K.M.

(2013). MicroRNA-mediated regulation of T helper cell differentiation and plasticity. Nat Rev Immunol, 13, 666–678.

Berschneider

, Ellwanger

D.C.

, Baarsma

H.A.

, Thiel

, Shimbori

, White

E.S.

, et al. (2014). miR-92a regulates TGF-beta1-induced WISP1 expression in pulmonary fibrosis. Int J Biochem Cell Biol, 53, 432–441.

Bhattacharyya

, and Varga

(2015). Emerging roles of innate immune signaling and toll-like receptors in fibrosis and systemic sclerosis. Curr Rheumatol Rep, 17, 474.

Chin

C.H.

, Chen

S.H.

, Wu

H.H.

, Ho

C.W.

, Ko

M.T.

, and Lin

C.Y.

(2014). cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol, 8 Suppl 4, S11.

Christmann

R.B.

, Sampaio-Barros

, Stifano

, Borges

C.L.

, de Carvalho

C.R.

, Kairalla

, et al. (2014). Association of interferon- and transforming growth factor beta-regulated genes and macrophage activation with systemic sclerosis-related progressive lung fibrosis. Arthritis Rheumatol, 66, 714–725.

Christmann

R.B.

, Wooten

, Sampaio-Barros

, Borges

C.L.

, Carvalho

C.R.

, Kairalla

R.A.

, et al. (2016). miR-155 in the progression of lung fibrosis in systemic sclerosis. Arthritis Res Ther, 18, 155.

D'Angelo

W.A.

, Fries

J.F.

, Masi

A.T.

, and Shulman

L.E.

(1969). Pathologic observations in systemic sclerosis (scleroderma). A study of fifty-eight autopsy cases and fifty-eight matched controls. Am J Med, 46, 428–440.

10.

Davidson-Moncada

, Papavasiliou

F.N.

, and Tam

(2010). MicroRNAs of the immune system Roles in inflammation and cancer. Ann N Y Acad Sci, 1183, 183–194.

11.

Elhai

, Meune

, Avouac

, Kahan

, and Allanore

(2012). Trends in mortality in patients with systemic sclerosis over 40 years: a systematic review and meta-analysis of cohort studies. Rheumatology (Oxford), 51, 1017–1026.

12.

Hoffmann

, Wilhelm

, Marsh

L.M.

, Ghanim

, Klepetko

, Kovacs

, et al. (2014). Distinct differences in gene expression patterns in pulmonary arteries of patients with chronic obstructive pulmonary disease and idiopathic pulmonary fibrosis with pulmonary hypertension. Am J Respir Crit Care Med, 190, 98–111.

13.

Hoffmann-Vold

A.M.

, Weigt

S.S.

, Palchevskiy

, Volkmann

, Saggar

, Li

, et al. (2018). Augmented concentrations of CX3CL1 are associated with interstitial lung disease in systemic sclerosis. PLoS One, 13, e0206545.

14.

Huang Da

, Sherman

B.T.

, and Lempicki

R.A.

(2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc, 4, 44–57.

15.

Jafarinejad-Farsangi

, Farazmand

, Gharibdoost

, Karimizadeh

, Noorbakhsh

, Faridani

, et al. (2016). Inhibition of microRNA-21 induces apoptosis in dermal fibroblasts of patients with systemic sclerosis. Int J Dermatol, 55, 1259–1267.

16.

Jafarinejad-Farsangi

, Farazmand

, Mahmoudi

, Gharibdoost

, Karimizadeh

, Noorbakhsh

, et al. (2015). MicroRNA-29a induces apoptosis via increasing the Bax:Bcl-2 ratio in dermal fibroblasts of patients with systemic sclerosis. Autoimmunity, 48, 369–378.

17.

John

, Enright

A.J.

, Aravin

, Tuschl

, Sander

, and Marks

D.S.

(2004). Human microRNA targets. PLoS Biol, 2, e363.

18.

Khalil

, Xia

, Bodempudi

, Kahm

, Hergert

, Smith

, et al. (2015). Pathologic regulation of collagen I by an aberrant protein phosphatase 2A/histone deacetylase C4/microRNA-29 signal axis in idiopathic pulmonary fibrosis fibroblasts. Am J Respir Cell Mol Biol, 53, 391–399.

19.

Kim

H.H.

, Yoon

J.W.

, P, Y.E., Cho

E.G.

, Sohn

J.K.

, Kim

T.K.

, et al. (2006). Cryopreservation of potato cultivated varieties and wild species: critical factors in droplet vitrification. Cryo Lett, 27, 223–234.

20.

Kohl

, Wiese

, and Warscheid

(2011). Cytoscape: software for visualization and analysis of biological networks. Methods Mol Biol, 696, 291–303.

21.

Kouri

F.M.

, Hurley

L.A.

, Daniel

W.L.

, Day

E.S.

, Hua

, Hao

, et al. (2015). miR-182 integrates apoptosis, growth, and differentiation programs in glioblastoma. Genes Dev, 29, 732–745.

22.

Lappi-Blanco

, Lehtonen

S.T.

, Sormunen

, Merikallio

H.M.

, Soini

, and Kaarteenaho

R.L.

(2013). Divergence of tight and adherens junction factors in alveolar epithelium in pulmonary fibrosis. Hum Pathol, 44, 895–907.

23.

Lepparanta

, Sens

, Salmenkivi

, Kinnula

V.L.

, Keski-Oja

, Myllarniemi

, et al. (2012). Regulation of TGF-beta storage and activation in the human idiopathic pulmonary fibrosis lung. Cell Tissue Res, 348, 491–503.

24.

Lewis

B.P.

, Burge

C.B.

, and Bartel

D.P.

(2005). Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell, 120, 15–20.

25.

Makino

, and Jinnin

(2016). Genetic and epigenetic abnormalities in systemic sclerosis. J Dermatol, 43, 10–18.

26.

Marasini

, Cossutta

, Selmi

, Pozzi

M.R.

, Gardinali

, Massarotti

, et al. (2005). Polymorphism of the fractalkine receptor CX3CR1 and systemic sclerosis-associated pulmonary arterial hypertension. Clin Dev Immunol, 12, 275–279.

27.

Milano

, Pendergrass

S.A.

, Sargent

J.L.

, George

L.K.

, McCalmont

T.H.

, Connolly

M.K.

, et al. (2008). Molecular subsets in the gene expression signatures of scleroderma skin. PLoS One, 3, e2696.

28.

Oliver

J.R.

, Kushwah

, and Hu

(2012). Multiple roles of the epithelium-specific ETS transcription factor, ESE-1, in development and disease. Lab Invest, 92, 320–330.

29.

Pandit

K.V.

, Corcoran

, Yousef

, Yarlagadda

, Tzouvelekis

, Gibson

K.F.

, et al. (2010). Inhibition and role of let-7d in idiopathic pulmonary fibrosis. Am J Respir Crit Care, 182, 220–229.

30.

Rajasekaran

, Vaz

, and Reddy

S.P.

(2012). Fra-1/AP-1 transcription factor negatively regulates pulmonary fibrosis in vivo. Plos One, 7, e41611.

31.

Reimers

, and Carey

V.J.

(2006). Bioconductor: an open source framework for bioinformatics and computational biology. Methods Enzymol, 411, 119–134.

32.

Roth

M.D.

, Tseng

C.H.

, Clements

P.J.

, Furst

D.E.

, Tashkin

D.P.

, Goldin

J.G.

, et al. (2011). Predicting treatment outcomes and responder subsets in scleroderma-related interstitial lung disease. Arthritis Rheum, 63, 2797–2808.

33.

Sakamoto

, Oda

, Iwamoto

, and Tsuneyoshi

(1999). A comparative study of fibrous dysplasia and osteofibrous dysplasia with regard to expressions of c-fos and c-jun products and bone matrix proteins: a clinicopathologic review and immunohistochemical study of c-fos, c-jun, type I collagen, osteonectin, osteopontin, and osteocalcin. Hum Pathol, 30, 1418–1426.

34.

Servettaz

, Agard

, Tamby

M.C.

, Guilpain

, Guillevin

, and Mouthon

(2006). [Systemic sclerosis: pathophysiology of a multifaceted disease]. Presse Med, 35, 1903–1915 (Article in French).

35.

Shannon

, Markiel

, Ozier

, Baliga

N.S.

, Wang

J.T.

, Ramage

, et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res, 13, 2498–2504.

36.

Smyth

G.K.

(2004). Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol, 3, Article3.

37.

Szklarczyk

, Franceschini

, Kuhn

, Simonovic

, Roth

, Minguez

, et al. (2011). The STRING database in 2011: functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res, 39, D561–D568.

38.

Takada

, Gresh

, Bozec

, Ikeda

, Kamiya

, Watanabe

, et al. (2011). Interstitial lung disease induced by gefitinib and Toll-like receptor ligands is mediated by Fra-1. Oncogene, 30, 3821–3832.

39.

Tanaka

, Suto

, Ikeda

, Sanayama

, Nakagomi

, Iwamoto

, et al. (2013). Alteration of circulating miRNAs in SSc: miR-30b regulates the expression of PDGF receptor beta. Rheumatology, 52, 1963–1972.

40.

Tyndall

A.J.

, Bannert

, Vonk

, Airo

, Cozzi

, Carreira

P.E.

, et al. (2010). Causes and risk factors for death in systemic sclerosis: a study from the EULAR Scleroderma Trials and Research (EUSTAR) database. Ann Rheum Dis, 69, 1809–1815.

41.

Vlachos

I.S.

, Kostoulas

, Vergoulis

, Georgakilas

, Reczko

, Maragkakis

, et al. (2012). DIANA miRPath v.2.0: investigating the combinatorial effect of microRNAs in pathways. Nucleic Acids Res, 40, W498–W504.

42.

Wang

, Huang

, Zhang

, Chen

, and Zhao

(2018). Molecular pathogenesis involved in human idiopathic pulmonary fibrosis based on an integrated microRNA mRNA interaction network. Mol Med Rep, 18, 4365–4373.

43.

Wong

, and Wang

(2015). miRDB: an online resource for microRNA target prediction and functional annotations. Nucleic Acids Res, 43, D146–D152.

44.

, Duan

, Cao

, Field

, Newnham

C.M.

, Koehler

D.R.

, et al. (2008). Regulation of epithelium-specific Ets-like factors ESE-1 and ESE-3 in airway epithelial cells: potential roles in airway inflammation. Cell Res, 18, 649–663.

45.

Young

N.A.

, Valiente

G.R.

, Hampton

J.M.

, Wu

L.C.

, Burd

C.J.

, Willis