Correlation Network Analysis Reveals Relationships between MicroRNAs, Transcription Factor T-bet, and Deregulated Cytokine/Chemokine-Receptor Network in Pulmonary Sarcoidosis

Sarcoidosis is an inflammatory granulomatous disease with unknown etiology driven by cytokines and chemokines. There is limited information regarding the regulation of cytokine/chemokine-receptor network in bronchoalveolar lavage (BAL) cells in pulmonary sarcoidosis, suggesting contribution of miRNAs and transcription factors. We therefore investigated gene expression of 25 inflammation-related miRNAs, 27 cytokines/chemokines/receptors, and a Th1-transcription factor T-bet in unseparated BAL cells obtained from 48 sarcoidosis patients and 14 control subjects using quantitative RT-PCR. We then examined both miRNA-mRNA expressions to enrich relevant relationships. This first study on miRNAs in sarcoid BAL cells detected deregulation of miR-146a, miR-150, miR-202, miR-204, and miR-222 expression comparing to controls. Subanalysis revealed higher number of miR-155, let-7c transcripts in progressing (n = 20) comparing to regressing (n = 28) disease as assessed by 2-year follow-up. Correlation network analysis revealed relationships between microRNAs, transcription factor T-bet, and deregulated cytokine/chemokine-receptor network in sarcoid BAL cells. Furthermore, T-bet showed more pronounced regulatory capability to sarcoidosis-associated cytokines/chemokines/receptors than miRNAs, which may function rather as “fine-tuners” of cytokine/chemokine expression. Our correlation network study implies contribution of both microRNAs and Th1-transcription factor T-bet to the regulation of cytokine/chemokine-receptor network in BAL cells in sarcoidosis. Functional studies are needed to confirm biological relevance of the obtained relationships.


Introduction
Pulmonary sarcoidosis is an inflammatory disorder of unknown etiology characterized by the accumulation of activated Th1/Th17 lymphocytes and macrophages in the alveoli and subsequent granuloma formation [1][2][3]. The key role in the pathogenesis of sarcoidosis is played by proinflammatory cytokines and chemokines, molecules crucially involved in the activation of immune and inflammatory cells and their trafficking to the site of disease [4]. However, there is still limited information about the regulation of the cytokine/chemokine-receptor network in pulmonary sarcoidosis and its phenotypes.
There is a growing body of evidence that the regulation of inflammatory response is a very complex process involving coordinated participation of multiple regulation systems, such as an integrated network of microRNAs (miR-NAs) and transcription factors [5,6]. The emerging role of 2 Mediators of Inflammation miRNAs, a class of single-stranded noncoding RNAs of [19][20][21][22][23][24][25] nucleotides in length, in regulation of inflammatory response has been already reported in chronic pulmonary diseases such as asthma [7] and chronic obstructive pulmonary disease [8]. In sarcoidosis, altered miRNA pattern has been reported in lung tissues [9], peripheral blood mononuclear cells [9][10][11], and serum [10]. However, there is no information regarding the miRNA pattern in bronchoalveolar lavage (BAL) cells and their regulatory capability related to cytokine/chemokine-receptor network in pulmonary sarcoidosis. Also, a Th1-transcription factor T-bet has emerged as key regulator of crucial immune genes such as interferon-(IFN-) and chemokine receptor CXCR3 in sarcoid inflammation [12][13][14] as well as in other inflammatory conditions [15][16][17]. However, no information about the possible cooperation of this Th1-transcription factor and inflammation-related microRNAs in regulation of cytokine/chemokine-receptor network in BAL cells in sarcoidosis exists yet.
In this study, we aimed to investigate the gene expression pattern of candidate inflammation-related miRNAs in BAL cells obtained from sarcoidosis patients and control subjects. In order to assess the possible contributions of miRNAs as posttranscriptional regulators and T-bet as a driver Th1transcription factor on sarcoid inflammation, we searched for relationships between miRNAs and T-bet with sarcoidosisassociated cytokine/chemokine-receptor network in BAL cells obtained from patients with sarcoidosis and subgroups with progressing and regressing disease as assessed by 2-year follow-up. We believe that understanding the transcriptional and posttranscriptional regulation of cytokine/chemokinereceptor network could shed light on the cause and progression of pulmonary sarcoidosis and other inflammatory and autoimmune diseases and eventually lay the groundwork for therapeutic options.

