Rewired functional regulatory networks among miRNA isoforms (isomiRs) from let-7 and miR-10 gene families in cancer

Keywords: ACC, adrenocortical carcinoma, BLCA, bladder urothelial carcinoma, BRCA, breast invasive carcinoma, CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma, CHOL, cholangiocarcinoma, COAD, colon adenocarcinoma, ESCA, esophageal carcinoma, GBM, glioblastoma multiforme, HNSC, head and neck squamous cell carcinoma, KICH, kidney chromophobe, KIRC, kidney renal clear cell carcinoma, KIRP, kidney renal papillary cell carcinoma, LAML, acute myeloid leukemia, LIHC, liver hepatocellular carcinoma, LGG, brain Lower grade glioma, LUAD, lung adenocarcinoma, LUSC, lung squamous cell carcinoma, MESO, mesothelioma, OV, ovarian serous cystadenocarcinoma, PAAD, pancreatic adenocarcinoma, PCPG, pheochromocytoma and paraganglioma, PRAD, prostate adenocarcinoma, READ, rectum adenocarcinoma, SARC, sarcoma, SKCM, skin cutaneous melanoma, STAD, stomach adenocarcinoma, TGCT, testicular germ cell tumors, THCA, thyroid carcinoma, THYM, thymoma, TSG, tumor suppressor gene, UCEC, uterine corpus endometrial carcinoma, UCS, uterine carcinosarcoma, UVM, uveal melanoma, MicroRNA (miRNA), IsomiR, Let-7, miR-10, Network, Function

Authors: Tingming Liang, Leng Han, Li Guo

Affiliations: Jiangsu Key Laboratory for Molecular and Medical Biotechnology, School of Life Science, Nanjing Normal University, Nanjing 210023, China; Department of Biochemistry and Molecular Biology, The University of Texas Health Science Center at Houston McGovern Medical School, Houston, TX 77030, USA; Department of Bioinformatics, Smart Health Big Data Analysis and Location Services Engineering Lab of Jiangsu Province, School of Geographic and Biologic Information, Nanjing University of Posts and Telecommunications, Nanjing, China

License: © 2020 The Author(s) CC BY 4.0 This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

Article links: DOI: 10.1016/j.csbj.2020.05.001  |  PubMed: 32542110  |  PMC: PMC7280754

Relevance: Relevant: mentioned in keywords or abstract

Full text: PDF (2.5 MB)

Introduction

MicroRNA (miRNA, about 22-nt) plays important regulatory roles in multiple biological processes by negatively regulating gene expression at the transcriptional or posttranscriptional levels ref. [1], ref. [2]. These small non-coding RNA molecules are also promising biomarkers for the diagnosis of cancer. Now, an increasing number of studies suggest that there are multiple varieties of miRNAs, termed isomiRs, which have divergent sequences and expression patterns. These diverse miRNA isoforms can be generated by nucleotidyl transferases, RNA editing during the miRNA maturation process ref. [3], which simultaneously leads to sequence variation (especially at the 3ʹ end) and complex small RNA and mRNA interactions. The specificity of miRNAs is largely the result of binding of the miRNA seed region (nucleotides 2–8) to the 3ʹ UTR of the targeted mRNAs ref. [4], ref. [5]. The variation introduced by miRNA isoforms contributes to a dynamic microRNAome and rewiring of the coding-non-coding RNA regulatory network ref. [6]. Therefore, the expression, function, and evolution of isomiRs is closely correlated with their clustered and/or homologous miRNA genes, implicating a potential synergistic effect of the diverse isomiRs in the regulatory network of a miRNA locus.

miRNAs from the same gene family or gene cluster have been studied and shown to have a similar function, such as those from the let-7 gene family ref. [7], ref. [8] and the miR-17-92 gene cluster ref. [9], ref. [10]. Several studies have focused on the multiple isomiRs that act as regulatory molecules ref. [11], ref. [12]. The most discriminatory isomiRs happen to also be differentially expressed between normal tissue and cancer, and they can successfully classify datasets from 32 cancers ref. [13]. Furthermore, some isomiRs exhibit cancer heterogeneity ref. [14] and gender-dependent expression profiles ref. [15]. The abundance profiles of isomiRs and tRNA-derived fragments associate with various molecular phenotypes, metastatic disease, and patient survival in uveal melanoma ref. [16]. It is still unclear whether regulatory networks are disrupted by the varied miRNA isoforms. Based on the potential expression and functional correlations among isomiRs that are similar to their homologous and/or clustered miRNAs, it is important to investigate potential relationships between isomiRs and homologous and/or clustered miRNAs. These relationships will be crucial in exploring miRNA diversity, the biological roles of isomiRs, and the importance of isomiR formation during the miRNA maturation process.

