Genome-Wide Expression Difference of MicroRNAs in Basal Cell Carcinoma

Wei, Hai-Peng; Zhan, Song; Zhu, Qing-An; Chen, Zhen-Juan; Feng, Xian; Chen, Jun-Yuan; Zhang, Qi-Lin; Zhao, Jingjie; Meng, Lingzhang

doi:https://doi.org/10.1155/2021/7223500

Journal of Immunology Research

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusions Abbreviations Data Availability Conflicts of Interest Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Special Issue

Extracellular Matrix and Immune Niches in Human Disease 2021

View this Special Issue

Research Article | Open Access

Volume 2021 | Article ID 7223500 | https://doi.org/10.1155/2021/7223500

Genome-Wide Expression Difference of MicroRNAs in Basal Cell Carcinoma

Hai-Peng Wei,¹Song Zhan,¹Qing-An Zhu,¹Zhen-Juan Chen,¹Xian Feng,¹Jun-Yuan Chen,²Qi-Lin Zhang,^1,3Jingjie Zhao ,⁴and Lingzhang Meng⁵

Academic Editor: Zenghui Teng

Received13 May 2021

Accepted22 Jul 2021

Published04 Aug 2021

Abstract

Distinct expression of the miRNAs has rarely been explored in basal cell carcinoma (BCC) of skin, and the regulatory role of miRNAs in BCC development remains quite opaque. Here, we collected control tissues from adjacent noncancerous skin (; control group) and tissues at tumor centers from patients with cheek BCC (; BCC group) using punch biopsies. After six small RNA sequencing- (sRNA-seq-) based miRNA expression profiles were generated for both BCC and controls, including three biological replicates, we conducted comparative analysis on the sRNA-seq dataset, discovering 181 differentially expressed miRNAs (DEMs) out of the 1,873 miRNAs in BCCs. In order to validate the sRNA-seq data, expression of 15 randomly selected DEMs was measured using the TaqMan probe-based quantitative real-time PCR. Functional analysis of predicted target genes of DEMs in BCCs shows that these miRNAs are primarily involved in various types of cancers, immune response, epithelial growth, and morphogenesis, as well as energy production and metabolism, indicating that BCC development is caused, at least in part, by changes in miRNA regulation for biological and disease processes. In particular, the “basal cell carcinoma pathways” were found to be enriched by predicted DEM targets, and regulatory relationships between DEMs and their targeted genes in this pathway were further uncovered. These results revealed the association between BCCs and abundant miRNA molecules that regulate target genes, functional modules, and signaling pathways in carcinogenesis.

1. Introduction

Basal cell carcinoma (BCC) is the most ordinary human solid tumor in the skin, accounting for approximately 80% of nonmelanoma skin cancers [1, 2]. Most BCC patients usually exhibit translucent lesions in superficial and easily develop symptoms such as bleeding and crusting. In general, BCC primarily occurs in the scalp, cheek, forearm, thorax, popliteal, and preauricular regions, with the cheek being the most common region [3, 4]. From the first case described by Krompecher in 1900, several investigations have been conducted to uncover the molecular pathogenesis of BCC. For example, several tumor-promoting signal transduction pathways have been reported, such as the sonic hedgehog (Shh) and mitogen-activated protein kinase/extracellular signal-regulated kinase (Mapk/Erk) pathways [1, 2, 5]. Additionally, several genes/proteins have also been reported to be associated with the development of BCC. In particular, the patched (Ptch1) receptor binding to the Shh ligand casts impacts on the smoothened protein (Smo) on the cell membrane, which results in the suppression of the zinc finger transcription factor gene (GLI1) [2]. A mechanism as such counteracts cell-cycle arrest and leads to the development of skin hair follicle BCCs [6–8]. Moreover, functional loss of the Ptch1 receptor caused by gene mutations has been detected in approximately 70% of BCCs [2, 8]. In addition, gain-of-function mutations of SMOH have been found in approximately 15% of BCCs [2]. These data indicates that variations in dimensions of expression and sequence (including pathways) of genes are vital for the progress of BCCs.

MicroRNAs (miRNAs) are a rather novel type of short RNAs which are known to be involved in posttranscriptional gene expression [9]. As it has been evaluated that 30–60% of the entire human gene expression is controlled by miRNAs, a possible function of miRNAs in the molecular etiopathogenesis of BCC needs to be investigated [1]. As a class, miRNAs function across numerous processes of the skin physiologically and pathophysiologically, as well as across other human tissues [10, 11]. The miRNA maturing enzymes DICER and DROSHA, partial miRNA microprocessor complex, along with the miRNA effector RNA-caused silencing complex have been discovered to be dysregulated in epithelial skin cancers which includes BCC [12]. Currently, several studies have revealed miRNA expression profiles from BCC and nonlesional epithelial skin samples of clinical patients by using small RNA sequencing (sRNA-seq) and miRNA microarrays [1, 2]. These miRNA-based researches have been primarily focused on identifying the differentially expressed miRNAs (DEMs) in BCCs by comparing with control tissues, but the molecular etiopathogenesis of BCCs with regard to miRNAs has not been broadly investigated, especially on the regulatory functions (i.e., target gene function and signaling pathways) and relationships (i.e., networks of miRNA-targeted genes and miRNA-targeted genes-pathways) of miRNAs. Hence, it is essential to systematically recognize the miRNAs participating in BCCs and analyze the potential genes and pathways that are possibly regulated by them. Though numerous reports have emphasized the biological/pathological role of miRNAs played in the development of BCC, few studies have been found to illustrate the landscape of miRNA in BCC biopsies.

