Identification of the Core MicroRNAs and Potential Molecular Mechanismsin Sarcoidosis Using Bioinformatics Analysis

Sarcoidosis is a systemic heterogeneous inflammatory disease; however, the etiology and pathogenesis of sarcoidosis are still unknown. Herein, we investigated the core microRNAs and potential molecular mechanisms in sarcoidosis. The DE-miRNAs were diagnosed using the LIMMA software package. DIANA-mirPath was employed to perform pathway and GO enrichment analysis of the DE-miRNAs. PPI networks and miRNA-target gene regulatory networks were used to obtain insight into the actions of DE-miRNAs. Expression of the hub genes along with miRNAs was validated in clinical specimens. Overall, 266 DE-miRNAs were screened. Among these DE-miRNAs, hsa-miR-144, hsa-miR-126, as well as hsa-miR-106a were the upmost upregulated miRNAs; hsa-miR-151-3p, hsa-miR-320d, and hsa-miR-324-3p were the top downregulated miRNAs. NR3C1, ZBTB7A, NUFIP2, BZW1, ERGIC2, and VEGFA were mapped as the most targeted hub genes in the upregulation of miRNAs, and MCL1 and SAE1 were the most targeted hub genes in the downregulation of miRNA. VEGFA and NR3C1 were selected and potentially modulated by hsa-miR-20b, hsa-miR-126, and hsa-miR-106a. In sarcoidosis pathological tissue, hsa-miR-126 was highly expressed, and VEGFA and NR3C1 were overexpressed. In conclusion, our results revealed the dysregulation of hsa-miR-126 and a potential regulatory mechanism for pathogenesis in sarcoidosis.


INTRODUCTION
Sarcoidosis is a systemic heterogeneous inflammatory disease that usually leads to a multisystemic granulomatous disorder (Iannuzzi et al., 2007;Jain et al., 2020). It can affect any organ; however, at least 90% of cases occur in intrathoracic lymph nodes, as well as the lungs (Ungprasert et al., 2019). The annual incidence of sarcoidosis changes with race and ethnicity and is highest in African Americans (17 to 35/100,000) and lowest in Asians (1/100,000). However, the incidence may be underestimated, as almost 50% of the sarcoidosis patients are asymptomatic, especially those in the early stage of the disease. Epidemiologic studies have reported that the maximum incidence of sarcoidosis occurs in female non-smokers with an average age of 30 (Jain et al., 2020).
The characteristic pathologic lesion of pulmonary sarcoidosis is a non-necrotizing granuloma (Iannuzzi et al., 2007). The etiology and pathogenesis of sarcoidosis are still unknown, and genetic susceptibility, autoimmunity, and environmental factors are all considered to be involved in the development of sarcoidosis.
MicroRNAs (miRNAs) are small, single-stranded, noncoding RNAs approximately 22 nucleotides long, some of which have been discovered to perform significant modulatory functions in animals through targeting the messages of protein-coding genes for translational inhibition. Research has documented that changes in the expression levels of miRNAs participate in multiple pathologies, as well as diseases. A previous study evaluated the expression trend of multiple miRNAs in sarcoidosis, and selected miRNAs were identified. However, regulatory miRNAs, their target genes, and the correlated protein expression have seldom been explored (Pattnaik et al., 2020). In this regard, we examined the microarray data of the gene expression profile of GSE26409 by a series of biological informatics approaches. The potential target genes of the upmost up-regulated, as well as down-regulated miRNAs were prognosticated by miRTarBase, and their potential functions were evaluated by DIANA-mirPath. Then, a protein-protein interaction (PPI) and miRNA-gene modulatory network were established by Cytoscape. Moreover, we validated the expression status of these hub genes and their modulated microRNAs in sarcoidosis samples by immunohistochemistry (IHC) along with fluorescence in situ hybridization (FISH). Based on a gene chip array, we used biological information technology to reveal the potential pathway in the etiology of sarcoidosis and yield additional information for its diagnosis, prognosis and treatment.

Microarray Data
The microarray GSE26409 cohort was selected for further study and abstracted from the National Center for Biotechnology Information GEO (Gene Expression Omnibus) website 1 , which is available online. This data set was based on the GPL9040 platform (febit Homo Sapiens miRBase 13.0) and included 45 samples from sarcoidosis patients and 55 from normal controls.