In this study, we selected the let-7 gene family to discuss potential relationships between isomiRs and their homologous and/or clustered miRNAs. Firstly, the let-7 gene family is an ancient family of miRNAs initially discovered as a heterochronic gene in Caenorhabditis elegans ref. [17]. Let-7 miRNAs have been widely and comprehensively studied because of their important biological roles. For example, they can regulate developmental timing in C. elegans, control cardiomyocyte metabolism, and adjust cell size, among other functions ref. [8], ref. [17], ref. [18], ref. [19], ref. [20], ref. [21], ref. [22], ref. [23]. Secondly, let-7 is the largest among all human miRNA gene families, with nine homologous miRNA genes, and has a close relationship with the mir-10 gene family (some members of the two families cluster together, with an inter-miRNA distance of <10,000 bp), which we analyzed simultaneously in this study. Finally, based on the range of experimentally validated and predicted target mRNAs, let-7 miRNAs have a broader mRNA interaction network compared with other gene families, indicating that it is a crucial gene family. Therefore, we selected the let-7 gene family and its related miR-10 gene family to analyze expression and functional networks at the miRNA and isomiR levels, and to comprehensively discuss miRNA diversity based on the isomiRs and miRNAs involved in human cancers.

Materials and methods

Data source

The sequences of miRNAs and pre-miRNAs from the let-7 and miR-10 gene families were obtained from the miRBase database (version 22.1) ref. [24]. The corresponding relationships among miRNAs, including physical relationships that mainly presented as miRNA gene clusters, were determined based on gene distributions (inter-miRNA distance < 10,000 bp).