Subjects.
Patients were further subdivided according to the disease development as assessed by 2-year followup. BAL was performed according to a standard protocol [18] in 48 patients with pulmonary sarcoidosis (S) and 14 control subjects (C) of Czech origin. The diagnosis of sarcoidosis was determined in compliance with the criteria from the Statement on Sarcoidosis [19]. No patient received corticosteroid treatment before BAL. Patients were further subdivided according to the disease development as assessed by the 2-year follow-up: (i) patients with progressing disease (Prog, = 20) and (ii) those where the regression was achieved (Reg, = 28). The control group consisted of 14 subjects undergoing BAL as a part of clinical investigation for psychogenic cough, cough associated with gastroesophageal reflux disease and lung hypertension. For clinical and laboratory characteristics of enrolled patients and control subjects, see Table 1.
All patients gave their informed consent for the use of BAL, taken primarily for diagnostic evaluation, for the purpose of this study. The local ethical committee of Palacký University and Faculty Hospital, Olomouc, approved the study.

BAL Sample Processing, miRNA/mRNA Isolation, and
Reverse Transcription. BAL cells were separated from the fluid by centrifugation and total RNA was isolated from unseparated BAL cells with mirVana miRNA kit (Ambion, Austin, USA); RNA quality and quantity were measured by 2100 Bioanalyzer using RNA 6000 Nano assays (Agilent Technologies, Palo Alto, USA). The reverse transcription of miRNA was performed with TaqMan microRNA Reverse Transcription kit (Applied Biosystems, Foster City, CA) using stem-loop RT primers, ensuring RT efficiency and specificity [20] according to the manufacturer's instructions. The reactions were incubated for 30 min at 16 ∘ C, 30 min at 42 ∘ C, and 5 min at 85 ∘ C. The resulting cDNA was stored at −20 ∘ C until use. The reverse transcription of mRNA was performed with Transcriptor First Strand cDNA Synthesis Kit with anchored dT primers (Roche Applied Science, Indianapolis, USA) according to the manufacturer's recommendation. The reactions were incubated for 10 min at 65 ∘ C, 60 min at 50 ∘ C, and 5 min at 85 ∘ C. The resulting cDNA was stored at −20 ∘ C until use.