In this study, we identified miRNAs in BCCs using sRNA-seq and bioinformatics approaches, while discovering a panel of DEMs in BCC tissues by comparing genome-wide miRNA transcriptional profiles of BCC with control groups and by conducting TaqMan quantitative real-time PCR (qPCR). Subsequently, based on the reference gene sets of humans, the primary regulatory functions and relationships of these miRNAs were further elucidated. The regulatory roles of miRNAs have been unveiled by our research, and their fundamental biological role will facilitate the comprehension of the molecular pathogenesis of BCCs.

2. Materials and Methods

2.1. BCC Materials

Our research was granted by the Ethical Review Board of the Beihai People’s Hospital Guangxi Province, China (Approval No. Beiyi (lunshen) 2020003). As we adopted the dermatological surgery section, the present research was conducted according to Ethical Principles for Medical Research Involving Human Sub in the light of the Declaration of Helsinki (version 2013). Every subject of this research has provided written informed consent.

The BCC and control samples in the cheeks of five patients were collected as per abovementioned methods [1, 2, 4]. To be specific, during the micrographic controlled surgery with cold steel by regional anaesthesia, 4 mm punch biopsies were both taken from cancer cores (BCC groups) and adjacent noncancerous epithelial skin (control groups). Each sample was conducted three times independently ( patients per collection; patients across three times) with three biological replicates. Overall, 15 patients ( female and male; mean age 67.4 years) with cheek BCC were included in our research (Table 1). The specimens were immediately put in RNAlater (Qiagen, Hilden, Germany) and reserved until further use at -80°C.

2.2. RNA Extraction, Library Construction, and Small RNA Sequencing

Entire RNA was extracted from each BCC and control specimens and treated with the TRIzol reagent (Invitrogen, CA, USA) as per the manual from the supplier. Residual genomic DNA was extracted via RNase-free DNase Set (Qiagen, Germany), as per the manual from the supplier. The purity was evaluated by using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, CA, USA), while the RNA structural completeness was validated using agarose gel electrophoresis at 1.5%. To eliminate biological variations by different levels of gene expression in each biological replicate, the entire RNA from each sample was diluted to 200 ng/μL using RNase-free water. Afterwards, an equal volume of the RNA from all five BCC samples of each biological replicate was pooled together to develop a sRNA-seq library, and RNA from five control samples was collected likewise. According to the MGIEasy Small RNA Library Prep Kit V2.0 (Illumina, CA, USA) protocol, approximately 2 μg of total pooled RNA was utilized to prepare a small RNA library. To generate three biological replicates, the experiments were independently repeated three times. Under rigorous quality control standards, we carried out single-end sequencing (45 bp) of six sRNA libraries on an BGISEQ-500 at the Beijing Genomics Institute (BGI, Shenzhen, China).

2.3. Preprocessing Analysis, Annotation of miRNAs, and Differentially Expressed miRNAs

