Bioinformatics analysis of the microRNA-mRNA network in sebaceous gland carcinoma of the eyelid.

Sebaceous gland carcinoma (SGC) of the eyelid is an uncommon aggressive tumor with a relatively high rate of local recurrence and a poor prognosis following metastasis. However, the molecular mechanisms underlying the pathogenesis of SGC remain unclear. The purpose of the present study was to clarify microRNA (miRNA) expression profiles in SGC and to explore novel miRNA‑mRNA networks of SGC. A small RNA‑sequencing analysis was performed to identify miRNAs differentially expressed between SGC and sebaceous adenoma control samples. Bioinformatics analyses were conducted to reveal biological functions, canonical pathways and molecular interaction networks using integrated miRNA‑mRNA datasets, including mRNA expression profiles of SGC from our previous study. The present results demonstrated that 16 upregulated miRNAs and 516 downregulated mRNAs were associated with loss of lipid metabolism function and enriched in cholesterol biosynthesis pathways. By contrast, 29 downregulated miRNAs and 194 upregulated mRNAs were mainly associated with the promotion of cell survival and proliferation in addition to enrichment of DNA damage‑induced cell cycle‑regulation pathways. Furthermore, network analyses revealed that the upregulated miRNAs, miR‑130a‑3p and miR‑939‑5p, and the downregulated miRNAs, miR‑146a‑5p, miR‑149‑3p, miR‑193a‑3p, miR‑195‑5p and miR‑4671‑3p, could be upstream regulators related to these functional changes of SGC. These results improved the understanding of molecular mechanisms of SGC and may help to improve the diagnosis of SGC.

Introduction

Sebaceous gland carcinoma (SGC) of the eyelid is a highly malignant tumor that most frequently arises from the meibomian gland and Zeis gland in the periocular region (1). According to clinical reports, it is uncommon in Caucasians populations, accounting for less than a few percent of cases of malignant eyelid tumors, but relatively common in Asian populations, accounting for ~30% of malignant eyelid tumors (2,3). Early diagnosis of SGC is difficult because the disease can mimic benign inflammatory conditions such as chalazion, unilateral conjunctivitis, blepharitis, tarsitis and blepharoconjunctivitis (1,4). SGC tends to be histopathologically misdiagnosed as squamous cell carcinoma or basal cell carcinoma (5,6), but immunohistological stains for adipophilin (7,8) and androgen receptor (9,10) and Oil Red O staining can be helpful to confirm a diagnosis. First-line treatments are surgery and cryotherapy followed by chemotherapy and radiotherapy (11), but some patients have poor prognosis. Delayed diagnosis may result in metastasis to lymph nodes and other organs, leading to metastasis-related mortality in ~6–9% of cases (12–14). Therefore, early and accurate diagnostic markers unique to SGC are needed to improve the prognosis.

The molecular mechanisms underlying the pathogenesis and progression of SGC remain to be fully elucidated. As is the case with other cancers, most SGC have point mutations in the p53 tumor suppressor gene (15) and overexpressed anti-apoptotic proteins including X-linked inhibitor of apoptosis (XIAP) (16) and BAG cochaperone 3 (BAG3) (17). High expression of growth factor receptors such as vascular endothelial growth factor receptor-2, epidermal growth factor receptor, and platelet-derived growth factor receptor are also known as clinicopathologic features of SGC (18). Moreover, expression levels of prognosis factors such as zinc finger E-box binding homeobox 2 (ZEB2) (19), human epidermal growth factor receptor 2 (20) and aldehyde dehydrogenase 1 (21) were known to be used for the prognosis prediction. In our previous study, we revealed the mRNA expression profiles of SGC and identified the gene network consisting of cell cycle related genes including cyclin dependent kinase inhibitor 2A (CDKN2A), cyclin dependent kinase 1 (CDK1) and cyclin E1 (CCNE1) (22). To date, although a number of studies have been conducted to explore novel therapeutic targets, no specific protein expression patterns have been identified for either primary or metastatic lesions of SGC (18).