Measurement of miRNA/mRNA Expression by Quantita-
tive RT-PCR. The gene expression for each miRNA and mRNA in BAL cells was investigated by qRT-PCR using specific primers and probes (see Table S1 and Table S2, in Supplementary Material available online at http://dx.doi.org/10.1155/ 2015/121378). qRT-PCR for each individual miRNA was performed in a 20 L reaction mixture that included 1.3 L of diluted RT product, 1 L of 20X TaqMan Individual microRNA assay, 10 L of 2X TaqMan Universal PCR Master Mix, No AmpErase UNG (Applied Biosystems), and 7.7 L of nuclease-free water. qRT-PCR for each individual mRNA expression was performed as described previously [21]. All reactions were performed on RotorGene3000 system (Qiagen Inc., Valencia, CA, USA); the reaction steps were as follows: 95 ∘ C for 10 min, followed by 40 cycles at 95 ∘ C for 15 s and 60 ∘ C for 1 min. The relative miRNA and mRNA expression levels were calculated by a second-derivative method (RotorGene Software 6.1.81, Corbett, Sydney, Australia); cDNA from human universal reference RNA (Stratagene, La Jolla, CA, USA) was used as a calibrator. A reference gene for miRNA analysis was endogenous control Mammalian U6 and for mRNA a reference gene PSMB2 [21]. Changes in expression levels are presented as mean relative expression with 95% confidence interval (CI).

Selection of Candidate miRNAs and Identification of
Binding Sites. The miRNA pathway analysis web server DIANA-mirPath v.3 (http://www.microrna.gr/miRPathv3/) [22] was used to nominate the candidate inflammationrelated miRNAs based on their possible involvement in cytokine-chemokine interaction pathway. For the identification of binding sites between candidate miRNAs and mRNAs, we used the mirSystem [23] (Table S3, Supplementary Material). This web-based tool integrates the seven most often used target gene prediction algorithms: DIANA, miRanda, miRBridge, PicTar, PITA, rna22, and TargetScan. Moreover, it contains validated data of interactions between particular 2.0 ± 0.5 (0.5-3.7) BALF: bronchoalveolar lavage fluid; : number of patients; n: no; y: yes; ex: ex-smoker; CXR: chest X-ray; -: not relevant. Data are presented as mean ± SD (minimum and maximum in parentheses). † Data were not available for four patients. # % of CD3, CD4, CD8, and CD19 refers to total lymphocyte counts. miRNA and its target genes from TarBase and miRecords databases.
2.5. Clustering Using Kohonen Self-Organizing Neural Networks. Kohonen self-organizing neural network (self-organizing maps, SOM), a clustering tool, was applied to find clusters of input data that are very close to each other [24]. The input data (the whole data set including miRNA and mRNA expression of all studied molecules) for each sarcoidosis patient were transformed to vectors, which were recorded in the neural network. Neurons in the cortex were organized in 2D; only the adjacent neurons were interconnected. If the data clustered to neurons in SOM show the internal clustering structure corresponding to the analyzed subgroups (in our case progressing and regressing disease), the data set has high potential to contain relevant biological relationships.
2.6. miRNA/mRNA Correlation Analysis. Correlation matrices were computed to investigate the relationships between miRNA expression and mRNA expression of T-bet and members of sarcoidosis-associated cytokine/chemokine-receptor network represented by Spearman correlation coefficient. The correlation matrices were graphically presented using heat maps, where a hierarchical agglomerative clustering analysis was performed to show the relationships between groups of miRNA and mRNA in the form of clusters. The colour of each cell of the heat map corresponds to the value of Spearman correlation coefficient between given miRNA and corresponding mRNA.

Circos Diagrams for Visualization of miRNA/mRNA Relationships.
To graphically represent the relationships between miRNAs, mRNAs, and T-bet, a Chord diagram generated in Circos system [25] (http://mkweb.bcgsc.ca/tableviewer/) was applied for (i) sarcoidosis patients and subgroups of patients with (ii) progressing sarcoidosis and (iii) regressing sarcoidosis. In the diagram, only significant correlations are depicted ( < 0.05). The intensity of correlations ( ) between individual miRNA-mRNA pairs corresponds to absolute value of Spearman correlation coefficient ( ), for which the mathematical transformation was performed to accentuate the differences in the intensities of correlations. For the transformation, the following formula was used: where represents the threshold for the highest significant value.

miRNA/mRNA/T-bet Network Analysis.
For schematic representation of the relationships between miRNAs, Tbet, and cytokines/chemokines/receptors, we constructed weighted gene coexpression network [26]. Further, an algorithm based on analysis of the nearest neighbors between the studied molecules (represented as vertices) was applied [27]. The larger vertices (spheres) in the network have more nearest neighbors based on correlation analysis than smaller vertices. Only edges connecting the nearest neighbors (pairs with the highest correlations) were preserved. In other words, the size of the vertices (spheres) and connection among vertices show relationships between the investigated molecules in the network.
2.9. Statistical Analysis. Data analysis was performed using GraphPad Prism 5.01 (GraphPad Software, La Jolla, CA, USA) and SPSS 16.0 (SPSS Inc., Chicago, IL, USA). Differences in miRNA and mRNA levels between the cohorts were assessed by nonparametric Kruskal-Wallis one-way analysis-of-variance-by-ranks test; nonparametric Mann-Whitney test was utilized to determine significant differences between two groups. Spearman correlation coefficient and its corresponding values were computed using R statistical software package (http://www.r-project.org/). A value < 0.05 was considered significant.

Analysis of miRNA Expressions in Sarcoidosis Patients
with Progressing/Regressing Disease. In order to investigate the miRNA expression pattern in sarcoidosis patients subdivided according to disease outcome after 2-year followup, the expression levels of a set of candidate miRNAs were determined in BAL cells from progressing and regressing sarcoidosis. To exclude the notion that variations in the levels of expression of transcripts might be related to changes in cell populations, we investigated the distribution of the absolute and relative number of BAL macrophages and lymphocytes and revealed no difference between subgroups of patients with progressing and regressing sarcoidosis ( > 0.05).
When comparing to regressing sarcoidosis, elevated expression of miR-155 ( = 0.017) and let-7c ( = 0.039) (Table 2, Figure 2(b)) was detected in progressing disease. No difference in expression levels of other studied miRNAs was observed between patients with progressing and regressing sarcoidosis (Figure 1(a)).

Self-Organizing Neural Networks for Sarcoidosis Patients.
In order to assess the inner data structure in our data sets, we applied the Kohonen SOM. The unsupervised clustering analysis classified the sarcoidosis patients into the following clusters represented by neurons (hexagons): (1) green Mediators of Inflammation 5    Color key and histogram Regressing sarcoidosis Progressing sarcoidosis neurons including both classes visualized as interpolation between red and green colors (Figure 3). This analysis supports the presence of inner structures in patient subgroups as assessed by disease outcome, thus supporting the fact that the subgroups differ in miRNA/mRNA profiles. Three clusters exactly matched patients with regressing sarcoidosis (25%) and one cluster patients with progressing disease (10.4%). Two clusters more or less correspond to patients with regressing sarcoidosis (27.1%), whereas one cluster more or less corresponds to patients with progressing disease (16.7%). On the other hand, two clusters of patients (20.8%) do not correspond to any class. Interestingly, the patients with regression formed more clusters than progressing patients, thus indicating that this subgroup is more heterogeneous regarding the expression profiles of studied miRNAs and mRNAs comparing to progressing disease.

Correlation Analysis of miRNAs, T-bet, and Cytokine/Chemokine-Receptor Network.
In order to assess the relationships between studied miRNAs and mRNAs in BAL cells, we performed correlation analysis. To distinguish the degree of correlation between miRNA-mRNA pairs, we set up cut-off for Spearman correlation coefficient corresponding to significant value ( < 0.05). Since no significant correlations between the studied molecules were detected    (Figure 4(a)).

Weighted Gene Coexpression Network Analysis for miRNA-mRNA Relationships.
In order to better characterize the relationships between studied miRNAs, mRNAs, and T-bet, we performed weighted gene coexpression network analysis.
Regarding control group, the investigated molecules were not closely related to each other (small vertices) and their relationship pattern differed from sarcoidosis-associated pattern ( Figure 5(d)).

Discussion
This is the first study investigating the suggested contribution of both transcriptional and posttranscriptional control to deregulated cytokine/chemokine-receptor network in pulmonary sarcoidosis. Our correlation network analysis revealed that both inflammatory-related microRNAs CXCR4 C X C R 3 C X C L 9 C X C L 2 C X C L 1 6 C X C L 1 2 C X C L 1 1  CXCR7  CXCL16  CXCR4  CXCR3  CXCR6  IL15RA  IL2RB  CCR1  CCR2B  CXCL3  CCR2A  abs ly  IFNG  CCL5  IL2RA  IL2  CXCL9  CXCL11  CXCL10  CCL7  CCL2  CCR5  CCL4  CCL3  CCL8  CXCL2  CCL13  CCL19   and a key Th1-transcription factor T-bet may contribute to deregulated cytokine/chemokine-receptor network in sarcoid BAL cells, whereas contributions of transcriptional and posttranscriptional regulation differ between progressing and regressing disease as assessed by 2-year follow-up. We also for the first time investigated the gene profile of inflammationrelated miRNAs in BAL cells and revealed altered miRNA pattern in BAL cells obtained from patients with sarcoidosis comparing to control subjects and also from patients with progressing versus regressing disease.
In sarcoidosis pathogenesis, the crucial role is attributed to inflammatory cytokines and chemokines, which direct circulating leukocytes to the sites of inflammation and control leukocyte activation and cytokine production, angiogenesis, and Th cell polarization [4,28,29]. In view of the crucial role of cytokines, chemokines, and their receptors in the pathogenesis of sarcoidosis, a fundamental question arises: how the cytokine/chemokine-receptor network is regulated. There is first evidence that cytokines and chemokines are regulated by both transcriptional and posttranscriptional mechanisms  allowing for rapid and transient production in response to a variety of stimuli [30,31]. In line with this observation, we hypothesized that the regulation of cytokine/chemokinereceptor network in BAL cells also involves coordinated participation of multiple regulatory systems, including miRNAs and a key sarcoidosis-associated Th1-transcription factor Tbet.
However, miRNAs pattern in BAL cells has not been investigated yet. We therefore investigated for the first time expression pattern of inflammation-related miRNAs in BAL cells and revealed downregulated miR-202 and miR-204 and upregulated miR-146a, miR-150, and miR-222 expression in sarcoidosis comparing to controls. When comparing progressing to regressing sarcoidosis, a higher number of miR-155 and let-7c transcripts were detected in progressing sarcoidosis. All detected miRNAs associated with sarcoidosis were already associated with regulation of inflammatory response. Inflammatory miR-155 and miR-146a were found to be upregulated by NF-B [32,33], miR-150 involved in pulmonary inflammation [34], miR-202 in rheumatoid arthritis [35], and miR-204 in the production of inflammatory mediators and apoptosis [36]. MiR-222 contributed to inflammation-mediated neovessel formation [37], and let-7c regulated allergic airway inflammation [38]. Moreover, miR-204 was also found to be deregulated in sarcoid lung tissues reported by Crouser et al. [9]. Interestingly, when comparing the miRNA profile in sarcoid BAL cells obtained in this study with profile reported in sarcoid lung tissues [9], we obtained different pattern of top-deregulated miRNAs (miR-92a, miR-125a, miR-126, miR-181a, miR-206, and miR-302) found in sarcoid lung tissues, thus supporting the fact that miRNA expression signatures vary depending on the tissues and cell types examined.
However, it is still unclear how miRNAs regulate cytokines, because 3 -untranslated region (3 -UTR) of most cytokines lacks critical "seed region" for binding to mRNA [39]. Current studies showed that cytokine and other immune genes may be regulated by other elements (e.g., ARE machinery components) [40] or through highly biologically functional "nonseed" miRNA target sites, containing single mismatches, GU wobbles, insertions, or deletions in the seed-match regions [41][42][43]. These "nonseed" miRNA target sites represent up to 50% of all miRNA-mRNA interactions and remain uncovered, despite their functional significance, using current target prediction algorithms [39]. Importantly, miRNA regulation may not be directly at the cytokine/chemokine or receptor level but rather at upstream or downstream steps in the pathway. These facts may explain our observation of correlation between numerous miRNAs and cytokine expression (e.g., CXCR4 and miR-24), where no binding site was discovered by prediction algorithms as well as numerous positive correlations between miRNAs and cytokine and chemokine expression (e.g., CCL2 and miR-214) as reported also in other studies [44,45]. Moreover, the unique role of miRNAs in inflammation is that the miRNAs modulate the expression of target genes to an optimum level rather than participating in on/off decisions [40].
Besides posttranscriptional regulation by miRNAs, also transcriptional control by transcription factors may contribute to deregulated cytokine/chemokine-receptor network in pulmonary sarcoidosis. We recently reported crucial role for a Th1-transcription factor T-bet, which has emerged as key regulator of IFNG and chemokine receptor CXCR3 and other immune genes in sarcoid inflammation [12]. In the present study, we searched for evidence of possible participation of transcriptional and posttranscriptional regulation by miRNAs and T-bet in controlling the inflammation in sarcoid BAL cells. Our correlation analysis revealed relationship of T-bet and sarcoid inflammation, namely, its association with key sarcoidosisassociated cytokine IFNG, and receptors IL2RB, IL15RA, CXCR3, and CXCR6, regardless of the disease outcome.
Importantly, the relationships between T-bet and particular miRNAs regarding the chemokine network changed between subgroups of patients, who differed in disease outcome after 2-year follow-up. In progressing sarcoidosis, T-bet correlated with CCL5, CXCR4, and CXCR7. We also observed positive correlations between numerous miRNAs and expression of CCL2, CCL7, and CXCL12 and negative correlation with CCL5, CXCR3, CXCR4, and CXCR7 expression in progressing disease. In regressing disease, T-bet showed less pronounced correlation with IL2 and receptors CXCR4, CXCR7, and IL2RA than in progression. Also miRNA pattern changed; in regression, the correlation analysis revealed positive correlation between miR-146a and CCL19 as well as miR-148a and CXCR4. Our analysis therefore nominated CCL2 and CCL5, both chemoattractants of mononuclear and mast cells [46], as well as CCL19, a chemokine implicated in T-lymphocyte recruitment [47], as key candidate chemokines associated with prognosis of pulmonary sarcoidosis. Similarly, CXCR4 and CXCR7 and their ligand, CXCL12, playing a role in regulation of leukocyte mobilization and trafficking [48,49], deserve further investigation.
To better characterize the network of inflammationrelated miRNAs, T-bet, and cytokine/chemokine-receptor network in sarcoid BAL cells, we performed weighted correlation network analysis. This analysis further supported the crucial relationships between T-bet, CCL5, IFNG, IL2RB, and IL15RA, together with let-7d and miR-202 in progressing sarcoidosis. In regressing disease, T-bet was more closely related to CXCR3 and less to IFNG and CCL5 than in progressing disease. Regarding miRNAs in regressing disease, miR-212, miR-146a, and let-7d were highlighted. Given the complexity of miRNA-mRNA-transcription factor network, the understanding of the contribution of network members and regulatory mechanisms is intricate and requires further studies which may reveal the functional outcome of the interactions between the associated molecules.
Our study has several limitations. First, this exploratory study on miRNA and mRNA expression pattern and their relationships was performed in unseparated BAL cells. To obtain complete picture of the regulatory processes ongoing in sarcoid inflammation, this study should continue by analysing the miRNA/mRNA expressions in distinct cell subpopulations. We are also aware that our set of inflammationrelated miRNAs and chemokine/receptor genes, selected based on their involvement in the pathogenesis of sarcoidosis, does not cover the whole potential network of molecules/pathways contributing to this disease. Further, we did not perform the functional studies needed to confirm the biological relevance of the obtained relationships between Tbet, miRNAs, and chemokine network. However, we believe that our study may contribute to nomination of interesting relationships for future studies of regulatory mechanisms involved in the pathogenesis of pulmonary sarcoidosis and its phenotypes.

Conclusions
Our correlation network analysis implies both microR-NAs and Th1-transcription factor T-bet in the regulation of cytokine/chemokine-receptor network in BAL cells in sarcoidosis. Future functional studies are however needed to confirm the biological relevance of the obtained relationships. The correlation analysis showed more pronounced regulatory capability of T-bet to sarcoidosis-associated chemokine receptors and cytokines than miRNAs, which function rather as "fine-tuners" of cytokine/chemokine gene expression. Moreover, we reported altered miRNA pattern in BAL cells obtained from sarcoidosis patients as well as from subgroups of patients with progressing and regressing sarcoidosis as assessed by 2-year follow-up. The knowledge and understanding of the regulatory mechanisms underlying the abnormal inflammatory response in sarcoid lungs could shed light on the cause and progression of sarcoidosis and many inflammatory and autoimmune diseases.