Screening for DE-miRNAs
First, normalization of the downloaded data was performed using the Normalize Between Array function in the Bioconductor R package "LIMMA" (available online: http://www.bioconductor. org/). Then, the unpaired Student's t-test was employed in comparing the sarcoidosis group with the control group. The cut-off criterion for screening differentially expressed miRNAs (DE-miRNAs) were p < 0.05 along with | fold change (FC)| > 1. Additionally, the GEO2R online analytic tool based on GEO data platform was employed further to validate the DE-miRNAs.  (Vlachos et al., 2015). The functional, as well as pathway enrichment of the candidate miRNAs were explored and annotated by mirPath v2.0 based on the database Tarbase v7. The GO (Gene Ontology) was annotated through the DIANA-mirPath online tool on the selected DE-miRNAs. The KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis of DE-miRNAs was also performed by using DIANA-mirPath. p < 0.05 signified statistical significance.

PPI Network Integration
To analyze the connection among proteins, the identify target genes of the upmost 5 most up and downregulated DE-miRNAs were uploaded to the STRING data resource 2 , and the results were visualized in Cytoscape 3.7.1 (Damian et al., 2015). Furthermore, we screened the significant nodes or hub genes based on degree, and then the miRNA-hub gene networks were mapped with Cytoscape 3.7.1.

Immunohistochemistry (IHC)
Immunohistochemistry was used to further validate the paraffin sections of mediastinal lymph node tissues from pulmonary sarcoidosis patients and tuberculosis patients (control group), which came from the Second Affiliated Hospital of Xi'an Jiaotong University, with ethical approval. All specimens were obtained by endobronchial ultrasound-guided transbronchial needle aspiration (EBUS-TBNA), and all the patients in the study group and the control group have been screened, and those with cardiovascular and cerebrovascular diseases, endocrine, infectious diseases, tumor diseases, etc., have been excluded. The diagnosis was made based on guidelines (Costabel and Hunninghake, 1999). The tissue paraffin sections were first deparaffinized with dimethylbenzene and ethyl alcohol for 15 min, repaired by boiling for 8 min in the repairing solution with EDTA (pH 9.0), digested for 10 min by 3% H 2 O 2 at 37 • C, washed with PBS solution and blocked for 30 min with 3% BSA at room temperature (RT). Afterward, the sections were overnightincubated with the primary antibody (NR3C1, VEGFA) at 4 • C (Abcam, Cambridge, United States). The next day, the tissues were inoculated with the horseradish-peroxidase (HRP) and incubated at RT for 60 min and incubated with diaminobenzidine (DAB) color liquid after washing. Afterward, the tissues were counterstained using Harris hematoxylin, 1% hydrochloric acid alcohol, and lithium carbonate blue sequentially, and then visualized using an optical microscope (10×, 40×, Nikon, Japan) after routine cleaning and drying.

Identification of DE-miRNAs and Prediction of Target Genes
The microarray GSE26409 cohort was acquired from the GEO website, and this cohort consisted of 45 cases of sarcoidosis (GSM648250-GSM648294) and 55 normal controls (GSM648195-GSM648249). After normalization (Figure 1), processing of the data was conducted via unpaired t-test, with the cutoff criterion of P < 0.05 along with | log2FC| > 1). Overall, 266 DE-miRNAs were finally screened (Figures 2, 3), which constituted 132 and 134 upregulated and downregulated miRNAs, respectively. The top ten most upregulated, as well as downregulated miRNAs are indicated in Table 1. On the basis of the fold change (FC), hsa-miR-144, hsa-miR-126, and hsa-miR-106a were the upmost 3 upregulated miRNAs and hsa-miR-151-3p, hsa-miR-320d, and hsa-miR-324-3p were the upmost 3 downregulated miRNAs. For the three upregulated miRNAs, 580 potential target genes were investigated, and 543 genes were predicted for the 3 downregulated miRNAs via the miRTarBase platform.

Functional Enrichment Analysis
To determine the enriched pathways and process enrichment of these target genes, we afterward carried out a KEGG pathway enrichment and GO functional annotation. The KEGG pathway analysis included non-small cell lung cancer, proteoglycans in cancer, viral carcinogenesis, etc. (Figure 4A, Table 2). We selected 3 GO categories consisting of biological pathway (BP), cellular component (CC), along with molecular function (MF) for functional annotation. The top ten GO terms of the target genes of the 10 upmost up-regulated DE-miRNAs are indicated in Figures 4B-D and Table 3, including the responses to stress, catabolic process, gene expression, etc. in the BP category; cellular component, nucleoplasm, cytosol, etc. in the CC category; and enzyme binding, nucleic acid binding transcription factor activity, RNA binding, etc. in the MF category. Subsequently, the target genes of the ten