MicroRNAs (miRNAs) are small non-coding RNAs that bind to complementary sequences of multiple target mRNAs, resulting in post-transcriptional inhibition of gene expression. In humans, more than 2,000 miRNAs controlling complex cellular processes such as proliferation, apoptosis, development and differentiation have been identified (23). In many cancers, miRNAs play roles as upstream regulators of tumorigenesis by contributing to alterations in the gene expression of oncogenes and tumor suppressor genes (24,25) and exploring miRNA-mRNA interactions will thus be critically important to improve our understanding of pathogenesis. Expression patterns of miRNAs vary critically depending on the types of cancers, and a number of clinical trials are currently underway to examine the use of circulating miRNAs as molecular biomarkers for cancer diagnosis (26,27). A few previous studies have shown expression changes in only a limited number of miRNA in SGC samples (28–30), and whole picture of the miRNA-mRNA network of SGC are not fully understood.

In the present study, a small RNA-sequencing analysis were performed to reveal the miRNA expression profiles of SGC and to identify differentially expressed miRNAs common to the tumor samples from three patients with SGC compared to a sebaceous adenoma control sample. In addition, we conducted integrated bioinformatics analyses to identify biological functions, canonical pathways and miRNA-mRNA networks of SGC using the data of mRNA expression profiles obtained from the same tumor sample sets in our previous study (22).

Materials and methods

Patient and tissue samples

This study was performed with the approval of the internal review board of the University of Toyama (no. 27-51), and the procedures conformed to the tenets of the World Medical Association's Declaration of Helsinki. Written informed consent was obtained from all patients prior to enrollment in the present study. Tissues were obtained by surgical excision of tumors from four patients: a 74-year-old woman (case 1), an 83-year-old woman (case 2) and a 58-year-old woman (case 3) with SGC of the eyelid; and a 92-year-old man with sebaceous adenoma of the eyelid (control) for comparison. Tissue samples were immediately frozen and stored at −80°C after sampling for RNA extraction.

RNA extraction and quality control

Total RNA including miRNA was extracted from tissue samples using a NucleoSpin miRNA kit (Macherey-Nagel GmbH & Co.) following the manufacturer's instructions. The quality and quantity of the miRNA were analyzed using a Bioanalyzer 2100 with an RNA 6000 Nano kit (Agilent Technologies) (31).

Library preparation and small RNA-sequencing

Small RNA libraries were prepared using a NEBNext Multiplex Small RNA Library Prep Set for Illumina Set 1 (New England BioLabs). In brief, 1 µg of total RNA per sample was ligated with 3′ and 5′ adaptors and reverse transcribed into first-strand cDNA. Each library was labeled with indexed primers by 15 cycles of PCR amplification and cleaned up using a QIAquick PCR purification kit (Qiagen). Appropriate fractions of 140–150 bp were size-selected by polyacrylamide electrophoresis on the Novex TBE PAGE gel 6% (Invitrogen; Thermo Fisher Scientific, Inc.) and then the purity and concentration were checked using a Bioanalyzer 2100 with a High Sensitivity DNA kit (Agilent Technologies). The pooled libraries were sequenced (2×150 bp) on the HiSeq X Ten platform (Illumina) by Genewiz, Inc. All low sequence data analyzed in this study were deposited in the DNA Data Bank of Japan database under the accession number DRA009187.

Row read data processing

The first 50 bp sequences were extracted from the raw 150 bp sequence reads using Seqkit. Adaptor sequences were trimmed from 50 bp reads using cutadapt. Low-quality (less than Q20) and short-length (less than 10 bp) sequences were removed from processed reads using the FASTX-tool kit. The filtered reads were mapped with hg19, and miRNA annotation was performed using Strand NGS Ver. 3.3.

Microarray data and miRNA-mRNA interaction analyses