To understand isomiR expression profiles in the two related gene families, we selected small RNA sequencing data from The Cancer Genome Atlas (TCGA, https://www.cancer.gov/tcga) to perform expression analysis using the R package TCGAbiolinks ref. [25]. The sequence information of isomiR types across different cancers in the study are showed in Table S1. The target mRNAs of abnormally-expressed isomiRs were further queried for association with published cancer hallmark genes (http://software.broadinstitute.org/gsea/msigdb/) ref. [26], Cancer Gene Census (CGC) (http://cancer.sanger.ac.uk/census) ref. [27], core essential genes (essential genes from Hart et al. ref. [28], Blomen et al. ref. [29], and Wang et al. ref. [30]), and tumor suppressor genes (TSG) and oncogenes ref. [31]. The regulation of transcription factors and miRNAs was predicted using TransmiR version 2.0 ref. [32].

Sequence and evolutionary analysis

Firstly, to understand the distribution of genes from the two gene families, particularly that of relevant gene clusters, the relative genomic locations of miRNA genes were estimated using chromPlot ref. [33]. Secondly, relevant homologous miRNAs and their pre-miRNA sequences were aligned using multiple alignment program ClustalX 2.1 ref. [34], while phylogenetic trees were generated based on the sequences of the miRNAs from the two gene families using the neighbor-joining (NJ) method in MEGA 7.0 ref. [35], and the Neighbor-Net method in SplitsTree 4.14 ref. [36]. Phylogenetic networks for the miRNAs and some isomiRs were constructed using Network software via the median-joining (MJ) method. Finally, to understand sequence features among multiple isomiRs and miRNAs from different animal species, sequence logos were constructed using WebLogo version 3.6.0 ref. [37].

Expression analysis of the small RNAs and target mRNAs

To determine the expression profiles of abnormal isomiRs and mRNAs and examine potential interactions between them, we used DESeq2 ref. [38] to obtain deregulated RNAs. Small RNA or mRNAs were classified as abnormal if the |log₂FC| value was >1.5 and the adjusted p-value was <0.05. Further functional analyses were then performed for these abnormal isomiRs. Genes were eliminated from further analyses if they were detected in <50% of all samples. Repeated isomiRs from multicopy pre-miRNAs or homologous miRNAs were also eliminated from further sequence and function analyses.

Prediction and validation of target mRNAs

To understand the potential biological functions of isomiRs from the two gene families, abnormal isomiRs and annotated miRNAs (based on their novel seeds and annotated seeds, respectively) were used to predict their target mRNAs using TargetScan 7.0 ref. [39]. The predicted targets were used in further analyses if the context score values were less than −0.40. The abnormally-expressed target mRNAs were then screened based on their deregulated mRNA expression profiles to identify isomiR:mRNA pairs where the two members showed opposing expression patterns in specific types of cancer (i.e., one was up-regulated while the other was down-regulated). Finally, the selected isomiR:mRNA pairs were further queried for expression correlations by co-expression analysis (n ≥ 50 for both the tumor and normal control groups), and candidate target mRNAs were identified based on a negative Spearman correlation coefficient and a false discovery rate (FDR) value of <0.05.

To further identify novel seeds with significantly different target mRNAs (based on candidate target mRNAs), we used a hypergeometric test to distinguish potential differences between the target mRNAs of novel and annotated seeds:

where N_x,y was the number of mRNAs shared by novel and annotated seeds, N was the total number of involved target mRNAs, and N_x and N_y were the number of targets of the novel and canonical seeds, respectively. Differences were considered significant at p < 0.05, and identified novel seeds were further analyzed to examine potential functional differences compared with their annotated seeds.

Functional, survival, and drug response analyses

To further validate predicted functional differences between the target mRNAs of the novel and annotated seeds, candidate mRNAs were subjected to functional enrichment analysis using DAVID version 6.8 (The Database for Annotation, Visualization and Integrated Discovery) ref. [40], and queried for association with cancer hallmark ref. [26], clinically-actionable genes ref. [41], CGC ref. [27], and core essential genes ref. [28], ref. [29], ref. [30]. Interaction networks of isomiR:mRNA pairs were presented using Cytoscape 3.6.0 ref. [42].

Survival analysis was performed to identify potential differences between isomiRs with length and sequence divergence. A log-rank test was used to assess differences between two groups that were divided based on the median value of a specific isomiR type, with significance accepted at p < 0.05.

Finally, to further understand difference in drug response as a result of the diverse isomiRs, drug response analysis was performed using data available from the Genomics of Drug Sensitivity in Cancer (GDSC) database ref. [43]. Drug sensitivity or resistance was accepted at |DF| > 0.1 and p < 0.05.

Statistical analysis

Hypergeometric, wilcoxon rank-sum test and trend tests were used to estimate significant differences. All statistical analyses were performed using R version 3.4.3.

Results

Overall expression patterns of isomiRs from the let-7 and miR-10 families

Although let-7 and miR-10 miRNAs were always clustered together on chromosomes, they displayed divergent evolutionary patterns and differences in distribution were observed across different animal species. To further understand the expression patterns of the various homologous miRNAs between the two related gene families, we performed a comprehensive expression analysis using public human high-throughput sequencing data (TCGA data).

Small RNAs were analyzed based on a detailed miRNA locus. Abnormally-expressed isomiRs from the two miRNA families in tumor samples were identified, and diverse expression patterns were detected for each miRNA locus and cancer type (Fig. 1A and Supplementary Fig. S1). IsomiR types and seed numbers varied across different loci and across cancer types (Supplementary Fig. S2A). No significant differences in the number of isomiRs were detected between the tumor samples and normal samples using a Wilcoxon rank-sum test (p = 0.4273 for the let-7 family and p = 0.4260 for the miR-10 family), indicating that isomiRs are stably transcribed and processed. The results also showed that isomiRs from the let-7 gene family were expressed relatively uniformly across the cancer types. In comparison, isomiRs from the miR-10 family, namely including miR-99b-5p, miR-125a-5p, miR-10b-5p, miR-10a-5p, showed inconsistent expression patterns across the cancer types (Fig. 1A and Supplementary Fig. S1).

IsomiR types and seed numbers also varied across different loci in breast invasive carcinoma (BRCA) patients (Supplementary Figs. S2B). miRNA let-7b-5p showed the largest number of isomiR types and seed types in the let-7 family, while for the miR-10 family, miR-10a-5p and miR-10b-5p showed an unusually large number of isomiR types and seed numbers compared with other miRNAs. The variations among loci were related to total enrichment levels at the specific locus, while variations among cancer types showed tissue-specificity and identified potential disease-associated small RNAs that may aid in screening and validating disease-associated isomiRs. Based on analysis of individual miRNA genes (pre-miRNA), we found that homologous miRNAs showed divergent expression levels in BRCA patients (n = 1,207) (p < 2.20e−16 for the let-7 family and p = 4.96e−05 for the miR-10 family, Supplementary Fig. S2C). Briefly, let-7b showed the highest level of expression (>10%), while miR-98 had the lowest level (<5%). Overall, miR-10-family miRNAs had higher average expression levels compared with let-7-family miRNAs (Supplementary Fig. S2C). Furthermore, the distribution of expression of homologous miRNAs of the let-7 and miR-10 families revealed that let-7a-3 and let-7b were the most highly expressed miRNAs from the let-7 family in BRCA. miR-10a, miR-10b, and miR-99 were the most highly expressed among the miR-10-family miRNAs (Supplementary Fig. S2D). These results indicated that miRNA enrichment levels are regulated by their precursor molecules.

Overall, reduced let-7 and miR-10 expression levels were frequently observed in 16 cancer types. We also identified a correlation between transcription factors (TF) and miRNAs from the two gene families in BRCA and LUAD patients. The results indicated that over-expression of a TF was always negatively correlated with the expression of its target miRNAs, while decreased expression of a TF was positively correlated with the expression of its target miRNAs (Fig. 1B and C). These findings may explain the reduced expression of let-7 and miR-10 family members in to some extent. In addition, opposing expression patterns in different cancer types (as was observed for let-7f-5p, miR-10a-5p, and miR-10b-5p) is indicative of cancer-specific expression of the small RNAs, while isomiRs with varied sequences but identical seed and functional regions would further complicate cross-talk between small RNAs and mRNAs.

Characterization of isomiR sequences and potential phylogenetic relationships

Typically, only annotated miRNAs are included in analyses, with other isomiRs tending to be ignored. Here, we analyzed the divergence between canonical miRNAs and other isomiRs to identify the potential roles of the ignored isomiRs. As expected, non-annotated isomiRs showed larger distribution intervals compared with their corresponding canonical miRNAs (p = 0.0003, Fig. 2A; p = 0.0088, Fig. S3A, based on isomiRs in Fig. 1A), with the isomiRs showing both sequence and length heterogeneities. Further analysis detected similar length distribution differences between all human isomiRs and all annotated canonical miRNAs (p = 0.8686 for the let-7 family and p = 0.7147 for the miR-10 family, Fig. 2A). Sequence analysis also showed similar sequence features between human isomiRs and annotated miRNAs, especially the highly-conserved miR-5p (Fig. 2B and Supplementary Fig. S3B). These similar distribution patterns were mainly caused by well-conserved features among the small RNAs, including among homologous miRNAs in the same species (such as homologous miRNAs in a gene family) and miRNAs across different animal species.

To determine the phylogenetic relationships among multiple human isomiRs, we attempted to construct a phylogenetic network using the MJ method. Compared with the network for isomiRs from the let-7 gene family, the network for isomiRs from the miR-10 family contained a greater number of hypothetical medians, and the corresponding network of annotated miRNAs also showed large genetic distances between miR-125, miR-99, and miR-100 and miR-10 (Fig. 2C and Supplementary Fig. S3C).

Varied seed regions regulate the diversity of target mRNAs

Our analysis showed that some isomiRs contained novel seed regions that arose via seed shifting events such as alternative and imprecise cleavage during the miRNA maturation process. These changes predominantly involved shifts in the 5ʹ or 3ʹ ends of the seed region, particularly the latter (Fig. 3A and Supplementary Fig. S4A). These novel isomiRs with shifted seeds showed various expression patterns, with some miRNA loci showing significant differences among shifted seeds based on log₂FC values (p < 0.05 based on trend test, Fig. 3A). The two related gene families showed dominant expression at 5p loci, with fewer 3p loci showing abundant enrichment levels. Interestingly, we found that some novel seed sequences were consistent with other miRNAs. For example, the novel seed sequence of let-7-5p was identical to that of annotated miRNA miR-196-5p (Supplementary Fig. S4A).

To understand divergence in the targets of novel seeds versus annotated seeds, we further screened deregulated seeds to perform functional analyses. For the let-7 family, well-conserved homologous genes were only associated with annotated seeds in dominant 5p loci. In comparison, for the miR-10 family, only three annotated seeds in dominant 5p loci from three clusters were identified (Fig. 3B). These novel seeds had several specific target mRNAs that differed from those of their annotated seeds, despite also sharing some common targets. Similar distribution patterns could be detected among different seeds (distribution in let-7, Fig. 3B; distribution in miR-10, Supplementary Fig. S4B), indicating that novel seeds may target novel mRNAs compared with their annotated seeds. Moreover, similar nucleotide compositions in the seed region were identified in all annotated human seeds and all deregulated novel seeds, although these novel seeds were shifted at the 5ʹ or 3ʹ ends (Fig. 3B). This finding indicates that novel seeds have a similar nucleotide distribution, which further ensures their potential interactions with target mRNAs.

Regulatory network rewiring showing function correlations in signaling pathways

Based on the predicted target mRNAs of the seed regions of the isomiRs, we further screened candidate targets based on opposite abnormal expression and correlation analysis. A hypergeometric test was then conducted to distinguish potential differences in target mRNAs between annotated seeds and novel seeds. Our results showed that the targets of the novel seeds from let-7-5p and miR-10-5p were significantly different from those of their annotated seeds (FDR < 0.05, Fig. 4A and B). Pairwise comparisons of these novel seeds revealed gain or loss of target mRNAs, although some seed sequence alterations only involved a single nucleotide difference (Fig. 4B). Gain or loss of target mRNAs among novel seeds, and between novel and annotated seeds, may contribute to further functional variation. For example, enriched GO terms and KEGG pathways were inconsistent among the annotated seeds and novel seeds, with many terms or pathways being gained or lost (Fig. 4C and D). Fewer GO terms and KEGG pathways were shared by annotated seeds and novel seeds, with many being specific to either the annotated or novel seeds. Novel seeds in miR-10-5p were associated with a higher number of genes or pathways associated with biological processes compared with the annotated seed through the acquisition of more specific target mRNAs. In comparison, novel seeds in let-7-5p resulted in decreased molecular function through the loss of some targets compared with their annotated seed (Fig. 4C and D). These functional variations strongly support the hypothesis that the previously ignored novel seeds will contribute to novel biological pathways through the gain and/or loss of target mRNAs.

Further analysis of the abnormal isomiRs showed variations in patient survival as a result of the changes, despite some of the isomiRs having identical seed sequences to those of their annotated miRNAs (Fig. 4E), validating the notion that diverse isomiRs should be treated as independent small RNAs. We found that increases and decreases in the expression of annotated let-7a-5p (location 52–73) had no significant impact (p = 0.8602) on the survival of colon adenocarcinoma (COAD) patients. However, significant differences in survival were associated with abnormal let-7a-5p isomiRs with altered chromosomal locations (location 52–71, p = 0.018; location 53–74, p = 0.0292). On the contrary, a more consistent survival trend was observed for miR-10b-5p compared with the annotated and abnormal isomiRs in kidney renal clear cell carcinoma (KIRC) patients (Fig. 4E).

Rewired regulatory networks of novel seeds compared with annotated seeds

Compared with canonical miRNAs, isomiRs display diverse sequences, length distributions, and expression patterns, along with potential changes in target mRNAs, all of which may lead to further functional variation. Some specific isomiRs have the potential to be used as diagnostic biomarkers of human disease. Among annotated target genes, including actionable genes, CGC, essential genes, oncogenes, and TSG, we identified obvious differences in the targets of novel seeds and annotated seeds (Supplementary Fig. S4C and S4D). Some target mRNAs were associated with specific cancer hallmarks, and the predominant hallmarks sometimes varied among seeds. However, the top three associated hallmarks were consistent between the two families, and included self-sufficiency in growth signals, tissue invasion and metastasis, and insensitivity to antigrowth signals (Fig. 5A). Because of the common targets associated with cancer hallmarks and genes with special annotations, we also generated interaction networks to show the detailed relationships among targets of different seeds based on cancer hallmarks and KEGG pathways, respectively (Fig. 5B and Supplementary Fig. S4E–G). Analysis of the specific targets of each seed showed an inconsistent number of genes with special annotations, gene distributions of associated cancer hallmarks, and KEGG pathways. For example, the interaction network for let-7-5p based on the cancer hallmarks revealed almost equal numbers of hallmark genes for the annotated and novel seeds. Among the novel seeds of let-7-5p, we found that let-7-5p + 1 had the greatest number of target hallmark genes (n = 16) (Fig. 5B). A similar finding was obtained from the rewired interaction network based on the KEGG pathway (Supplementary Fig. S4E).

Finally, to explore the potential differences in drug responses for the isomiRs, we performed an analysis using data from the GDSC database. Based on the expression patterns of the multiple isomiRs of let-7-5p, we selected six isomiRs containing two different seeds for the drug response analysis. Our results confirmed that the various isomiRs induced differing responses across the different cancer types (Fig. 5C and Supplementary Fig. S5). Interestingly, isomiRs that shared the same seed sequence showed diverse drug sensitivities despite only containing alterations at the 3ʹ end of the seed sequence. Further, some isomiRs with different seeds were associated with opposite drug responses. For example, isomiRs located at hg38:chr22:46112751–46112773 (DR = 0.1012, p = 0.0305) and hg38:chr22:46112752–46112774 (DR = − 0.1112, p = 0.0277) caused opposing responses to Elesclomol in KIRP patients, while hg38:chr22:46112751–46112774 (DR = 0.2302, p = 0.0163) and hg38:chr22:46112752–46112772 (DR = − 0.1890, p = 0.0313) caused opposite responses to AICAR drug in CHOL patients (Fig. 5C).

Discussion

This study indicates that the production of multiple highly similar isomiRs is a standard outcome of the evolutionary process. We examined the repertoire of small RNAs belonging to the let-7 and miR-10 gene families, which contain multiple isomiRs with length, sequence, and expression heterogeneities, and assessed potential associations among homologous and/or clustered miRNA loci. The various isomiRs generated from a single miRNA locus in a specific species is a miniature of canonical miRNA repertoires across animal species, and the presence of multiple sequences allows the flexibility to select optimal regulatory molecules. Decreased expression of isomiRs from the let-7 and miR-10 gene families was observed in various cancers (Fig. 1A). This is consistent with previous experimental ref. [44]. The expression of let-7 and miR-10 miRNAs was regulated at various stages during biogenesis, and also depended on the cell type. Many factors controlled the dynamic expression. With a focus on BRCA, we searched for TFs that regulated let-7 and miR-10 members via promoter interactions using TransmiR version 2.0. The expression of some TFs, including EGR1, IL-6, EBF1, and FOXO1, was positively correlated with the expression of isomiRs in BRCA patients. However, other highly expressed TFs, such as APOBEC3B, EZH2, E2F1, GATA3, and MYB, were negatively correlated with isomiRs from the two selected gene families. Therefore, these dysregulated TFs likely contribute to the reduced expression of isomiRs.

Numerous oncogenes and signaling pathways were demonstrated to be targets of let-7 and miR-10 miRNAs in the current study. A number of groups have reported that miRNA variants involved in regulating distinctive target genes play an important biological role in controlling miRNA-mediated gene expression ref. [3], ref. [45], ref. [46]. Tan et al. used luciferase assays to confirm predicted differences between miRNA and isomiR pairs in three out of six cases (including PTEN, CDH1, DNMT3B, NCAM2, HMGA2, and BTG1)ref. [47]. Ma et al. showed that the expression of metabolism-related proteins, including PEPCK, G6Pase, FAS, and CPT1A, was significantly suppressed by canonical miR-27b. In contrast, the expression of these proteins was only slightly inhibited by isomiR-27b-1 or isomiR-27b-2 following transfection into AML-12 cells ref. [48]. It is likely that comprehensive analysis would reveal that this is not a rare phenomenon. Thus, variations in the expression of multiple isomiRs and cross-talk among the isomiRs provide hidden information for the coding-non-coding RNA interaction networks. Our study confirms that isomiRs with novel seeds disturb coding-non-coding RNA regulatory networks and pathways and lead to functional alterations through the gain and/or loss of target mRNAs. Altered interactions and rewired networks strongly support the hypothesis that isomiRs will result in functional variation.

Further, the observed differences in survival and drug response among the different isomiRs confirm that the various isomiRs should not be ignored. The isomiRs of let-7a-5p may have different clinical relevance in terms of survival despite having the same functional regions (Fig. 4E). The coding-non-coding RNA regulatory network based on cancer hallmark genes (Fig. 5B) indicate a novel and more complex network than the single network considered the canonical miRNA only. This is because novel seeds are involved in the coding-non-coding network and altered the original interactions. The rewiring network enrich and complicate the coding-non-coding RNA regulatory network at the isomiR level by adding more seeds and additional target mRNAs. The isomiRs also have the differences in analysis of drug responses (Fig. 5C and Fig. S5), and isomiRs with novel seeds may cause opposite results despite being generated from the same miRNA locus and showing highly similar sequence relationships. Similar results could be extended to multiple isomiRs from homologous miRNAs and canonical miRNAs, although homologous miRNAs are usually studied together because they often have the same functional regions. IsomiRs containing sequence variations may have additional biological functions that should be the focus of further studies. These observations support the hypothesis that isomiRs exhibiting sequence divergence are functional regulatory molecules, and that they may contribute to biological processes via coordinated interactions in regulatory networks. Changes at the miRNA isoform level likely occur in tandem with disease progression, meaning that further attention should be given to cancer-related miRNAs not only at the classical miRNA expression level, but also at the isoform level, when considering disease-ameliorating therapies.

In summary, this study identifies a series of isomiRs of the let-7 and miR-10 gene families that display significant differences from their annotated miRNAs. Our findings also confirm that the isomiR repertoire is closely associated with homologous and/or clustered miRNA genes. IsomiRs with different sequences may have diverse clinical relevance, indicating their potential diagnostic role. Further studies should examine nucleotide variation in small RNAs, and focus on the biological functions of the isomiRs using experimental validation.

Funding

This work was supported by the 10.13039/501100001809National Natural Science Foundation of China [grant number 61771251]; the Key Project of Social Development in Jiangsu Province [grant number BE2016773]; the National Natural Science Foundation of Jiangsu [grant number BK20171443]; the Qinglan Project in Jiangsu Province; Nanjing Normal University Overseas Studies [grant number NJNU-2017]; NUPTSF [grant numbers NY220041]; and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).