Development and Analysis of the PPI Network and miRNA-Target Network
The STRING online data resource was used to develop the PPI network complex of the filtered target genes of the 5 upmost upregulated and downregulated miRNAs. By using Cystoscope software, the hub nodes were screened out and further developed the miRNA-hub gene network (Supplementary Figures 1, 2). After simplification, a concise relationship between miRNAs and hub genes was achieved. As shown in Figure 6A, NR3C1, ZBTB7A, NUFIP2, BZW1, ERGIC2, VEGFA, BTG2, and BCL2L11 were the eight most targeted hub genes in the upregulation of miRNAs. Among these, four hub genes could potentially be modulated by the upregulation of hsa-miR-144. Six hub genes and eight hub genes could potentially be modulated by hsa-miR-126 and hsa-miR-106a, respectively. Eight and five hub genes could potentially be modulated by hsa-miR-20b and hsa-miR-363, respectively. Moreover, MCL1 and SAE1 were the two most targeted hub genes in the downregulation of miRNA. MCL1 could potentially be targeted by hsa-miR-30d, hsa-miR-320d, and hsa-miR-151-3p; SAE1 could potentially be regulated by hsa-miR-324-3p, hsa-miR-30d, as well as hsa-miR-151-3p ( Figure 6B).
To further identify reliable hub genes, the targeted hub genes of miRNAs were screened in the gene nodes of the PPI network, and only two overlapping genes (NR3C1, VEGFA) were found with a relatively high degree (19 and 18, respectively) in the up and downregulated hub genes. After combining the above results, NR3C1 and VEGFA were selected and potentially modulated by hsa-miR-126, hsa-miR-20b, and hsa-miR-106a. Of note, these three miRNAs might be potential regulators in manipulating the pathogenesis and development of sarcoidosis.