In our previous study, mRNA expression profile data of the same sample set were obtained using the GeneChip system with Clariom S human arrays (Affymetrix) (22). Briefly, the raw intensity data (Gene Expression Omnibus; accession no. GSE125582) were normalized using GeneSpring GX 14.9 software (Agilent Technologies). To examine the molecular functions and interaction networks of differentially expressed miRNA and mRNA, the combined data from our present and previous studies were analyzed using Ingenuity Pathways Analysis (IPA) software (Ingenuity Systems).

Results

Identification of differentially expressed miRNAs and mRNAs

To reveal the miRNA expression profiles of SGC samples, a total of 280,241,626 raw reads were obtained in this study, including at least 20,000,000 reads for each sample. Raw sequencing reads were quality checked, and the low-quality sequences and adaptors were removed; the reads were then aligned against the human miRBase using Strand NGS. We obtained read counts of over 2,600 miRNAs, and then identified miRNAs that were at least 2.0-fold differentially expressed compared with the control sample. As shown in the Venn-diagrams in Fig. 1, 16 upregulated and 29 downregulated miRNAs were common to all three tumor samples; these 45 miRNAs are listed in Tables SI and SII. Hierarchical clustering showed that there were distinct expression profiles of miRNAs in SGC and the control sebaceous adenoma sample (Fig. 2). Similarly, 194 upregulated and 516 downregulated mRNAs with at least 2.0-fold change were also identified from our previous study (22).

Figure 1.

Identification of differentially expressed miRNAs in sebaceous gland carcinoma samples. The Venn diagrams show the number of (A) upregulated or (B) downregulated miRNAs in three tumor samples compared with the sebaceous adenoma control sample. miRNA, microRNA.

Figure 2.

Hierarchical clustering representing expression patterns of the microRNAs that were differentially expressed between the three sebaceous gland carcinoma tumor samples and the sebaceous adenoma control sample.

Functional analyses of differentially expressed miRNAs and mRNAs

To explore the biological functions and canonical pathways involving the miRNAs and mRNAs that were differentially expressed between the SGC and control samples, two integrated miRNA-mRNA data sets including 16 upregulated miRNAs with 516 downregulated mRNAs and 29 downregulated miRNAs with 194 upregulated mRNAs were analyzed by IPA software. The top 5 biological functions with positive z-scores in the differentially expressed miRNAs and mRNAs are summarized in Tables I and II. Most biological functions that were significantly enriched in the 16 upregulated miRNAs with 516 downregulated mRNAs were related to the decreasing of the lipid metabolism (i.e., ‘Synthesis of lipid’ and ‘Fatty acid metabolism’). In contrast, 29 downregulated miRNAs with 194 upregulated mRNAs were associated with the increasing of the cell survival and proliferation (i.e., ‘Cell viability of tumor cell lines’, ‘Cell viability’ and ‘Cell proliferation of tumor cell lines’). As shown in Figs. S1 and S2, the top 10 canonical pathways that were differentially activated or suppressed in two integrated miRNA-mRNA data sets are presented, respectively. Multiple annotations associated with cholesterol biosynthesis (i.e., ‘Superpathway of Cholesterol Biosynthesis’ and ‘Cholesterol Biosynthesis I’ to ‘Cholesterol Biosynthesis III’) were redundantly listed in 16 upregulated miRNAs with 516 downregulated mRNAs. Additionally, 29 downregulated miRNAs with 194 upregulated mRNAs were mainly involved in DNA damage-induced cell cycle regulation pathways (i.e., ‘DNA damage-induced 14-3-3δ Signaling’, ‘GADD45 Signaling’ and ‘Cell Cycle: G2/M DNA Damage Checkpoint Regulation’) with no activity pattern available.

Construction of molecular interaction networks of differentially expressed miRNAs and mRNAs