Authorship contribution statement

Tingming Liang: Conceptualization, Visualization, Methodology, Software, Writing – original draft. Leng Han: Data curation, Software, Supervision. Li Guo: Conceptualization, Visualization, Investigation, Supervision, Writing – original draft.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

1. H. Guo, N.T. Ingolia, J.S. Weissman, D.P. Bartel. Mammalian microRNAs predominantly act to decrease target mRNA levels. Nature, 2010. [PubMed]
2. J. Krol, I. Loedige, W. Filipowicz. The widespread regulation of microRNA biogenesis, function and decay. Nat Rev Genet, 2010. [PubMed]
3. S.L. Ameres, P.D. Zamore. Diversifying microRNA sequence and function. Nat Rev Mol Cell Biol, 2013. [PubMed]
4. G.C. Shukla, J. Singh, S. Barik. MicroRNAs: processing, maturation, target recognition and regulatory functions. Mol Cell Pharmacol, 2011. [PubMed]
5. D.P. Bartel. MicroRNAs: target recognition and regulatory functions. Cell, 2009. [PubMed]
6. C.T. Neilsen, G.J. Goodall, C.P. Bracken. IsomiRs–the overlooked repertoire in the dynamic microRNAome. Trends Genet, 2012. [PubMed]
7. B. Brueckner, C. Stresemann, R. Kuner, C. Mund, T. Musch. The human let-7a-3 locus contains an epigenetically regulated microRNA gene with oncogenic function. Cancer Res, 2007. [PubMed]
8. K.T. Kuppusamy, D.C. Jones, H. Sperber, A. Madan, K.A. Fischer. Let-7 family of microRNA is required for maturation and adult-like metabolism in stem cell-derived cardiomyocytes. Proc Natl Acad Sci USA, 2015. [PubMed]
9. W. Tan, Y. Li, S.G. Lim, T.M. Tan. miR-106b-25/miR-17-92 clusters: polycistrons with oncogenic roles in hepatocellular carcinoma. World J Gastroenterol, 2014. [PubMed]
10. J. Chen, Z.P. Huang, H.Y. Seok, J. Ding, M. Kataoka. mir-17-92 cluster is required for and sufficient to induce cardiomyocyte proliferation in postnatal and adult hearts. Circ Res, 2013. [PubMed]
11. O. Salem, N. Erdem, J. Jung, E. Münstermann, A. Wörner. The highly expressed 5’isomiR of hsa-miR-140-3p contributes to the tumor-suppressive effects of miR-140 by reducing breast cancer proliferation and migration. BMC Genomics, 2016. [PubMed]
12. A. Bhardwaj, H. Singh, C.M. Trinidad, C.T. Albarracin, K.K. Hunt. The isomiR-140-3p-regulated mevalonic acid pathway as a potential target for prevention of triple negative breast cancer. Breast Cancer Res, 2018. [PubMed]
13. A.G. Telonis, R. Magee, P. Loher, I. Chervoneva, E. Londin. Knowledge about the presence or absence of miRNA isoforms (isomiRs) can successfully discriminate amongst 32 TCGA cancer types. Nucleic Acids Res, 2017. [PubMed]
14. A.G. Telonis, P. Loher, Y. Jing, E. Londin, I. Rigoutsos. Beyond the one-locus-one-miRNA paradigm: microRNA isoforms enable deeper insights into breast cancer heterogeneity. Nucleic Acids Res, 2015. [PubMed]
15. P. Loher, E.R. Londin, I. Rigoutsos. IsomiR expression profiles in human lymphoblastoid cell lines exhibit population and gender dependencies. Oncotarget, 2014. [PubMed]
16. E. Londin, R. Magee, C.L. Shields, S.E. Lally, T. Sato. IsomiRs and tRNA-derived fragments are associated with metastasis and patient survival in uveal melanoma. Pigm Cell Melanoma R, 2019
17. B.J. Reinhart, F.J. Slack, M. Basson, A.E. Pasquinelli, J.C. Bettinger. The 21-nucleotide let-7 RNA regulates developmental timing in Caenorhabditis elegans. Nature, 2000. [PubMed]
18. X. Jiang, J.S. Hawkins, J. Lee, C.O. Lizama, F.L. Bos. Let-7 microRNA-dependent control of leukotriene signaling regulates the transition of hematopoietic niche in mice. Nat Commun, 2017. [PubMed]
19. F. Yu, H. Yao, P. Zhu, X. Zhang, Q. Pan. let-7 regulates self renewal and tumorigenicity of breast cancer cells. Cell, 2007. [PubMed]
20. J. Takamizawa, H. Konishi, K. Yanagisawa, S. Tomida, H. Osada. Reduced expression of the let-7 microRNAs in human lung cancers in association with shortened postoperative survival. Cancer Res, 2004. [PubMed]
21. K.A. Worringer, T.A. Rand, Y. Hayashi, S. Sami, K. Takahashi. The let-7/LIN-41 pathway regulates reprogramming to human induced pluripotent stem cells by controlling expression of prodifferentiation genes. Cell Stem Cell, 2014. [PubMed]
22. S. Manier, J.T. Powers, A. Sacco, S.V. Glavey, D. Huynh. The LIN28B/let-7 axis is a novel therapeutic pathway in multiple myeloma. Leukemia, 2017. [PubMed]
23. L. Yan, J. Zhou, Y. Gao, S. Ghazal, L. Lu. Regulation of tumor cell migration and invasion by the H19/let-7 axis is antagonized by metformin-induced DNA methylation. Oncogene, 2015. [PubMed]
24. A. Kozomara, M. Birgaoanu, S. Griffiths-Jones. miRBase: from microRNA sequences to function. Nucleic Acids Res, 2019. [PubMed]
25. A. Colaprico, T.C. Silva, C. Olsen, L. Garofano, C. Cava. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res, 2016
26. A. Subramanian, P. Tamayo, V.K. Mootha, S. Mukherjee, B.L. Ebert. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA, 2005. [PubMed]
27. P.A. Futreal, L. Coin, M. Marshall, T. Down, T. Hubbard. A census of human cancer genes. Nat Rev Cancer, 2004. [PubMed]
28. T. Hart, M. Chandrashekhar, M. Aregger, Z. Steinhart, K.R. Brown. High-resolution CRISPR screens reveal fitness genes and genotype-specific cancer liabilities. Cell, 2015. [PubMed]
29. V.A. Blomen, P. Majek, L.T. Jae, J.W. Bigenzahn, J. Nieuwenhuis. Gene essentiality and synthetic lethality in haploid human cells. Science, 2015. [PubMed]
30. T. Wang, K. Birsoy, N.W. Hughes, K.M. Krupczak, Y. Post. Identification and characterization of essential genes in the human genome. Science, 2015. [PubMed]
31. B. Vogelstein, N. Papadopoulos, V.E. Velculescu, S. Zhou, L.A. Diaz. Cancer genome landscapes. Science, 2013. [PubMed]
32. Z. Tong, Q. Cui, J. Wang, Y. Zhou. TransmiR v2. 0: an updated transcription factor-microRNA regulation database. Nucleic Acids Res, 2019. [PubMed]
33. K.Y. Orostica, R.A. Verdugo. chromPlot: visualization of genomic data in chromosomal context. Bioinformatics, 2016. [PubMed]
34. M.A. Larkin, G. Blackshields, N.P. Brown, R. Chenna, P.A. Mcgettigan. Clustal w and clustal x version 2.0. Bioinformatics, 2007. [PubMed]
35. S. Kumar, G. Stecher, K. Tamura. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol, 2016. [PubMed]
36. D.H. Huson. SplitsTree: analyzing and visualizing evolutionary data. Bioinformatics, 1998. [PubMed]
37. G.E. Crooks, G. Hon, J.M. Chandonia, S.E. Brenner. WebLogo: a sequence logo generator. Genome Res, 2004. [PubMed]
38. M.I. Love, W. Huber, S. Anders. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol, 2014. [PubMed]
39. V. Agarwal, G.W. Bell, J.W. Nam, D.P. Bartel. Predicting effective microRNA target sites in mammalian mRNAs. Elife, 2015
40. W.H. Da, B.T. Sherman, R.A. Lempicki. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nature Protoc, 2008
41. A. Subramanian, P. Tamayo, V.K. Mootha, S. Mukherjee, B.L. Ebert. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA, 2005. [PubMed]
42. P. Shannon, A. Markiel, O. Ozier, N.S. Baliga, J.T. Wang. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res, 2003. [PubMed]
43. W. Yang, J. Soares, P. Greninger, E.J. Edelman, H. Lightfoot. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res, 2013. [PubMed]
44. J. Schultz, P. Lorenz, G. Gross, S. Ibrahim, M. Kunz. MicroRNA let-7b targets important cell cycle molecules in malignant melanoma cells and interferes with anchorage-independent growth. Cell Res, 2008. [PubMed]
45. A. Hinton, S.E. Hunter, I. Afrikanova, G.A. Jones, A.D. Lopez. sRNA-seq analysis of human embryonic stem cells and definitive endoderm reveals differentially expressed microRNAs and novel IsomiRs with distinct targets. Stem Cells, 2014. [PubMed]
46. K.C. Vickers, P. Sethupathy, J. Baran-Gale, A.T. Remaley. Complexity of microRNA function and the role of isomiRs in lipid homeostasis. J Lipid Res, 2013. [PubMed]
47. G.C. Tan, E. Chan, A. Molnar, R. Sarkar, D. Alexieva. 5′ isomiR variation is of functional and evolutionary importance. Nucleic Acids Res, 2014. [PubMed]
48. M. Ma, Z. Yin, H. Zhong, T. Liang, L. Guo. Analysis of the expression, function, and evolution of miR-27 isoforms and their responses in metabolic processes. Genomics, 2019. [PubMed]