Verification of the Expression of Hub Genes and miRNAs in Clinical Specimens
To validate the expression of VEGFA and NR3C1, we performed IHC to examine 12 mediastinal lymph node tissues from pulmonary sarcoidosis patients and 13 mediastinal lymph node tissues from tuberculosis patients. Ten of the twelve sections showed moderate to strong dark brown staining of VEGFA in both the epithelioid cells and lymphocytes (Figures 7C,D), and negative to very weak staining for control group (Figures 7A,B).
Nine of the twelve sections presented weak to moderate staining expression of NR3C1 (Figures 7G,H) comparing the control sections which present negative to very weak staining (Figures 7E,F).
The expression contents of hsa-miR-126, hsa-miR-20b, and hsa-miR-106a were further analyzed via FISH on the 12 mediastinal lymph node tissues of pulmonary sarcoidosis patients versus the 13 control tissues. Typical results were shown in Figure 8. The staining revealed strong enrichment of hsa-miR-126 along the edge of the non-caseating granulomas (Figure 8B), presenting higher expression than the controls (Figure 8A), while hsa-miR-20b (Figures 8D,E) and hsa-miR-106a had low expression (Figures 8G,H) in both groups. Quantification of fluorescence for hsa-miR-126, hsa-miR-20b, and hsa-miR-106a were measured, and it demonstrated that hsa-miR-126 significantly increased in mediastinal lymph node tissues of pulmonary sarcoidosis patients (Figures 8C,F,I).
At the molecular level, prediction of the miRNA target genes was performed, then GO annotation and associated cascades were analyzed. Furthermore, by constructing the protein-protein interaction (PPI) networks, NR3C1, ZBTB7A, NUFIP2, BZW1, ERGIC2, VEGFA, BTG2, and BCL2L11 were mapped as the most targeted hub genes in the upregulation of miRNAs, and MCL1 and SAE1 were the most targeted hub genes in the downregulation of miRNA. NR3C1 and VEGFA were selected and are potentially modulated by hsa-miR-126, hsa-miR-20b, and hsa-miR-106a. To validate the bioinformatics analysis data, immunohistochemistry along with fluorescence in situ hybridization were employed to assess gene expression in clinical specimens. Among the identified miRNAs involved in sarcoidosis, hsa-miR-126 was highly expressed, and the two target genes VEGFA and NR3C1 were also highly expressed in lymph node samples of sarcoidosis. The GO annotation and KEGG pathway analysis have identified more cancer related pathways that were closely related to upregulated DE-miRNAs in sarcoidosis, which including proteoglycans in cancer, viral carcinogenesis, non-small cell lung cancer, pancreatic cancer, colorectal cancer, glioma, etc. The relationship between sarcoidosis and cancer is complex, and meta-analyses have concluded that an increased risk of neoplasia in sarcoidosis (El Jammal et al., 2020). Fatty acid biosynthesis and metabolism also participated in the biological functions of miRNA in sarcoidosis. The mismatch of myocardial fatty acid metabolism and myocardial perfusion under positron emission tomography computed tomography (PET-CT) could diagnosis and detection of myocardial injury in cardiac sarcoidosis (Momose et al., 2015).
VEGFA (Vascular endothelial growth factor-A) is a PDGF/VEGF growth factor family member, and by triggering the proliferation, as well as the migration of vascular endothelial cells, it is vital for physiological and pathological of angiogenesis (Tamura et al., 2019). A previous clinical study confirmed that the serum levels of VEGF were higher in sarcoidosis subjects than in normal controls (Mirsaeidi et al., 2018). Further research evaluated the VEGFA level in serum, as well as bronchoalveolar lavage (BAL) and found similar results: VEGFA was remarkably higher in individuals with sarcoidosis relative to the healthy controls (p = 0.0002). The results also showed that the immunoexpression of VEGFA was higher in serum in contrast to the BAL fluid. However, the levels of VEGFA had no correlation with the clinical phenotype (acute onset or insidious onset), stage, or lung function (Piotrowski et al., 2015). Immunohistochemistry analyses on VATS (video-assisted thoracoscopic surgery) lung biopsy samples of sarcoidosis lesions were very consistent with our findings that VEGF was upregulated compared to the control group (Tzouvelekis et al., 2012). NR3C1 (Nuclear receptor subfamily 3 group C member-1) codes for the glucocorticoid receptor and extensively participates in cellular proliferation, differentiation of target tissues, as well as inflammatory responses. Mutations that occurred in NR3C1 are related to generalized glucocorticoid resistance. Glucocorticoids (GCs) are the most common therapeutic agents in clinical use, and due to the fundamental regulatory roles in suppression of inflammation, GCs were widely used for autoimmune and inflammatory diseases (Gharesouran et al., 2018). GWAS (Genome-wide association studies) have found that SNPs (single nucleotide polymorphisms) located in NR3C1 were positively replicated in the corticosteroid drug pathway of asthma (Keskin et al., 2019). After evaluating the SNPs of NR3C1 in in vivo and in vitro models, the results showed that it was closely associated with hypersensitivity to glucocorticoids in asthma (Cuzzoni et al., 2012). In this study, NR3C1 was upregulated in sarcoidosis and expressed in pathological specimens, which may be related to the sensitivity of sarcoidosis to glucocorticoid therapy. Simultaneously, we found that NR3C1 is not always highly expressed, which may be related to the insensitivity of some sarcoidosis patients to glucocorticoid treatment in clinical practice. It can finally be seen that the DNA polymorphism reflects the heterogeneity of the disease.
At present, the role of miRNA in sarcoidosis remains unclear. Kiszałkiewicz et al. investigated the expression of ten screened miRNAs in 94 pulmonary sarcoidosis patients and 50 controls, and remarkable differences were discovered between the sarcoidosis patients and healthy controls in the expression of diverse detected miRNAs in peripheral blood lymphocytes (miR-27b, miR-let-7f, miR-15b, miR-192, miR-130a miR-221, miR-222). Three miRNAs were merely reported to differ in BALF cells (miR-15b, miR-192, as well as miR-221) (Kiszałkiewicz et al., 2016). A subgroup analysis further observed several correlations among the levels of expression of miRNA, lung function parameters (FVC, FEV 1 /FVC), and selected laboratory molecular markers (CD4+/CD8+, serum Ca 2+ concentration). The results showed only gender was not related to the expression of miRNA in BALF (Kiszałkiewicz et al., 2016). Nevertheless, in our study, we did not evaluate the relationship of some clinical indicators to miRNAs expression. We summarized the following reasons. Firstly, the expression of miRNA in peripheral blood is affected by various factors throughout the body, so is the alveolar lavage fluid. For patients with lung infection, COPD, asthma, and tumors, the detected miRNA itself is not accurate enough, subsequently, error in correlation analysis is even greater. Secondly, for early stage sarcoidosis, the lung function is hardly affected, especially FVC, FEV1/FVC, and diffusion function parameters. And, changes in lung function cannot reflect the true situation of patients with sarcoidosis. Thirdly, our study is a retrospective study, not all included patients had lung function tests and samples of alveolar lavage fluid, thus, some of these parameters could not be tested.
To elucidate the function of miRNAs in sarcoidosis, Jazwa et al. (2015) compared the expression levels of selected inflammatory miRNAs in patients with sarcoidosis and healthy individuals. After screening the miRNA transcriptome, miR-34a was isolated from sarcoidosis patients in peripheral blood mononuclear cells, which suggested that it may be involved in the pathology of sarcoidosis. Additionally, miR-155 was shown to be one of the miRNAs that are differentially expressed in macrophages, which are considered to have a role in formation of granuloma in sarcoidosis (Zhang et al., 2013). A study on the mechanism of miRNA in pulmonary sarcoidosis proved that miR-34a along with miR-155 enhance the secretion of interferon-gamma and negatively modulate SIRT-1, therefore inhibiting NF-κB transcription, as well as deacetylate the p53 protein, ultimately inactivating  Barna et al. (2016) found that elevated miRNA-33 promoted the formation of granuloma through the inhibition of the alveolar macrophage lipid transporters ABCA1 and ABCG1. In our study, the lymphatic tissue of tuberculosis patients was used as a control group. This is why we raise the  question from clinic and return to clinic through this research. Pulmonary sarcoidosis and tuberculosis are the most difficult to distinguish in clinical practice. Even if they have the gold standard of pathological diagnosis, it still needs experienced pathologists to give a definite diagnosis. Sometimes, for the small tissue, it can interfere with the definitive diagnosis of sarcoidosis. At the same time, the treatments for these two diseases are completely opposite. For treatment of sarcoidosis requires glucocorticoids, but the use of glucocorticoids for tuberculosis often causes the spread of disease. More significantly, the mediastinal lymph node of healthy people is not enlarged, and it is not available or difficult to obtain, which also violates ethical principles. Therefore, we used lymphatic tissue of tuberculosis as control group in this study.
There are still several limitations in this study. First of all, in the section of clinical specimen verification, the sample size we included is still small. Although immunohistochemistry and FISH have obtained positive results, a large sample study is still needed to further corroborate our results. Secondly, the specific regulatory mechanism of miR-126 in the pathogenesis and development of sarcoidosis still needs to be explored in vitro and in vivo. Finally, the incidence of sarcoidosis in different races is slightly different. In this study, FIGURE 8 | FISH methods were conducted to detect the expression levels of hsa-miR-126, hsa-miR-20b, as well as hsa-miR-106a in mediastinal lymph node tissues of pulmonary sarcoidosis (SP group, n = 12, 400×) and tuberculosis (CRTL group, n = 13, 400×). (A-C) hsa-miR-126 (+); (D-F) hsa-miR-20b (+); (G-I): hsa-miR-106a (+). Data have been represented as mean ± SD, **p < 0.01. MFI: mean fluorescence; ns: no significant.
we only included Asian individuals. This may bring certain limitations to our results.
From our results, we found that VEGFA and NR3C1 were overexpressed in the pathological tissue of sarcoidosis, which was mainly regulated by hsa-miR-126. To our knowledge, miRNAs are small non-coding RNAs that function as negative gene expression regulators. However, emerging evidences have shown that its function was switching from repression to activation, and its gene-activation function can up-regulate translation in nuclear (Vasudevan et al., 2007;Xiao et al., 2017). The GO annotation of our study comprised the molecular function on current biological knowledge of miRNA in sarcoidosis, which showed nucleic acid binding transcription factor activity has a pronounced gene count. This novel finding was also worth exploring for the study of the role of miRNAs in the regulation of sarcoidosis.
In summary, this work sheds further light on understanding the pathogenesis of sarcoidosis and might indicate novel diagnostic markers and treatment targets in the clinical practice of sarcoidosis. Nevertheless, further studies should be conducted to clarify the role of miRNAs and their related pathways in the pathogenesis of sarcoidosis.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Biomedical Ethics Committee of Second Affiliated Hospital of Xi'an Jiaotong University. The patients/participants provided their written informed consent to participate in this study.