To further understand the regulatory interaction of differentially expressed miRNAs and mRNAs in SGC samples, target prediction analyses were conducted based on the miRNA target filter tool on the IPA software. We identified two miRNA-mRNA networks including 7 upregulated miRNAs that downregulated 8 mRNAs related to decreasing of the synthesis of lipid (Fig. 3) and 5 downregulated miRNAs that upregulated 9 mRNAs related to increasing of the cell proliferation of tumor cell lines (Fig. 4). Based on the results shown in Fig. 3, miR-130a-3p and miR-939-5p were identified as the key hub nodes connected with 3 target mRNAs in the network. In addition, the results shown in Fig. 4 demonstrated that downregulation of miR-146a-5p, miR-149-3p, miR-193a-3p, miR-195-5p and miR-4671-3p played regulatory roles in the promotion of cell proliferation.

Figure 3.

Molecular networks among 7 upregulated miRs and 8 downregulated mRNAs associated with the function of synthesis of lipids. Up- and downregulated molecules are presented as red and green nodes, respectively. miR, microRNA.

Figure 4.

Molecular networks among 5 downregulated miRs and 9 upregulated mRNAs associated with the function of cell proliferation of tumor cell lines. Up- and downregulated molecules are represented as red and green nodes, respectively. miR, microRNA

Discussion

SGC of the eyelid is a rare aggressive tumor with a relatively high rate of metastasis and mortality. One of the treatments for SGC is a complete surgical resection, but the disease occasionally recurs with poor prognosis (6,32). The pathogenesis of SGC remains unclear; therefore, a detailed understanding of the molecular mechanisms will be crucial for improvement of the disease diagnosis, treatment and prognosis. In the present study, we determined the relevant miRNA expression profiles and identified 16 miRNAs that were upregulated and 29 miRNAs that were downregulated in SGC samples compared with the control sebaceous adenoma sample. We then explored the biological functions and canonical pathways and miRNA-mRNA networks related to the clinicopathological characteristics of SGC using integrated miRNA-mRNA data sets. To the best of our knowledge, there are only three previous studies about miRNA expression of SGC. In the first, Bhardwaj et al (28) showed that underexpression of miR-200c and miR-141 were correlated with clinicopathological parameters in SGC, but in our present study we did not observe the downregulation of either of these miRNAs. In the second, Bladen et al (29) demonstrated miRNA expression profiles of SGC using miRNA arrays with 800 probe sets and identified overexpression of miR-16-5p and miR-34a-5p, which were downregulated in the present study. These inconsistencies might be attributable to the control samples used in each experiment: the previous studies used normal tissues (adjacent normal epidermis and tarsal plate) as the control samples while we used sebaceous adenoma. In the third study, Tetzlaff et al (30) examined the expression of 387 miRNAs in SGC by real-time polymerase chain reaction techniques and found that miR-486-5p and miR-184 were overexpressed, while miR-211 and miR-195 were downregulated. They used formalin-fixed paraffin-embedded tissue of sebaceous adenoma as the control samples, and their finding that miR-195 was downregulated was in agreement with our present results. Unlike these previous studies, our present analysis revealed comprehensive expression profiles of over 2,600 miRNAs in SGC using next-generation sequence techniques, and thus our results could provide novel findings of molecular mechanisms of SGC.

One of the interesting findings of this study is that the 16 upregulated miRNAs and 516 downregulated mRNAs in SGC samples were highly associated with the downregulation of lipid metabolism functions and enriched in cholesterol biosynthesis pathways. The main origins of SGC tumors, the meibomian gland and Zeis gland, produce oily substances to protect the periocular regions. As mentioned above, lipid accumulation in the cytoplasm, which is detectable by immunohistochemical staining for adipophilin and Oil Red O staining, is a practical pathological marker of SGC (7,8), which supports the idea that a malfunction in lipid metabolism of the sebaceous glands is involved in the pathogenesis of SGC. It is noted that the gene expressions of both thyroid hormone responsive spot14 (THRSP) and MID1 interacting protein 1 (MID1IP1), a ligand/receptor pair that regulates fatty acid synthesis in non-hepatic cells (33) and lipogenic cancer cells (34,35), were downregulated in this study. In addition, downregulated genes including low-density lipoprotein receptor (LDLR) and glycerol-3-phosphate dehydrogenase 1 (GPD1) were regulatory factors in the synthesis of cholesterol and triglycerides (36,37). Our network demonstrated that miR-130a-3p and miR-939-5p were upstream regulators controlling the expression of these genes related to lipid metabolism. In particular, miR-130 suppressed adipogenesis in human adipocytes in association with a decrease in the gene expression of peroxisome proliferator-activated receptor γ (PPARγ) (38), implying that upregulation of miR-130a-3p could be an index marker related to loss of lipid metabolism functions.