The raw tags were filtered by eliminating defective sequences, sequences with poly-A tail and adaptors, while the remaining sequences were retained as clean tags. The length distributions of clean sequences were then summarized, followed by the removal of the highly repetitive clean tags via the RepeatMasker tool (http://www.repeatmasker.org/). In order to avoid nonhuman contamination, the clean tags were mapped to the human reference genome (NCBI version: GRCh38.p13) using the Bowtie2 [13] tool. Subsequently, these clean tags were searched within the GenBank and Rfam databases by using Bowtie2 to determine the annotation information of non-miRNA sRNAs (e.g., rRNAs, tRNAs, and snoRNAs). Eventually, the non-miRNA sequences were removed for the following analysis.

To identify known miRNAs, the obtained clean tags were then annotated using the miRBase 22.1 database (http://www.mirbase.org/), while the completely matched (no mismatch base) clean tags with mature miRNAs in the miRBase were recognized as miRNAs. To acquire dependable outcomes by strict standards, miRNAs that were shared by the 3 biological duplicates were reserved as the miRNA set for control and BCC samples, respectively. In order to calculate the mapped tags and quantify the miRNA expression in downstream analyses, all clean tags were aligned to the miRNA sequences using Bowtie2.

The expression of miRNAs across the entire 6 libraries was quantified as per the obtained count of mapped tags for each miRNA via the quantification method of transcripts per million kilobase (TPM) values [14]. Pearson’s correlation () between two libraries was computed to determine the reproducibility within the 3 biological duplicates. Certain miRNAs (more than two mapped tags) were either not detected or only detected in BCC groups, which suggests that the expression of these miRNAs was affected by the BCC. In order to include these potential candidates for the BCC process, the union of recognized miRNAs in the control and BCC samples was utilized to identify DEMs. DESeq v1.36.0 [15], an R program suitable for the identification of differentially expressed genes from high-throughput sequencing data, helped identify DEMs between the control and BCC samples. The miRNAs with and values (Wald test in DESeq) rectified by a were recognized as DEMs.

2.4. Prediction and Functional Enrichment Analysis of miRNA Target Genes

Available human gene sets downloaded from NCBI were utilized as the reference geneset to predict miRNA targets. Furthermore, gene functions were annotated by searching within the NCBI nonredundant protein (NR) databases by using a Blastx tool with default parameters. To determine the controlling role of DEMs in BCCs, the potential target genes of miRNAs were predicted by adopting TargetScan v7.2 [16] and miRanda both of which are appropriate for forecasting the target genes of mammalian miRNAs [17]. The forecast shared by both methods was reserved as dependable targets.

In order to elucidate the biological role of DEM-regulated target genes and signaling pathways, the functional analyses of miRNA target genes were implemented by using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) v6.8 [18, 19]. By adopting the Functional Annotation Clustering tool, Gene Ontology (GO) terms, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, the annotated functional terms of DEMs were clustered, while each clustered set with an FDR- was taken as remarkably enriched. Furthermore, the superfluous GO terms were thoroughly removed using GO Trimming [20].

2.5. Validation of sRNA-seq Data Using TaqMan qPCR

In order to validate the reliability of sRNA-seq results, relative quantification of the TaqMan qPCR approach was utilized to determine the expression of 15 casually picked DEMs, involving 8 upregulated and 7 downregulated DEMs, by utilizing the identical total RNA specimens which were employed in sRNA-seq. The TaqMan® Advanced miRNA cDNA Synthesis Kit (Applied Biosystems, Darmstadt, Germany) was adopted to conduct the poly(A) tailing reaction, adaptor ligation reaction, reverse transcription (RT) reaction, and miR-Amp reaction (miRNA preamplification) in order, as per the manufacturer’s instructions. Afterwards, RNase-free water was added in order to dilute miR-Amp products into 10 ng/μL. The TaqMan® Advanced miRNA Assays kits for qPCR reaction were made-to-order for each miRNA, including the 2x TaqMan Fast Advanced Master Mix, TaqMan probes, and PCR primers. The U6 snRNA was selected as a reference gene [1], while product information for each miRNA is listed in Supplementary Table 1. Each qPCR reaction system at an eventual reaction volume of 10 μL was implemented as per the manufacturer’s instructions of kits in an ABI 7500 Fast Real-Time PCR System (Applied Biosystems, Darmstadt, Germany). The following experimental instructions were adopted for qPCR: enzyme activation at 95°C for 20 s, subsequently 40 cycles of denaturation at 95°C for 3 s, and annealing and extension at 60°C for 30 s. A negative control without template was included, while 3 technical duplicates and 3 biological duplicates were conducted for the qPCR reaction of each miRNA.

Relative expression was computed by the 2^-ΔΔCT approach [21] in the MedCalc v19.6.1 statistics software, the outcomes of which were described as the . By using significance, calculation of the Pearson correlation coefficient between sRNA-seq and qPCR results was also implemented via this statistical software.

3. Results

3.1. Overview of Sequencing Data

For the 6 sRNA-seq libraries, about 88 million and 87 million raw tags were attained by sequencing the control and BCC samples, separately. The percentage of Q20 tags in the raw data of each library was greater than 99%, while the clean tags accounted for more than 80% of the raw tags in every sRNA-seq library after data filtering. Afterwards, the obtained clean tags were uploaded onto the NCBI BioProject (accession No. PRJNA597543, weblink: https://www.ncbi.nlm.nih.gov/bioproject/?term=597543). The percentage (approximately 92-96%) of mapped clean tags to the reference genome was almost identical among the six sRNA-seq libraries, which indicates no bias in data generation among libraries. After the annotation of various sRNA types, most clean tags were confirmed to be miRNAs, accounting for 60-65% of all clean tags in controls and 78-83% of all the clean tags in BCCs. Detailed information regarding the miRNAs is generalized in Table 2, while the length distribution of the most abundant clean tags is primarily concentrated from 21-23 bp in each library (Supplementary Figure 1). These results are consistent with the known size of sRNAs [22], indicating the reliability of sRNA-seq.

Based on the expression (TPM values) of each identified sRNA molecule, the Pearson correlation coefficient was more than 0.9 for in pair contrast of biological replicate libraries from the identical treatment (Figure 1(a)). The outcome suggests sampling effectiveness and reliable expression quantification within this study.

(a)

(b)

Figure 1

Overview of sRNA sequencing results. (a) Heatmap of Pearson correlation coefficients across diverse specimens. Greater values and darker-blue squares refer to higher similarity between the two sequencing libraries. BCC_1 to BCC_3 indicate sRNA-seq libraries of BCCs across three biological replicates. Control_1 to Control_3 indicate sRNA-seq libraries of controls across three biological replicates. (b) Volcano plot shows the distribution of 181 DEMs (BCC vs. control). The DEMs with and were considered significant ones.

3.2. Identification of miRNAs

For the control groups, 1,997, 1,954, and 1,901 miRNAs were identified as known miRNAs across the three biological replicates. In total, 1,534 known miRNAs were shared across these 3 libraries. For the BCC groups, 1,982, 1,999, and 1,987 known miRNAs were also identified across the 3 biological duplicates, 1,656 known miRNAs of which were shared by the entire 3 libraries. Eventually, when we integrated the control and BBC samples, a total of 1,873 known miRNAs (Supplementary Table 2) were identified and used for the identification of additional DEMs.

3.3. Differentially Expressed miRNAs and TaqMan qPCR Validation

Overall, 181 miRNAs were identified to be DEMs (Figure 1(b)) in this study, indicating a systematic identification for DEMs in BCCs. These DEMs accounted for 9.66% of all miRNAs that were identified in this study and comprised 53 upregulated and 128 downregulated miRNAs (Supplementary Table 3). Meanwhile, The miRNAs with the top 10 different expression levels are listed in Table 3. Seven miRNAs (hsa-miR-2113, hsa-miR-1269b, hsa-miR-135a-3p, hsa-miR-551a, hsa-miR-548y, hsa-miR-548i, and hsa-miR-9983-3p) were shown to be downregulated by greater than 16-fold in BCCs, while all top 10 miRNAs were downregulated by greater than 16-fold. The TaqMan qPCR analysis for all 15 DEMs represents a close correlation (Pearson’s correlation , ) in fold changes for DEMs between sRNA-seq and qPCR (Figure 2), which suggests the accuracy of sRNA-seq and DEM analyses. Among the upregulated DEMS, it is known that has-miR-9983-3p, has-miR-5695, and has-miR-941 contribute to the development of nasopharyngeal carcinoma, thyroid carcinoma, and laryngeal squamous cell carcinoma.

3.4. Function of Target Genes Regulated by Differentially Expressed miRNAs

In total, 1,356 DEM target genes that were shared by TargetScan and miRanda analyses were identified. Results from the GO enrichment method suggested that these target genes were significantly enriched within the 39 biological process terms (Supplementary Table 4). The “biological process” functions of DEM target genes were primarily involved in epidermal growth and morphogenesis, immune response, cell regulation and migration, metabolism processes, and skin development.

Moreover, 38 KEGG signaling pathways are remarkably enriched by target genes of DEMs (Table 4). Furthermore, 10 KEGG signaling pathways enriched by target genes of DEMs were identified as typical cancer-related pathways, including basal cell carcinoma (ko05217), breast cancer (ko05224), and hepatocellular carcinoma (ko05225). The remaining KEGG terms are primarily associated with ATP production (i.e., oxidative phosphorylation (ko00190) and citrate cycle (TCA cycle) (ko00020)) and metabolism (e.g., protein digestion and absorption (ko04974), histidine metabolism (ko00340), and mineral absorption (ko04978)).

3.5. Regulatory Roles of DEMs for Basal Cell Carcinoma Pathway

“Basal cell carcinoma” (ko05217) enriched by predicted target genes of DEMs is a typical signaling pathway involved in the development of BCC, and therefore, its gene members have a key role. The DEMs-targeting gene members in this pathway are considered to be key miRNAs (KmiRNAs) for BCCs. In order to further explore the regulatory roles of these KmiRNAs in BCCs, their target gene members in the “basal cell carcinoma” KEGG pathway were further extracted (Supplementary Table 5), according to the miRNA-target relationships predicted above, and regulatory networks between the KmiRNAs and these pathways were constructed (Figure 3). In total, 49 DEMs were predicted to be KmiRNAs, and 24 gene members (marked in a red box in Supplementary Figure 2) in the “basal cell carcinoma” pathway were shown to be regulated by these KmiRNAs. Notably, the majority of KmiRNAs (e.g., hsa-miR-629-3p, hsa-miR-3176, hsa-miR-3194-3p, and hsa-miR-4701-5p) that regulated the smoothened homolog precursor (SMO), an oncogene, are significantly downregulated in BCCs. Conversely, there is a remarkable upregulated expression of most KmiRNAs (e.g., hsa-miR-320b, hsa-miR-30c-1-3p, hsa-miR-30b-3p, hsa-miR-6815-5p, and hsa-miR-1247-3p) that regulate genes encoding tumor suppressors (i.e., PTCH1 and TP53) in BCCs.

In addition, six genes that are involved in tumorigenic processes, such as uncontrolled proliferation, increased survival of cancer cells, and genomic instability, have been predicted to be targets of KmiRNAs. These genes include cyclin-dependent kinase inhibitor 1A (CDKN1A), apoptosis regulator BAX (BAX), DNA damage-binding protein 2 (DDB2), growth arrest and DNA-damage-inducible protein (GADD45), bcl-2 homologous antagonist/killer 1 (BAK1), and DNA polymerase kappa (POLK). Furthermore, seven genes associated with cancer cell invasion were detected in the basal cell carcinoma pathway, including frizzled 4 (FZD4), segment polarity protein dishevelled 3 (DVL3), glycogen synthase kinase 3 beta (GSK3B), axin 1 (AXIN1), adenomatosis polyposis coli protein (APC2), catenin beta 1 (CTNNB1), and transcription factor 7 (TCF7). Moreover, seven genes that play a role in cancer cell proliferation were also found to be regulated by KmiRNAs, including suppressor of fused homolog (SUFU), zinc finger protein GLI3 (GLI3), bone morphogenetic protein 2 (BMP2), hedgehog interacting protein 1 (HIP1), GLI1, patched 2 (PTCH2), and wingless-type MMTV integration site family, member 9 (WNT). Notably, the majority of KmiRNAs that regulate these cancer/invasion/proliferation-promoting genes were found to be significantly downregulated in BCCs. Interestingly, numerous gene members in the basal cell carcinoma pathway were found to be regulated by multiple miRNAs, with the exception of POLK and KIF7; meanwhile, for many KmiRNAs, one of them (e.g., hsa-miR-320b, hsa-miR-18a-3p, and hsa-miR-370-3p) targets more than one gene members.

4. Discussion

For miRNAs, several studies have characterized DEMs in BCCs through the use of sRNA-seq, microarray, and RT-PCR miRNA profiling [1, 2]. However, to date, no studies have explored the regulatory function and relationships of miRNAs in BCCs. In the present research, we identified massive miRNAs in BCCs using sRNA-seq and bioinformatics. Furthermore, we also identified a group of DEMs, which is far more than the number of DEMs that were previously reported in BCCs (e.g., 26 DEMs [1] and 33 upregulated miRNAs [2]). Thus, it is clear that previous studies have largely underestimated the number and types of DEMs in BCCs. In addition, we found that most DEMs that were identified by previous studies were also detected in this study, including hsa-miR-223, hsa-miR-452-5p, and hsa-miR-941. This observation indicates that these miRNAs are common regulators in various types of BCCs and that they are basic miRNAs for the development of BCC. These common DEMs may represent a potential biomarker for early BCC detection.

Overall, we identified 181 miRNAs that were significantly differentially expressed in our cohort, indicating that these miRNAs participate in BCC processes. In particular, upregulated miRNAs were detected among the list of top 10 DEMs in BCC. The regulatory function of several of these DEMs has been reported in previous studies. For example, upregulated expression of hsa-miR-1269b has been found in hepatocellular carcinoma [23] and was considered a prognostic marker for hepatocellular carcinoma [24]. Deregulation of expression of hsa-miR-135a-3p was identified in urothelial carcinoma of the bladder, and this miRNA is a biomarker for predicting survival of patients [25]. Furthermore, hsa-miR-551a was found to significantly affect the prognosis of patients with pancreatic cancer and can be used as a biomarker of tumor prognosis [26]. Previous studies have shown that single nucleotide polymorphisms (SNPs) in the hsa-miR-548 family binding site at the glyceraldehyde-3-phosphate dehydrogenase alter susceptibility to breast cancer [27]. In combination with results from this current study, the hsa-miR-548 family (i.e., hsa-miR-548y and hsa-miR-548i) can regulate cancer processes at both the gene expression and SNP level. The survivin gene (baculoviral IAP repeat containing 5, BIRC5), a novel member of the inhibitor of apoptosis protein family, is only expressed in tumor and embryonic tissues and is tightly associated with the differentiation, proliferation, invasion, and metastasis of tumor cells. The expression of hsa-miR-542-3p was found to be negatively correlated to survivin levels in oral squamous cell carcinoma [28]. This study found upregulated expression of hsa-miR-542-3p in BCCs, which suggests that inhibition of survivin gene expression by increasing apoptosis likely promotes BCCs. For downregulated miRNAs listed in the list of top 15 DEMs in BCC, the expression of hsa-miR-516a-5p has been described as downregulated in the decidua and villus of recurrent miscarriage patients [29]. Previous studies have demonstrated reduced expression of hsa-miR-378i in stage II colon cancer, which affects cell proliferation, cell-cell interaction, and apoptosis [30]. Furthermore, it is worth noting that regulatory functions of several miRNAs (i.e., hsa-miR-2113, hsa-miR-9983-3p, hsa-miR-3619-3p, hsa-miR-550b-2-5p, and hsa-miR-523-3p) in the list of top 15 DEMs in BCC have not yet been explored, particularly in the context of cancer. Thus, despite these miRNAs being highly differentially expressed in BCCs in comparison with normal tissues, the regulatory role of these miRNAs has yet to be fully elucidated. The results presented here provide key miRNA candidates in the development of BCCs.

In order to further explore the global regulatory functions of miRNAs for BCC, functional enrichment analyses of potential mRNAs targets of DEMs are necessary for uncovering miRNA-mediated biological processes and pathways. In this current investigation, the predicted target genes of DEMs were found to be remarkably enriched according to various GO and KEGG terms. Many cases have reported that BBCs arise from abnormal epidermal cell differentiation levels/epidermal diseases (i.e., cyst, hyperplasia, and cell carcinomas) [31, 32]. These GO terms are primarily involved in skin development, epidermal growth, and morphogenesis, which can explain why BCCs are related to an abnormal epidermal epidermis. Furthermore, we identified GO terms that are primarily involved in the immune response and metabolism, which are expected as disturbances in the immune system and metabolism and have been widely demonstrated across tumors such as BCCs [33, 34]. Cell regulation and migration is a key step for BCC developments [35]. In this study, we found that BCCs were also enriched for GO terms that are associated with cell regulation and migration. These DEMs that target regulation and migration of BCC cells can represent potential therapeutic targets to control BCC invasion.

Several enriched KEGG pathways were identified in this study, and they were largely found to be immune-related. These outcomes unveiled partial vital immune-related pathways participating in BCCs. In particular, dysfunction of immune signaling pathways possibly has an impact on cancer risk and features of cancer progression [36]. In this study, we found miRNAs that target these immune signaling pathways, providing evidence for modulation of risk and progression of cancer caused by alterations in miRNA expression, with BCCs as an example. Furthermore, at least 10 significant pathways were primarily related to typical human cancers, such as breast cancer, gastric cancer, and small-cell lung cancer. Thus, the molecular pathogenesis of BCCs may be implicated in that of other cancers, and many cancer-related signaling pathways that are regulated by miRNAs are shared among various cancer types. Previous studies that are involved in the molecular pathogenesis of human breast cancer have also targeted genes of DEMs in tumor tissues that are enriched in pathways involved in many other cancers (i.e., prostate cancer, basal cell carcinoma, and acute myeloid leukemia) [37], further supporting the speculation. Overall, this study is the first to reveal signaling pathway candidates regulated by miRNAs and the relationship between these miRNAs and pathways in BCCs.

Notably, a typical BCC-related pathway “basal cell carcinoma” was enriched in this study, which led us to uncover miRNAs that regulate this pathway and construct regulatory networks of miRNA-gene members in this pathway. In the “basal cell carcinoma pathway” that was enriched in this study, almost all genes associated with these three primary BCC processes (i.e., tumorigenesis, proliferation of cancer cells, and invasion) can be detected as targets of KmiRNAs. These results indicate that miRNAs probably regulate all the key steps in BCC processes, similar to other cancer types including breast cancer [38], prostate cancer [39], and gastric cancer [40]. Interestingly, in general, the role of miRNAs seems to be gene regulation by targeting specific transcripts for the suppression of gene expression. This study has found downregulated expression in most KmiRNAs. Evidence as such suggests that relax/reduction of miRNA-mediated inhibition of genes may play a vital role in BCC development, which needs to be further studied by detecting expression changes in target gene members. The majority of KmiRNAs that have upregulated expression primarily target tumor suppressor genes in the BCC pathway. For example, mutation-based functional loss of PTCH1 receptor has been significantly detected in most BCCs [2, 8]. Upregulated expression of most KmiRNAs that regulate PTCH1 in this study suggests that suppression of PTCH1 also induces the development of BCCs and further demonstrates the function of PTCH1 as a tumor suppressor. In addition, the tumor suppressor gene TP53 is a key regulator of cell division and/or apoptosis, and overexpression of the TP53 protein may suppress its antiapoptotic activity in BCC cells [41]. Therefore, upregulated expression of all KmiRNAs that target TP53 indicates that miRNAs have an effect on BCC development by regulating TP53-mediated apoptosis. This study lacks the evidence from cell culture and/or animal experiments, to show that miRNAs could influence the biological features in BCC; however, these results provide strong correlations between miRNAs and typical BCC-related pathways. Certainly, regulatory relationships of multiple miRNAs-multiple target genes have been found in the basal cell carcinoma pathway, which shows extreme complexity in regulatory networks of miRNAs in BCCs.

5. Conclusions

For vital regulators of BCC molecular etiopathogenesis, the present research provides a dataset of DEMs/miRNA candidates which should be verified and examined functionally. Moreover, we are the first to explore the regulatory roles/functions of miRNAs and regulatory networks of miRNAs-target genes in BCC tissues by using sRNA-seq and bioinformatics. Thus, our research facilitates deeper comprehension of miRNA regulatory mechanisms which function in the BCC molecular etiopathogenesis. We found that certain DEMs (e.g., hsa-miR-2113, hsa-miR-9983-3p, and hsa-miR-3619-3p) cannot be directly linked to the BCC processes, as their functions were not previously reported in cancer and other human diseases, and only differences on expression levels do not completely verify their regulatory roles in BCCs, although those miRNAs indeed displayed a remarkable expression variation in this study. Hence, more researches using sRNA-seq and functional genomics are required to verify their roles.

Abbreviations

BCC:	Basal cell carcinoma
DEMs:	Differentially expressed miRNAs
sRNA-seq:	Small RNA sequencing
qPCR:	TaqMan quantitative real-time PCR
TPM:	Transcripts per million kilobase
NR:	NCBI nonredundant protein
DAVID:	Database for Annotation, Visualization, and Integrated Discovery
GO:	Gene Ontology
KEGG:	Kyoto Encyclopedia of Genes and Genomes
SNPs:	Single nucleotide polymorphisms.

Data Availability

The data used to support the findings of this study is available from the corresponding authors upon reasonable request.

Conflicts of Interest

The authors declare no competing financial interests.

Authors’ Contributions

Q-LZ, LM, and JZ conceived and designed the study. H-PW performed the surgical isolation of carcinoma biopsies from patients and performed data analysis. SZ, Q-AZ, and Z-JC performed RNA isolation and quantitative PCR. XF performed partial data analysis. J-YC composed the manuscript and performed data analysis. Every author has read, commented on, and accepted the manuscript.

Acknowledgments

Our research was supported by the Science and Technology Plan Projects of Beihai City (BKH201884029).

Supplementary Materials

Supplementary Figure 1: length distribution of sRNA tags in six sequencing libraries. Supplementary Figure 2: KEGG pathway of basal cell carcinoma and 24 miRNA-regulated gene members marked in a red box. Supplementary Table 1: information of TaqMan probes used in quantitative PCR. Supplementary Table 2: list of miRNAs identified in this study. Supplementary Table 3: differentially expressed miRNAs (DEMs) in the BCCs when compared with control. Supplementary Table 4: list of the enriched “biological process” GO terms of targeted genes of DEMs between the control and BCC groups. Supplementary Table 5: miRNA-targeted gene analysis in the basal cell carcinoma (ko05217) pathway. u: upregulation; d: downregulation. (Supplementary Materials)

References

M. Sand, M. Skrygan, D. Sand et al., “Expression of microRNAs in basal cell carcinoma,” British Journal of Dermatology, vol. 167, no. 4, pp. 847–855, 2012.
View at: Publisher Site | Google Scholar
M. Sand, F. G. Bechara, T. Gambichler et al., “Next-generation sequencing of the basal cell carcinoma miRNome and a description of novel microRNA candidates under neoadjuvant vismodegib therapy: an integrative molecular and surgical case study,” Annals of Oncology, vol. 27, no. 2, pp. 332–338, 2016.
View at: Publisher Site | Google Scholar
W. Gou, T. Yu, and W. Zhang, “Clinical pathological analysis of 364 cases of basal cell carcinoma,” Journal of Practical Dermatology, vol. 10, no. 5, pp. 280–283, 2017.
View at: Google Scholar
H. P. Wei, S. Zhan, and Q. A. Zhu, “Analysis of clinical characteristics of 113 cases of basal cell carcinoma of skin in Beihai city,” Youjiang Medical Journal, vol. 48, no. 4, pp. 272–275, 2020.
View at: Google Scholar
D. Goppner and M. Leverkus, “Basal cell carcinoma: from the molecular understanding of the pathogenesis to targeted therapy of progressive disease,” Journal of Skin Cancer, vol. 2011, Article ID 650258, 8 pages, 2011.
View at: Publisher Site | Google Scholar
A. E. Oro, K. M. Higgins, Z. Hu, J. M. Bonifas, E. H. Epstein, and M. P. Scott, “Basal cell carcinomas in mice overexpressing sonic hedgehog,” Science, vol. 276, no. 5313, pp. 817–821, 1997.
View at: Publisher Site | Google Scholar
M. Nilsson, A. B. Unden, D. Krause et al., “Induction of basal cell carcinomas and trichoepitheliomas in mice overexpressing GLI-1,” Proceedings of the National Academy of Sciences, vol. 97, no. 7, pp. 3438–3443, 2000.
View at: Publisher Site | Google Scholar
P. Boukamp, “Non-melanoma skin cancer: what drives tumor development and progression?” Carcinogenesis, vol. 26, no. 10, pp. 1657–1667, 2005.
View at: Publisher Site | Google Scholar
W. Filipowicz, S. N. Bhattacharyya, and N. Sonenberg, “Mechanisms of post-transcriptional regulation by microRNAs: are the answers in sight?” Nature Reviews Genetics, vol. 9, no. 2, pp. 102–114, 2008.
View at: Publisher Site | Google Scholar
K. Ono, Y. Kuwabara, and J. Han, “MicroRNAs and cardiovascular diseases,” FEBS Journal, vol. 278, no. 10, pp. 1619–1633, 2011.
View at: Publisher Site | Google Scholar
M. Sand, “Paper 8: microRNAs and the skin: tiny players in the body’s largest organ,” in in MicroRNAs in Malignant Tumors of the Skin: First Steps of Tiny Players in the Skin to a New World of Genomic Medicine, P. D. D. M. Sand, Ed., pp. 211–231, Wiesbaden: Springer Fachmedien Wiesbaden, 2016.
View at: Google Scholar
M. Sand, T. Gambichler, M. Skrygan et al., “Expression levels of the microRNA processing enzymes Drosha and dicer in epithelial skin cancer,” Cancer Investigation, vol. 28, no. 6, pp. 649–653, 2010.
View at: Publisher Site | Google Scholar
B. Langmead and S. L. Salzberg, “Fast gapped-read alignment with Bowtie 2,” Nature Methods, vol. 9, no. 4, pp. 357–359, 2012.
View at: Publisher Site | Google Scholar
Q.-L. Zhang, Q.-H. Zhu, F. Zhang, B. Xu, X.-Q. Wang, and J.-Y. Chen, “Transcriptome-wide analysis of immune-responsive microRNAs against poly (I:C) challenge inBranchiostoma belcheriby deep sequencing and bioinformatics,” Oncotarget, vol. 8, no. 43, pp. 73590–73602, 2017.
View at: Publisher Site | Google Scholar
S. Anders and W. Huber, “Differential expression analysis for sequence count data,” Genome Biology, vol. 11, no. 10, p. R106, 2010.
View at: Publisher Site | Google Scholar
V. Agarwal, G. W. Bell, J. W. Nam, and D. P. Bartel, “Predicting effective microRNA target sites in mammalian mRNAs,” eLife, vol. 4, 2015.
View at: Publisher Site | Google Scholar
D. Betel, M. Wilson, A. Gabow, and C. Sander, “The microRNA.org resource: targets and expression,” Nucleic Acids Research, vol. 36, no. database issue, pp. D149–D153, 2008.
View at: Publisher Site | Google Scholar
D. W. Huang, B. T. Sherman, and R. A. Lempicki, “Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists,” Nucleic Acids Research, vol. 37, no. 1, pp. 1–13, 2009.
View at: Publisher Site | Google Scholar
D. W. Huang, B. T. Sherman, and R. A. Lempicki, “Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources,” Nature Protocols, vol. 4, no. 1, pp. 44–57, 2009.
View at: Publisher Site | Google Scholar
S. G. Jantzen, B. J. G. Sutherland, D. R. Minkley, and B. F. Koop, “GO Trimming: systematically reducing redundancy in large Gene Ontology datasets,” BMC Research Notes, vol. 4, no. 1, 2011.
View at: Publisher Site | Google Scholar
K. J. Livak and T. D. Schmittgen, “Analysis of relative gene expression data using real-time quantitative PCR and the 2^{−ΔΔ _C_}_T method,” Methods, vol. 25, no. 4, pp. 402–408, 2001.
View at: Publisher Site | Google Scholar
M. J. Axtell, J. O. Westholm, and E. C. Lai, “Vive la différence: biogenesis and evolution of microRNAs in plants and animals,” Genome Biology, vol. 12, no. 4, p. 221, 2011.
View at: Publisher Site | Google Scholar
X. Ma, L. Zhou, and S. Zheng, “Transcriptome analysis revealed key prognostic genes and microRNAs in hepatocellular carcinoma,” PeerJ, vol. 8, no. 5, p. e8930, 2020.
View at: Publisher Site | Google Scholar
X. Wang, J. Gao, B. Zhou, J. Xie, G. Zhou, and Y. Chen, “Identification of prognostic markers for hepatocellular carcinoma based on miRNA expression profiles,” Life Sciences, vol. 232, p. 116596, 2019.
View at: Publisher Site | Google Scholar
T. Inamoto, H. Uehara, Y. Akao et al., “A panel of microRNA signature as a tool for predicting survival of patients with urothelial carcinoma of the bladder,” European Urology Supplements, vol. 16, 2018.
View at: Publisher Site | Google Scholar
J. Gao, X. Zhang, and C. Weu, “Involvement of differentially expressed miRNA in prognosis of pancreatic cancers,” Carcinogenesis, Teratogenesis and Mutagenesis, vol. 31, no. 3, pp. 173–179, 2019.
View at: Google Scholar
S. N. Fesharaki, S. N. Fesharaki, A. Esmaeili, M. Azadeh, and K. Ghaedi, “SNP rs1803622 in hsa-miR-548g binding site at _GAPDH_ alters susceptibility to breast cancer,” Gene Reports, vol. 21, p. 100845, 2020.
View at: Publisher Site | Google Scholar
N. Wu, Y. Lu, and J. Z. Liang, “Expression and correlation of survivin and hsa-miR-542-3p in patients with oral squamous cell carcinoma,” Shanghai Journal of Stomatology, vol. 25, no. 6, pp. 720–724, 2016.
View at: Google Scholar
J.-m. Wang, Y. Gu, Y. Zhang et al., “Deep-sequencing identification of differentially expressed miRNAs in decidua and villus of recurrent miscarriage patients,” Archives of Gynecology and Obstetrics, vol. 293, no. 5, pp. 1125–1135, 2016.
View at: Publisher Site | Google Scholar
C. Gungormez, A. H. Gumushan, N. Dilsiz, and E. Borazan, “Novel miRNAs as potential biomarkers in stage II colon cancer: microarray analysis,” Molecular Biology Reports, vol. 46, no. 4, pp. 4175–4183, 2019.
View at: Publisher Site | Google Scholar
J. Thélu, P. Rossio, and B. Favier, “Notch signalling is linked to epidermal cell differentiation level in basal cell carcinoma, psoriasis and wound healing,” BMC Family Practice, vol. 2, 2002.
View at: Publisher Site | Google Scholar
M. Tanaka, T. Terui, S. Sasai, and H. Tagami, “Basal cell carcinoma showing connections with epidermal cysts,” ournal of the European Academy of Dermatology and Venereology, vol. 17, no. 5, pp. 581-582, 2003.
View at: Publisher Site | Google Scholar
F. Festa, R. Kumar, S. Sanyal et al., “Basal cell carcinoma and variants in genes coding for immune response, DNA repair, folate and iron metabolism,” Mutation Research, vol. 574, no. 1-2, pp. 105–111, 2005.
View at: Publisher Site | Google Scholar
S. Rangwala and K. Y. Tsai, “Roles of the immune system in skin cancer,” British Journal of Dermatology, vol. 165, no. 5, pp. 953–965, 2011.
View at: Publisher Site | Google Scholar
Z. R. Li, Y. Jiang, J. Z. Hu, Y. Chen, and Q. Z. Liu, “SOX2 knockdown inhibits the migration and invasion of basal cell carcinoma cells by targeting the SRPK1-mediated PI3K/AKT signaling pathway,” Oncology Letters, vol. 17, no. 2, pp. 1617–1625, 2019.
View at: Publisher Site | Google Scholar
A. G. Kutikhin and A. E. Yuzhalin, “C-type lectin receptors and RIG-I-like receptors: new points on the oncogenomics map,” Cancer Management and Research, vol. 4, pp. 39–53, 2012.
View at: Publisher Site | Google Scholar
J. Zhu, Z. Zheng, J. Wang et al., “Different miRNA expression profiles between human breast cancer tumors and serum,” Frontiers in Genetics, vol. 5, 2014.
View at: Publisher Site | Google Scholar
M. Ma, M. He, Q. Jiang et al., “miR-487a promotes TGF-β1-induced EMT, the migration and invasion of breast cancer cells by directly targeting MAGI2,” International Journal of Biological Sciences, vol. 12, no. 4, pp. 397–408, 2016.
View at: Publisher Site | Google Scholar
B. Liu, W. Zhou, H. Jiang, Z. Xiang, and L. Wang, “miR-1303 promotes the proliferation, migration and invasion of prostate cancer cells through regulating the Wnt/β-catenin pathway by targeting DKK3,” Experimental and Therapeutic Medicine, vol. 18, no. 6, pp. 4747–4757, 2019.
View at: Publisher Site | Google Scholar
Y. Feng, T. Sun, Y. Yu, Y. Gao, X. Wang, and Z. Chen, “MicroRNA-370 inhibits the proliferation, invasion and EMT of gastric cancer cells by directly targeting PAQR4,” Journal of Pharmacological Sciences, vol. 138, no. 2, pp. 96–106, 2018.
View at: Publisher Site | Google Scholar
M. Adamkov, E. Halasova, J. Rajcani et al., “Relation between expression pattern of p53 and survivin in cutaneous basal cell carcinomas,” Medical Science Monitor International Medical Journal of Experimental and Clinical Research, vol. 17, no. 3, pp. BR74–BR80, 2011.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2021 Hai-Peng Wei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

409

Downloads

667

Citations