On the other hand, our results also indicated that functions related to cell survival and proliferation were activated in the 29 downregulated miRNAs and 194 upregulated mRNAs in SGC samples. Considering that DNA damage-induced cell cycle regulation pathways were also significantly enriched in the data set, there should be abnormalities in the G2/M cell cycle checkpoint resulting in cell-cycle progression in SGC. Several studies previously showed that high expressions of cell cycle regulatory proteins including p21, p27, cyclin E and p16 (39,40), but hypermethylation of promoter region of CDKN2A were found in half of the cases of SGC tumor (41). Our network showed that upregulation of cell cycle-related genes, including CCNE1, CCNE2 and cyclin dependent kinase inhibitor 3 (CDKN3), was caused by downregulation of miR146a-5p, miR-195-5p and miR-4671-3p. Most importantly, miR-195 is also known to inhibit cell proliferation in association with a decrease in protein expression of cyclin D1 in human cervical cancer cells (42,43), human glioma cells (44) and squamous cell lung cancer (45). In addition, miR-149 acts as a tumor suppressor miRNA controlling cell proliferation and invasion in medullary thyroid carcinoma (46) and renal cell carcinoma (47), and it also plays important roles in regulating the expression of multiple genes in SGC. Overexpression and point mutation of the p53 gene were detected in two-thirds of SGC samples (15,39), suggesting that dysregulation of the cell cycle with downregulation of these miRNA was one of the critical mechanisms in the tumorigenesis and development of SGC tumors.

In conclusion, the present study provides the first comprehensive description of the differentially expressed miRNAs and miRNA-mRNA interaction networks in SGC. We also identified several changes in the expression of miRNAs that control important functional alterations in SGC, including loss of lipid metabolism and promotion of the cell proliferation. These results could improve our understanding of the pathophysiological mechanisms of SGC and provide novel clues for earlier and more accurate diagnosis. Further studies will be needed to confirm the functional roles of these miRNA-mRNA networks in the pathogenesis of SGC.

Table I.

Top 5 biological functions of 16 upregulated microRNAs with 516 downregulated mRNAs.

Functional annotation	P-value	Predicted activation state	Activation z-score	No. of molecules
Transport of molecule	7.84×10−4	Decreased	−3.421	53
Synthesis of lipid	2.88×10−18	Decreased	−2.881	64
Proliferation of pancreatic cancer cell lines	1.23×10−2	Decreased	−2.781	10
Neoplasia of tumor cell lines	8.15×10−3	Decreased	−2.774	8
Fatty acid metabolism	4.07×10−9	Decreased	−2.666	35

Activation z-score, >2.0 or <-2.0 was considered to be significantly predictive.

Table II.

Top 5 biological functions of 29 downregulated microRNAs with 194 upregulated mRNAs.

Functional annotation	P-value	Predicted activation state	Activation z-score	No. of molecules
Cell viability of tumor cell lines	4.41×10−3	Increased	3.810	22
Cell viability	4.45×10−3	Increased	3.732	25
Cell proliferation of tumor cell lines	6.65×10−3	Increased	2.993	35
Cell movement of endothelial cells	2.50×10−3	Increased	2.375	10
Cell proliferation of breast cancer cell lines	5.66×10−4	–	1.843	16

Activation z-score, >2.0 or <-2.0 was considered to be significantly predictive.