Integrative expression network analysis of microRNA and gene isoforms in sacred lotus

Integrative expression network analysis of microRNA and gene isoforms in sacred lotus | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research article Integrative expression network analysis of microRNA and gene isoforms in sacred lotus Yue Zhang, Razgar Seyed Rahmani, Xingyu Yang, Jinming Chen, Tao Shi This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.2.19330/v3 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 25 Jun, 2020 Read the published version in BMC Genomics → Version 3 posted 12 You are reading this latest preprint version Show more versions Abstract Background: Gene expression is complex and regulated by multiple molecular mechanisms, such as miRNA-mediated gene inhibition and alternative-splicing of pre-mRNAs. However, the coordination of interaction between miRNAs with different splicing isoforms, and the change of splicing isoform in response to different cellular environments are largely unexplored in plants. In this study, we analyzed the miRNA and mRNA transcriptome from lotus ( Nelumbo nucifera ), an economically important flowering plant. Results: Through RNA-seq analyses on miRNAs and their target genes (isoforms) among six lotus tissues, expression of most miRNAs seem to be negatively correlated with their targets and tend to be tissue-specific. Further, our results showed that preferential interactions between miRNAs and hub gene isoforms in one coexpression module which is highly correlated with leaf. Intriguingly, for many genes, their corresponding isoforms were assigned to different co-expressed modules, and they exhibited more divergent mRNA structures including presence and absence of miRNA binding sites, suggesting functional divergence for many isoforms is escalated by both structural and expression divergence. Further detailed functional enrichment analysis of miRNA targets revealed that miRNAs are involved in the regulation of lotus growth and development by regulating plant hormone-related pathway genes. Conclusions: Taken together, our comprehensive analyses of miRNA and mRNA transcriptome elucidate the coordination of interaction between miRNAs and different splicing isoforms, and highlight the functional divergence of many transcript isoforms from the same locus in lotus. Epigenetics & Genomics microRNA transcript isoforms co-expression network sacred lotus Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Background The genetic central dogma only illustrates a portion of gene regulation since gene expression regulation is a multi-layer mechanism involving more processes such as alternative splicing of pre-mRNAs, and non-coding RNA regulation. Among non-coding RNAs, microRNAs (miRNAs) are one of the most important groups that can interact with the gene at the RNA level. In plants, microRNAs (miRNAs) are a class of small endogenous single-stranded noncoding RNAs ranging from 18 to 24 nucleotides in length [1]. The primary miRNAs are derived from MIRNA genes transcribed by RNA polymerase II and further processed by dicer-like 1 (DCL1) to yield the precursor-miRNAs (pre-miRNAs) [2, 3]. The pre-miRNAs are later diced into short miRNA duplexes containing one or two mature miRNAs. Given that many miRNAs are tissue or species-specific, much research has been conducted to explore the function of plant miRNAs indicating that the plant miRNAs play key roles in response to plant development, abiotic and biotic stresses through regulating their target genes [4-6]. The silencing or translational repression of genes containing miRNA binding sites is a post-transcriptional mechanism of gene regulation [7]. Several studies have suggested that a substantial amount of the miRNA targets are transcription factors or stress-response factors that are essential for biological processes. Lacking miRNA regulation, plants would face multiple developmental defects in many critical developmental stages [8-10]. High throughput small RNA sequencing is efficient and accurate to elucidate miRNA expression profiles and has been employed in many plant studies to uncover the roles of miRNAs in organ growth and response to the environmental stimuli [11-14]. Through differential expression analyses, studies found many differentially expressed miRNAs that participate in different processes and pathways such as auxin signal transduction during pollination of maize silks [15] and root development in Arabidopsis [16, 17]. RNA alternative splicing (AS) is another important post-transcriptional regulation mechanism, producing diverse transcript isoforms encoded by the same genes [18]. With the widespread application of full-length transcriptome sequencing technology, plenty of isoforms produced by alternative splicing events were identified in plants [19-21]. The structure variation in transcript isoforms can often result in proteins with altered physical characteristics and molecular functions [22]. In some cases, the presence or absence of the miRNA binding site in the isoform determines the possibility of its silencing by a complementary miRNA, allowing some isoforms to escape from being targeted due to lack of the miRNA binding site. This phenomenon of miRNA escaping through mRNA splicing has been identified in cotton and maize, indicating the gene regulation which can be interplayed by both miRNAs and AS [23, 24]. Nowadays, investigations on the regulated network of miRNA-mRNA interactions have been carried out in some model plants, such as Arabidopsis and rice, to identify the key genes related to abiotic stress [25, 26]. These studies focused on the regulation of miRNA on target gene expression, but the influence of miRNAs on the co-expression network of different splicing isoforms calls for further investigation in the plant. Besides, our understanding of expression and functional divergence of isoforms in response to different developmental and growth factors is impeded by the paucity of relevant case studies in plants [19-22]. Lotus or sacred lotus ( Nelumbo nucifera ) is an important aquatic plant with utility in horticulture, landscape, and medicine, which is widely cultivated in Asia. Our previous deep-sequencing of miRNAs in six different tissue samples uncovers the evolution and diversity of miRNAs in lotus [27]. Meanwhile, by combining the full-length transcriptome sequencing and RNA-seq dataset of lotus, we also identified a large amount of AS events showing tissue-specific regulatory manner [28]. However, the interactions between miRNAs and targets at the isoform level, and the impact of miRNAs on target gene and isoform expression profiles are still unclear. In this study, comparative analyses of expression profiles between miRNAs and their target genes (and isoforms) were carried out, aiming to explore the spatial and temporal regulation of miRNAs in lotus. Combining the identified full-length isoforms and small RNA-seq data, we also comprehensively investigated the interactions between miRNAs and their target isoforms by WGCNA (weighted gene co-expression network analysis) to uncover the impact of miRNAs on the expression and function of their target isoforms. Results Identification of microRNAs in the new lotus reference genome To obtain a more comprehensive miRNA profile, we reanalyzed sRNA-seq datasets from six lotus tissues including leaf, petiole, petal, anther, unpollinated carpel and pollinated carpel, based on an updated miRbase database and an improved chromosome-level genome assembly of lotus. A total of 22.2 million filtered reads were mapped to the known miRNAs in miRBase (Table 1). The ratio of filtered high-quality reads mapped to the miRBase is 0.33%, i.e. a total of 50 866 reads aligned to the reference genome (nelumbo.biocloud.net) (Table 1) [29]. After merging with previous lotus miRNAs [27] and removing the redundant (overlapping) hairpin loci, a total of 1103 potential mature miRNA and 104 miRNA-star sequences were identified, and these miRNAs are produced by 1416 pre-miRNAs (hairpin loci) (Fig. 1a)(Additional file 2: Table S1 and S2). The number of detected mature miRNAs is less than pre-miRNAs because many pre-miRNAs from distinct duplicate MIRNA genes can produce identical (short) mature miRNA sequence, which was also reported in other species (http://mirbase.org). Comparing the origin of the pre-miRNAs with transposable elements (TE) region in genome, 623 (43.99%) pre-miRNAs were found to be TE-related, suggesting that a substantial number of the miRNAs originate from TEs [30, 31]. In addition, only 444 (40.25%) of those mature miRNAs were identified as miRNA in the previous analysis [27]. Furthermore, 235 (19.46%) of miRNAs were known sequences in miRBase database and 528 (43.74%) are novel miRNAs identified in this study. Among these currently identified novel miRNAs, 348 (65.9% of novel) are potentially produced by TE-related MIRNA-likes genes. By length, the 24bp miRNAs are the most abundant while 388 (58.43%) of which are TE-related, supporting that the emerging of novel miRNAs from TEs [32, 33](Additional file 1: Fig.S1). Furthermore, we observed that the frequency of each nucleobase (A, U, C and G) in the miRNAs was close to 25% (Additional file 1: Fig.S2). However, we also determined the frequency of the base of the mature miRNAs, the result showed that the 20bp, 21bp, and 22bp miRNAs preferentially start with ‘U’ at the first base (46.96%, 55.37%, and 61.22%, respectively) (Additional file 1: Fig.S3), while 24bp miRNAs preferred ‘A’ (58.5%). Comparing with miRNA’s first nucleotide bias analysis in other species, we found the bias tendency in 21bp, 22bp and 24bp miRNAs is similar to Camellia japonica [34], pomegranate[35]. Expression dynamics of miRNAs and their target genes across different tissues Through differential regulation in different tissues or developmental stages, miRNAs play pivotal roles in diverse biological processes including development [4, 5]. To gain insight into the miRNA expression pattern across different lotus tissues, we first performed hierarchical clustering on the expression data from our identified mature miRNAs (Fig.1a). Interestingly, we found that the majority of miRNAs are preferentially expressed in specific tissues. Only 110 miRNAs are commonly expressed in all tissues; carpel has the most specific miRNAs, followed by anther (Fig.1b). A total of 1,003 differentially expressed miRNAs were identified. We identified differentially expressed miRNA in other tissues relative to pollinated carpel, and the up-regulated miRNAs outnumber the down-regulated miRNA in the pollinated carpel, indicating that there could be intensive activation of miRNAs in carpel after pollination (Fig.1c). The Pearson correlation correlation coefiicients among gene expression profiles generated by the RNA-seq analysis of biological replicates suggested the high reproducibility between replicates (ave r > 0.859, all p-value < 0.0001)(Additional file 1: Fig.S4). To explore the expression pattern of miRNA target genes among different tissues, pairwise comparisons of these six samples were conducted to identify differentially expressed genes (DEGs). A total of 28 701 DEGs were identified by using the edgeR package. The comparison between anther and petiole shows the most DEGs, whereas the comparison between pollinated carpel and unpollinated carpel reveals the least DEGs (Fig.2a). To explore whether differentially expressed miRNAs might escalate the expression difference of their target genes between tissue samples, we calculated the proportion of DEGs in the target genes of those differentially expressed miRNAs (DEMTGs) and compared it to DEGs in the genome background. The comparison between anther and petiole also exhibits the highest percentage 49.26% (740) of DEMTGs, while the comparison in pollinated carpel and unpollinated carpel has the lowest percentage of 5.07% (18) (Fig.2a). The proportion of DEGs in DEMTGs is generally higher than that of DEGs in all genes for most between-tissue comparisons, especially in the comparison between carpel and leaf, between carpel and petiole ( χ 2 test, all pvalue<0.01), except for the comparison between petiole and leaf. This indicates that the differentially expressed miRNAs among tissues might influence the expression of their targeted gene to some extent. To further explore how intensively the expression pattern of target genes was influenced by the miRNA, the expression correlation analyses between target genes and miRNAs across different tissue samples were carried out (Additional file 2: Table S3). In this study, the correlation coefficient ( r ) between miRNA and target gene is divided into six levels: strong negative correlation (-1 to -0.75), intermediate negative correlation (-0.75 to -0.25), weak negative correlation (-0.25 to 0), weak positive correlation (0 to 0.25), intermediate positive correlation (0.25 to 0.75) and strong positive correlation (0.75 to 1). The result showed a substantial bias toward negative correlations such that the negative correlations are about double comparing with positive correlations (Fig.2b). The intermediate negative correlations and weak negative correlations are prevalent, and the strong negative correlations are the least, suggesting that miRNAs still mainly repress their target genes (Fig.2b). We further investigated the expression level of targeted genes in different samples, which revealed that the expression of targeted genes is varied between samples possibly due to the expression difference of miRNAs between samples (Fig.2c). To validate the potential regulation of miRNA targets, we randomly selected 15 miRNA targeted genes to perform real-time qPCR experiments. We carried out correlation analyses between miRNAs expression and RT-PCR result of target genes and compared with corresponding correlation obtained from RNA-seq expression. Among 15 pairs of correlation between miRNA and target genes, 12 pairs (80%) showed the negative correlation based on both results from RT-PCR and RNA-seq, further revealing the complex regulatory relationships between miRNAs and target genes. (Fig.3, Additional file 1: Fig.S5). Differentially expressed miRNA and their target isoforms Taking advantage of transcript isoform analyses from our previous study [28], we further analyzed the miRNA-target isoforms instead of genes. A total of 10345 unique target isoforms were predicted (Additional file 2: Table S4). Most target isoforms (8850, 85.54%) contain only one miRNA target site; a small portion of isoforms (847, 8.18%) contain two miRNA target sites; the rest contain more than two miRNAs target sites (Additional file 1: Fig.S6a). Notably, the isoforms ‘Nn8g40904.1’ and ‘Nn8g40902.1’ can be bound by many miRNAs, with 38 and 31 homologous miRNAs from the family miR169, respectively. We also calculated the number of regulatory miRNA s per target gene, and expectedly the distributions of the number of regulatory miRNAs for miRNA-targeted genes and miRNA-targeted isoforms are similar (Additional file 1: Fig.S6b). Not all miRNA-targeted genes have all their corresponding isoforms being targeted by miRNAs--there are only 1637 target genes having all of their isoforms targeted by the specific miRNAs, such as ‘Nn3g21300’ ( AFB3 ) (Additional file 1: Fig.S7), whereas there are 2449 target genes with only a portion of their isoforms being targeted, such as ‘Nn3g21564’ (Additional file 1: Fig.S7). We further compared the expression level of miRNA-targeted isoforms and non-miRNA-targeted isoforms from the same genes. Interestingly, we found that miRNA targeted isoforms tend to have significantly higher expression level in all investigated tissue samples, suggesting that the isoforms containing miRNA binding sites are under miRNA-mediated expression tuning and buffering likely because of their high expression level representing the functional importance (Additional file 1: Fig.S8). The most miRNA target sites in gene bodies are on coding regions (CDSs) (74.76%), whereas the 5’-UTRs (9.59%) and 3’-UTRs (15.65%) regions have fewer target sites by miRNAs. Given that a substantial number of TE-related miRNAs were found in this study, it is essential to know if they also have a regulatory role in gene expressions. We found that 43.57% of TE-related miRNAs have a target gene while 50.28% of non-TE-related miRNAs have a target gene, suggesting that the TE-related miRNAs also play an important role in regulating genes (Additional file 2: Table.S2, S4). To understand the biological functions of miRNAs, especially those tissue-specific miRNAs, functional annotation based on gene ontology (GO) was used. We found that only 1979 out of 4086 miRNA target genes were annotated by GO categories (Additional file 2: Table S5; Additional file 1: Fig.S9). Among the most significantly enriched GO terms of target genes are “endonuclease activity,” “regulation of transcription, DNA-templated” and “Cul4-RING ubiquitin ligase complex,” indicating that the genes targeted by miRNA can regulate numerous key processes and many belonging to transcription factors [36, 37]. The specific miRNA may regulate specific genes being crucial in the different developmental stages, and therefore GO functional enrichment analysis was conducted for six samples (Additional file 1: Fig.S10). In anther, the most enriched GO terms are related to plant reproductive processes such as “microtubule organizing center,” “auxin-activated signaling pathway” and “endonuclease activity.” In petiole, the miRNA target genes are enriched in “chloroplast stromal thylakoid” and “leaf development.” Both in the pollinated and unpollinated carpel, the most enriched GO terms are the same, i.e. “sepal development,” “regulation of anthocyanin biosynthetic process” and “miRNA binding.” These results collectively revealed that the functions of the miRNA target genes are closely related to the tissue-specification. Functional differentiation of isoforms in the co-expression networks It is often assumed that the tightly connected genes in the co-expression network are likely participating in the same biological process, and therefore it provides a means to identify functional divergence between isoforms. Here, we performed WCGNA at the transcript isoform level. We found that some isoforms are exhibiting dramatic expression differences among different tissues. To explore the potential function of miRNA-targeted isoforms in different tissue, we first performed a hierarchical clustering analysis of total isoforms, and we found that a substantial portion of isoforms showed strong tissue-specificity (Additional file 1: Fig.S11). After filtering out the lowly expressed (FPKM <0.1) and universally expressed (C.V. of FPKMs across six tissue samples <2) isoforms, 56,583 isoforms were retained to construct a co-expression network by using WGCNA. A total of 10 modules were defined as clusters of major tree branches (Fig.4a), with the module size ranging from 766 to 13,309, and isoforms within the same cluster have high correlation coefficients among each other (Additional file 2: Table S6, Fig.4b). We further investigated correlations between the tissues and the 10 co-expression modules. Most modules are significantly (p <0.05) correlated with single tissue, except that the black module is significantly correlated with both pollinated carpel and unpollinated carpel. Basically, isoforms in each module are over-represented in the corresponding tissue, and the 150 candidate hub isoforms for each module were assigned (Additional file 1: Fig.S12). The correlation analysis between the modules revealed that black, cyan, green and pink module, which are significantly correlated with the three floral organs, also have high correlation among each other, proving the accuracy of the module clustering and the homology of differentiated floral organs (Additional file 1: Fig.S13). Because the leaf and petiole are both vegetative tissues, six modules are significantly correlated with leaf or petiole, respectively. To explore the influence of miRNAs on the co-expression network of isoforms, we calculated the content of miRNA-targeted isoforms and the number of hub isoforms in every module(Additional file 1: Fig.S14). Moreover, our further χ 2 test analysis at module level revealed that only the proportion of isoforms in brown modules being targeted by miRNAs (184/2260, 8.14%) is significantly lower than the corresponding proportion of isoforms in hubs (51/150, 34%) (χ 2 test, p <0.01) (Additional file 1: Fig.S14). This suggested that miRNAs preferentially target hub isoforms in the brown module, which is highly correlated with leaves. The isoforms from the same gene are often translated into protein variants with different structures and, hence, performing different functions [22]. To understand the scale of functional differentiation among isoforms from the same gene, we identified isoforms that were assigned to different modules in the co-expression network. Interestingly, among 11,302 genes with multiple isoforms being assigned to modules, 3029 genes have their isoforms being assigned into different modules (GIDDM). Moreover, 464 of these GIDDMs were targeted by miRNAs. This supports that substantial genes with multiple isoforms show functional divergence between isoforms. For example, “Nn5g29774”, annotated as ‘responding to salt stress’, produce a total of 41 isoforms, and 18 of them were clustered into five modules, including 12 in cyan, three in red, one in pink, one in black and one in brown (Additional file 1: Fig.S15). Among these 18 isoforms belonging to different modules, and five of them were regulated by two miRNAs, one by nnu-miR200 and one by miR-1655-3p. If the isoforms of the same gene are functionally divergent, we assume that these different isoforms might likely convert into different genes (duplicates) to play their independent functions during the long-term evolution. To validate this assumption, we searched the closest homologous isoform in rice and Arabidopsis, respectively, for each lotus isoform. After filtering out genes with only one isoform, the gene can be divided into three categories: different isoforms from the same lotus gene with their closest homologs being different genes in rice or Arabidopsis (I); different isoforms from the same lotus gene with their closest homologs being the same isoform from the same gene in rice or Arabidopsis (II); different isoforms from the same lotus gene with their closest homologs being different isoforms from the same gene in rice or Arabidopsis (III). The results show that the number of genes in category II is the most (Additional file 2: Table S7, Additional file 1: Fig.S16). However, interestingly, when only considered the GIDDM, the proportions of isoforms in category I, were largely increased by 10.3% and 9.8%, respectively in rice and Arabidopsis. The result further substantiates that different isoforms, belonging to different co-expression modules, from the same gene tend to evolve into more divergent sequence structures. Meanwhile, this somehow shows that these isoforms were more likely to convert into different duplicate copies during long-term evolution. MiRNA-targeted isoforms in plant hormone signaling pathways To further elucidate the functions of miRNAs and their target isoforms, we focused on the phytohormone pathways enrichment since they are essential in almost all biological processes in the plant. First, KEGG annotation found that a total of 397 miRNA target genes were assigned to 106 pathways. ‘Plant hormone signaling transduction’ was the third most enriched pathways and represented by 20 genes. These 20 genes are in auxin-, cytokinin-, gibberellin-, abscisic acid-, ethylene, brassinosteroid- and jasmonic acid-associated signaling pathways, targeted by 24 miRNAs (Additional file 1: Fig.S17). Among the 20 signaling genes, 16 of them were assigned to different modules in the co-expression network at isoform level and four were not assigned to any module, suggesting that miRNA target genes in the hormone pathways are mostly functionally relevant in different lotus tissues (Additional file 1: Fig.S18). In auxin signaling pathways, the auxin receptor TRANSPORT INHIBITOR RESPONSE1 ( TIR1 ), the auxin-responsive gene auxin/indole-3-acetic acid ( AUX/IAA ), the auxin response factor ( ARF ) the small auxin up RNA ( SAUR ) are targeted by “nnu-miR393b-1s”, “nnu-miR393b-2b”, “NmiRNA#40_469”, “nnu-miR102”, ”nnu-miR156c-1*”. Combining the expression of miRNA and their target isoforms, the result revealed that the altered expression levels of targeted isoforms were not always negatively correlated with their corresponding miRNAs (Fig.5). For example, the expression pattern of most isoforms from TIR1 was almost negatively regulated by miRNAs (Fig.5). However, the two isoforms of TIR1 , “Nn3g21300.4” and “Nn4g26020.5”, were highly expressed in leaf and petiole, same with the expressed pattern of the corresponding miRNAs “nnu-miR393b-1s” and “nnu-miR393b-2b”. In another example, we found that the low expression of “nnu-miR102” in petiole might be associated with high expression of most targeted isoforms “Nn1g04271” ( ARF ) in petiole except for “Nn1g04271.8” (Fig.5). Meanwhile, there are four isoforms, transcribed by Nn1g04271, clustered into red and magenta modules highly correlated with petiole, suggesting the important regulatory relationship of nnu-miR102- ARF in auxin signaling of the petiole. In addition to the auxin signaling genes, similar regulatory relationship, and expression patterns in abscisic acid signaling were also observed (Additional file 1: Fig.S19). For example, the phosphatase 2C ( PP2C ) is targeted by five miRNAs, which are lowly expressed in the unpollinated carpel. Isoforms from one of two miRNA-targeted genes, “Nn6g35319”, are highly expressed in the unpollinated carpel, whose five isoforms were clustered into the black module. Nevertheless, the miRNA “nnu-miR8” and its targeted genes, homologous to serine/threonine-protein kinase CTR1 , have high expression in anther, which, however, do not have a negative expression pattern between the miRNAs and target genes. Although we found many miRNAs and their target isoforms in plant hormone signaling pathways, their regulatory relationships are much more complex. Discussion Gene expression is complicated and regulated by multiple mechanisms, such as non-coding RNAs and transcript splicing. Increasing pieces of evidence indicate that miRNAs play a vital role in plant growth and development by regulating their target genes [4]. Previously, we have identified a great number of miRNAs in six tissues of lotus, and also unveiled the evolutionary patterns for miRNA families with different ages [27]. Further, our other study on lotus transcript isoforms has been facilitated by full-length transcript sequencing using a combination of PacBio and Illumina [28]. The strategy of the combination of miRNAs and mRNA by deep sequencing has been successfully applied in many plant species, such as soybean[38], peanut[39] and cotton[40]. Therefore, in this study, we focused on how the expression of miRNAs in different lotus tissues influences the expression dynamics of their target genes and especially targeted isoforms. Currently, based on an updated miRbase database and our new version of lotus genome assembly, we identified 1,207 unique mature miRNAs in these lotus tissues, which help us to discover some novel miRNAs which were missed ­by the previous study. Most of the novel identified miRNAs in our current study are related to transposable elements. Several studies found that many miRNAs, including Ta miR1123 that function in vernalization in wheat, 10 miRNAs in Arabidopsis and 38 in rice are derived from miniature inverted-repeat transposable elements (MITEs) [30, 41, 42]. As we found a great number of TE-related miRNA that can also have target genes in lotus, we demonstrated that TEs are an important source to give birth to novel functional miRNAs. Meanwhile, combining miRNA sequencing with the corresponding RNA-seq transcriptome profiles, our study revealed the importance of miRNA-mediated regulation in the growth of different lotus tissues. More interestingly, we elucidated the interactions between miRNA and different transcript isoforms. By building a weight isoform co-expression network, we evaluated the impact of miRNA on isoform expression pattern and uncovered the functional divergence of many isoforms originating from the same gene by partitioning into different co-expression modules. These findings will facilitate our current understanding of gene regulation by miRNAs and splicing isoforms in the plant. In this study, most miRNAs are found to be preferentially accumulated in specific lotus tissues, suggesting that most miRNAs play more specialized roles during plant growth and development. In many case studies in plants, the miRNA sequences also appear to be tissue-specific, supporting that the specialized functions for most miRNAs [11, 43]. In our study, the tissue-specific miRNAs are preferentially enriched in the reproductive tissues, such as carpel and anther. Functional enrichment of miRNA-targeted genes showed that the tissue-specific miRNAs were important regulators in anther development, such as the formation of spindle fiber during mitosis. For example, the miR393 was identified to regulate MPS1 , which regulates cell cycle function during anther development in cotton[44], and the homologous miRNA “nnu-miR393a” was also found to be specifically expressed in lotus anther, reflecting its conserved biological function. Besides, the miR172 family members, which regulate AP2, are involved in forming the sepal and petal primordia in Arabidopsis [45]. Several miR172 family members, such as “nnu-miR172b,” were specifically expressed in lotus anther to probably keep it from being developed into stamen petaloid. Our further miRNA differential expression analyses showed that there are more miRNAs induced in the pollinated carpel than in other tissues, indicating that miRNAs might be more active in tissues that undergo dramatic physiological changes, such as pollen-pistil interaction. One novelty of our study is to investigate the relationship between miRNAs and their target isoforms because there are few studies focus on this aspect. Studies in human elucidated the interaction between miRNA and hub gene in a gene network [46, 47], yet not at isoforms levels. However, the intensive interactions between miRNAs and the brown module hub genes from the co-expression network in our current lotus study highlight the vital regulatory roles of miRNAs in the leaf. More interestingly, due to the different regions of isoforms from the same gene caused by alternative splicing or alternative transcriptional start site, we found some miRNAs can only target a portion of isoforms from the same gene while they miss the target for the other isoforms. Another novelty of our study is that we further explored the functional relevance of target genes between tissues and functional divergence of isoforms by building an isoform-based co-expression network. Meanwhile, the reason why we selected WGCNA for co-expression analysis is its diverse functions and maturity of data processing which has been successfully applied in many gene network and functional genomic studies [38, 48, 49]. Our findings revealed the case that some isoforms from the same gene appear to be divergent in the isoform co-expression network and selectively targeted by miRNAs due to their difference in sequence structure, suggesting the coordination of both structural and expression alteration during gene regulation. Our further homology-based analysis uncovered that those isoforms from the same gene, partitioned into different modules, also tend to have their homologs matching different genes in another plant species. We speculated that gene duplication and splicing isoforms might be interchangeable during the long-term evolution of plants. This phenomenon of ‘isoform-duplicate conversion’ is also found during vertebrate evolution including humans [50]. Given that several studies suggest that miRNA has emerged as key regulators of plant hormone response pathways by affecting their metabolism, distribution, and perception[16], we also focused on miRNA-targeted-isoform relationships in these pathways. It was shown that the miR393 regulates the lateral root development in Arabidopsis and leaf morphogenesis in cucumber by down-regulating the receptor TIR1 during auxin signaling [51, 52], and the similar miRNA-target relationship was found in lotus, suggesting the conserved regulatory relationship during plant evolution. We also found that the SAUR gene related to auxin signaling pathway was targeted by miR156 (nnu-miR156c-1*), as this is also found in Arabidopsis [53]. Additionally, the PPC2 , abscisic acid signal pathway genes, were found to be targeted by miR166 and involved in regulating the plant height in Gossypium hirsutum [54]. Our study found that the PPC2 orthologous gene in lotus is targeted by five miRNAs, and the target isoforms were lowly expressed in the petiole, implying that the miR166- PPC2 might participate in the regulation of plant height in lotus. As plant hormones have pivotal roles in plant development, these hormone-related miRNAs found in our study might provide an important genetic basis for future molecular genetic studies and manipulation. Conclusion In the present study, we conducted systematic miRNA and mRNA transcriptome analyses in six lotus tissues and discovered more novel miRNAs. We showed that most miRNAs have tissue-specificity in expression and a negative correlation with the expression of their targeted genes. Meanwhile, the genes regulated by miRNA are involved in multiple biological processes, especially plant hormone response pathways. Our co-expression network at isoform level by using the WGCNA highlights the core regulatory role of miRNAs as the intensive interactions between miRNAs and hub isoforms. We also found that the isoforms from the same gene can be selectively targeted by miRNAs, and we further explored the functional divergence of these isoforms in both structure and expression. Collectively, this study on the interactions between the miRNAs and isoforms, and functional divergence of isoforms can facilitate our current understanding of the complexity of gene regulation in plants. Methods Plant material and RNA extraction The tissues including leaf, petiole, petal and mature anther were harvested from N. nucifera var. ‘China Antique’ in Wuhan Botanical Garden, CAS (114 o 30’E, 30 o 60’N). The unpollinated and pollinated carpel were collected on a blooming day according to our previous study [27]. The samples were immediately frozen in liquid nitrogen and the total RNAs were extracted from each sample using RNAprep pure Plant Kit (TIANGEN). And the RNA was only used for RT-PCR validation experiments in this study. Meanwhile, the corresponding RNA-seq datasets from our previous study were downloaded from the NCBI Short Read Archive (SRA) accession number PRJNA503979 and PRJNA492157 [59]. Small RNA-seq analysis The small RNA-seq dataset was downloaded from the SRA accession number SRX1591010 [27]. The latest plant miRNA sequences were downloaded from miRbase 22.0 for miRNA annotation. For each sample, the high-quality reads without adaptors were filtered to retain sRNAs reads length from 18bp to 29bp. Before miRNA identification, the retained reads were further processed using the miRDeep2 [55]. After quality filtering, the Bowtie was used to map the high-quality reads to the reference genome (nelumbo.biocloud.net) with zero mismatches [29, 56]. The mapped small RNAs were guided by the known precursor miRNA dataset (miRbase 22.0) to identify potential miRNA precursors in the lotus genome, allowing the duplicate (hairpin) loci less than five. The hairpin structures and the aligned small RNAs were processed as described in the miRDeep-P package [57]. The overlapping (redundant) precursors were removed. Subsequently, we identified TE-related miRNAs as those pre-miRNAs overlapping with transposable elements in the lotus genome. To quantify the expression level of miRNA, the read counts were converted into TPM [58]. Differential expression analysis was performed using edgeR. Differentially expressed miRNAs between two samples were defined as those with Benjamín-Hochberg false discovery rate (FDR) less than 0.05 (FDR < 0.05) and log2 of fold_change more than one. Isoform identification and quantification To identify the isoforms, the PacBio full-length transcripts were downloaded and mapped to the new lotus reference genome as previously described [59]. Illumina RNA-Seq from petiole, leaf, petal, anther, unpollinated carpel and pollinated carpel were also downloaded to quantify isoform expression. In brief, the high-quality reads were mapped to the reference by HISTA2 v2.1.0, and then the FPKM (fragments per kilobase per million) value of the genes and isoforms were calculated by StringTie v1.3.5 using a combination of Illumina and full-length transcript based annotations. The differentially expressed genes and isoforms were identified by the edgeR package using the same threshold used for miRNA differential expression analysis. Prediction of miRNA target genes and isoforms Comparing with the animal miRNA, most miRNAs and their target mRNAs in the plant have near-perfect sequence complementarity [60]. The miRNAs target isoforms were predicted by TargetFinder [61], with a strict screening criterion. Briefly, FASTA search was applied to find the potential targets, while a position-dependent mispair penalty system was used to score the target sequence after assessment of penalties for mismatches, bulges, gaps (+1 per position) and G:U pairs (+0.5 per position)[62]. Penalties were doubled if the mismatch, bulge, gap, or G:U pairs occurred at position 2 to 13 relative to the 5’ end of the miRNA. Only one single nucleotide bulge or singe nucleotide gap was allowed while other parameters were set at default. Finally, only predicted targets with scores of three or less were retained. Gene Ontology, KEGG pathway analyses, and ortholog identification The GOseq package based on Wallenius non-central hypergeometric distribution [63] was used to identify the significantly enriched GO terms and KEGG pathways for each gene set, including miRNA target genes and the differentially expressed genes. The KEGG pathway annotation was performed using the target gene as queries in GhostKOALA (https://www.kegg.jp/). The orthologs among Arabidopsis (TAIR11), rice (http://rice.plantbiology.msu.edu/, version 7.0) and lotus were identified using the OrthoMCL with an e-value < 1e-15 and an inflation parameter of 2.0. Co-expression network analysis of transcript isoforms by WCGNA To determine whether isoform expression is correlated with the tissues, the co-expression networks were constructed based on 15 RNA-seq samples by using WGCNA (v1.0) in R[64]. These 15 samples include five leaf samples, two petioles, two mature anthers, two pollinated carpels, and unpollinated carpels (PRJNA503979 and PRJNA492157). To filter out the silent or constantly expressed isoforms, those with an averaged FPKM > 0.1 and coefficient of variation (C.V.) of FPKM > 2 were retained for the subsequent pipeline of WGCNA. According to Topological Overlap Matrix [65], the transcripts were first clustered hierarchically. The transcripts were assigned to nine co-expression modules using a bottom-up algorithm known as the dynamic hybrid cut method, which was named after randomly assigned colors. Modules were identified under the minimum module size of 600, and the threshold of the module similarity was set at 0.5. The kME of isoform, which measures the correlation between isoforms and modules, was calculated. The top 150 isoforms, which were most correlated with each corresponding module, were defined as hub isoforms according to the WGCNA protocol. RT-qPCR analysis of sampled target genes The cDNA libraries were constructed using 1 µg total RNA from five samples and they were diluted to 100 µL before performing RT-PCR. The miRNA targeted gene-specific primers were designed by using Primer Premier (v.5.0) (Additional file 2: Table S8). The qRT-PCR reactions were performed on QuantStudio 6 Flex (Life Technologies, USA) in a final volume of 10 µg containing 1 µg cDNA. The reaction procedure was initiated at 95°C for 10 min, followed by 42 cycles of 95°C for 15 s, 60°C for 30 s, and 72°C for 30 s with fluorescence detection. The miRNA targeted gene expressions were analyzed using the 2 - △△ CT method, with the lotus β-Actin gene used as an internal standard. Declarations Ethics approval and consent to participate Not applicable Consent for pulication Not applicable Availability of data and materials All data that support the findings of this study are included in the manuscript and supplementary information. Competing interests The authors declare that they have no competing interests. Funding This work was supported by grants from the Strategic Priority Research Program CAS (No. XDB31000000), the National Natural Science Foundation of China (No. 31570220, No. 31870208 and No. 31700197), Youth Innovation Promotion Association of Chinese Academy of Sciences (No. 2019335) and Bureau of Landscaping and Forestry of Wuhan Municipality (No.2018-28). The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript. Authors’ Contributions TS and JC conceived and designed the research. YZ and XY participated in the preparation of the reagents and materials. YZ conducted the experiments. TS and YZ analyzed the data. YZ and RSR wrote the manuscript. All the authors read and approved the manuscript. Acknowledgements We thank all individuals who participated in this study. We thank our friend Owen Edward (USA) for proofreading the manuscript. Abbreviations AS: Alternative Splicing DEG: Differentially Expressed Gene DEMTG: Differentially Expressed miRNA Targeted Gene FDR: False Discovery Rate FPKM: Fragments Per Kilobase of transcript per Million fragments mapped GIDDM: Genes have their Isoforms Divided into Different Modules GO: Gene Ontology KEEG: Kyoto Encyclopedia of Gene and Genomes MIRNAs: miRNA genes TPM: Transcripts per million reads WGCNA: Weight Gene Co-expression Network Analysis References Reinhart BJ, Weinstein EG, Rhoades MW, Bartel B, Bartel DP. MicroRNAs in plants. Genes & development. 2002;16(13):1616-1626. Lee Y, Kim M, Han J, Yeom KH, Lee S, Baek SH, Kim VN. MicroRNA genes are transcribed by RNA polymerase II. EMBO journal. 2004;23(20):4051-4060. Cui J, You C, Chen X. The evolution of microRNAs in plants. Current opinion in plant biology . 2017;35:61-67. Li C, Zhang B. MicroRNAs in Control of Plant Development. Journal of cellular physiology . 2016;231(2):303-313. Ferdous J, Hussain SS, Shi BJ. Role of microRNAs in plant drought tolerance. Plant biotechnology journal. 2015;13(3):293-305. Yao S, Yang Z, Yang R, Huang Y, Guo G, Kong X, Lan Y, Zhou T, Wang H, Wang W, et al. Transcriptional Regulation of miR528 by OsSPL9 Orchestrates Antiviral Response in Rice. Molecular plant. 2019;12(8):1114-1122. Zhang S, Yan S, Zhao J, Xiong H, An P, Wang J, Zhang H, Zhang L. Identification of miRNAs and their target genes in Larix olgensis and verified of differential expression miRNAs. BMC plant biology . 2019;19(1):247. Song X, Li Y, Cao X, Qi Y. MicroRNAs and Their Regulatory Roles in Plant-Environment Interactions. Annual review of plant biology . 2019;70:489-525. Wang JW. Regulation of flowering time by the miR156-mediated age pathway. Journal of experimental botany . 2014;65(17):4723-4730. Chen X. A microRNA as a translational repressor of APETALA2 in Arabidopsis flower development. Science . 2004;303(5666):2022-2025. Pokoo R, Ren S, Wang Q, Motes CM, Hernandez TD, Ahmadi S, Monteros MJ, Zheng Y, Sunkar R. Genotype- and tissue-specific miRNA profiles and their targets in three alfalfa ( Medicago sativa L ) genotypes. BMC genomics . 2018;19(10):913. Celik O, Akdas EY. Tissue-specific transcriptional regulation of seven heavy metal stress-responsive miRNAs and their putative targets in nickel indicator castor bean ( R. communis L. ) plants. Ecotoxicology and environmental safety . 2019;170:682-690. Taller D, Balint J, Gyula P, Nagy T, Barta E, Baksa I, Szittya G, Taller J, Havelda Z. Expansion of Capsicum annuum fruit is linked to dynamic tissue-specific differential expression of miRNA and siRNA profiles. PloS one . 2018;13(7):e0200207. Chen Q, Deng B, Gao J, Zhao Z, Chen Z, Song S, Wang L, Zhao L, Xu W, Zhang C, et al. Comparative Analysis of miRNA Abundance Revealed the Function of Vvi-miR828 in Fruit Coloring in Root Restriction Cultivation Grapevine ( Vitis vinifera L. ). International journal of molecular sciences . 2019;20(16). Li XM, Sang YL, Zhao XY, Zhang XS. High-throughput sequencing of small RNAs from pollen and silk and characterization of miRNAs as candidate factors involved in pollen-silk interactions in maize. PloS one . 2013;8(8):e72852. Curaba J, Singh MB, Bhalla PL. miRNAs in the crosstalk between phytohormone signalling pathways. Journal of experimental botany . 2014;65(6):1425-1438. Jin D, Wang Y, Zhao Y, Chen M. MicroRNAs and their cross-talks in plant development. Journal of genetics and genomics . 2013;40(4):161-170. Filichkin S, Priest HD, Megraw M, Mockler TC. Alternative splicing in plants: directing traffic at the crossroads of adaptation and environmental stress. Current opinion in plant biology . 2015;24:125-135. Qiao D, Yang C, Chen J, Guo Y, Li Y, Niu S, Cao K, Chen Z. Comprehensive identification of the full-length transcripts and alternative splicing related to the secondary metabolism pathways in the tea plant ( Camellia sinensis ). Scientific reports . 2019;9(1):2709. Wang L, Jiang X, Wang L, Wang W, Fu C, Yan X, Geng X. A survey of transcriptome complexity using PacBio single-molecule real-time analysis combined with Illumina RNA sequencing for a better understanding of ricinoleic acid biosynthesis in Ricinus communis. BMC genomics . 2019;20(1):456. Wang J, Deng Y, Zhou Y, Liu D, Yu H, Zhou Y, Lv J, Ou L, Li X, Ma Y, et al. Full-length mRNA sequencing and gene expression profiling reveal broad involvement of natural antisense transcript gene pairs in pepper development and response to stresses. Plant journal. 2019;99(4):763-783. Satyawan D, Kim MY, Lee SH. Stochastic alternative splicing is prevalent in mungbean ( Vigna radiata ). Plant biotechnology journal. 2017;15(2):174-182. Wang M, Wang P, Liang F, Ye Z, Li J, Shen C, Pei L, Wang F, Hu J, Tu L, et al. A global survey of alternative splicing in allopolyploid cotton: landscape, complexity and regulation. New phytologist . 2018;217(1):163-178. Thatcher SR, Danilevskaya ON, Meng X, Beatty M, Zastrow-Hayes G, Harris C, Van Allen B, Habben J, Li B. Genome-Wide Analysis of Alternative Splicing during Development and Drought Stress in Maize. Plant physiology . 2016;170(1):586-599. Chen J, Li L. Multiple Regression Analysis Reveals MicroRNA Regulatory Networks in Oryza sativa under Drought Stress. International journal of genomics . 2018; 9395261. Wu C, Li X, Guo S, Wong SM. Analyses of RNA-Seq and sRNA-Seq data reveal a complex network of anti-viral defense in TCV-infected Arabidopsis thaliana. Scientific reports . 2016;6:36007. Shi T, Wang K, Yang P. The evolution of plant microRNAs: insights from a basal eudicot sacred lotus. Plant journal. 2017;89(3):442-457. Zhang Y, Nyong AT, Shi T, Yang P. The complexity of alternative splicing and landscape of tissue-specific expression in lotus ( Nelumbo nucifera ) unveiled by Illumina- and single-molecule real-time-based RNA-sequencing. DNA research. 2019;26(4):301-311. Shi T, Rahmani RS, Gugger PF, Wang M, Li H, Zhang Y, Li Z, Wang Q, Van de Peer Y, Marchal K et al. Distinct expression and methylation patterns for genes with different fates following a single whole-genome duplication in flowering plant. Molecular biology and evolution 2020. Li Y, Li C, Xia J, Jin Y. Domestication of transposable elements into MicroRNA genes in plants. PloS one . 2011;6(5):e19212. Lu S. De novo origination of MIRNAs through generation of short inverted repeats in target genes. RNA biology . 2019;16(6):846-859. Barrera-Figueroa BE, Gao L, Wu Z, Zhou X, Zhu J, Jin H, Liu R, Zhu JK. High throughput sequencing reveals novel and abiotic stress-regulated microRNAs in the inflorescences of rice. BMC plant biology . 2012;12:132. Roberts JT, Cardin SE, Borchert GM. Burgeoning evidence indicates that microRNAs were initially formed from transposable element sequences. Mobile genetic elements . 2014;4:e29255. Yin H, Fan Z, Li X, Wang J, Liu W, Wu B, Ying Z, Liu L, Liu Z, Li J. Phylogenetic tree-informed microRNAome analysis uncovers conserved and lineage-specific miRNAs in Camellia during floral organ development. Journal of experimental botany . 2016;67(9):2641-2653. Saminathan T, Bodunrin A, Singh NV, Devarajan R, Nimmakayala P, Jeff M, Aradhya M, Reddy UK. Genome-wide identification of microRNAs in pomegranate ( Punica granatum L. ) by high-throughput sequencing. BMC plant biology . 2016;16(1):122. Samad AFA, Sajad M, Nazaruddin N, Fauzi IA, Murad AMA, Zainal Z, Ismail I. MicroRNA and Transcription Factor: Key Players in Plant Regulatory Network. Frontiers in plant science . 2017;8:565. Koyama T, Sato F, Ohme-Takagi M. Roles of miR319 and TCP Transcription Factors in Leaf Development. Plant physiology . 2017; 175(2):874-885. Yu JY, Zhang ZG, Huang SY, Han X, Wang XY, Pan WJ, Qin HT, Qi HD, Yin ZG, Qu KX, et al. Analysis of miRNAs Targeted Storage Regulatory Genes during Soybean Seed Development Based on Transcriptome Sequencing. Genes . 2019;10(6). Chen H, Yang Q, Chen K, Zhao S, Zhang C, Pan R, Cai T, Deng Y, Wang X, Chen Y, et al. Integrated microRNA and transcriptome profiling reveals a miRNA-mediated regulatory network of embryo abortion under calcium deficiency in peanut ( Arachis hypogaea L. ). BMC genomics . 2019;20(1):392. Xie F, Zhang B. microRNA evolution and expression analysis in polyploidized cotton genome. Plant biotechnology journal. 2015;13(3):421-434. Piriyapongsa J, Jordan IK. Dual coding of siRNAs and miRNAs by plant transposable elements. RNA . 2008;14(5):814-821. Yu M, Carver BF, Yan L. TamiR1123 originated from a family of miniature inverted-repeat transposable elements (MITE) including one inserted in the Vrn-A1a promoter in wheat. Plant science . 2014;215-216:117-123. Wang R, Yang Z, Fei Y, Feng J, Zhu H, Huang F, Zhang H, Huang J. Construction and analysis of degradome-dependent microRNA regulatory networks in soybean. BMC genomics . 2019;20(1):534. Chen J, Su P, Chen P, Li Q, Yuan X, Liu Z. Insights into the cotton anther development through association analysis of transcriptomic and small RNA sequencing. BMC plant biology . 2018;18(1):154. Luo Y, Guo Z, Li L. Evolutionary conservation of microRNA regulatory programs in plant flower development. Developmental biology . 2013;380(2):133-144. Ju Q, Zhao YJ, Dong Y, Cheng C, Zhang S, Yang Y, Li P, Ge D, Sun B. Identification of a miRNA-mRNA network associated with lymph node metastasis in colorectal cancer. Oncology letters . 2019;18(2):1179-1188. Zhou XG, Huang XL, Liang SY, Tang SM, Wu SK, Huang TT, Mo ZN, Wang QY. Identifying miRNA and gene modules of colon cancer associated with pathological stage by weighted gene co-expression network analysis. OncoTargets and therapy. 2018;11:2815-2830. Tai Y, Liu C, Yu S, Yang H, Sun J, Guo C, Huang B, Liu Z, Yuan Y, Xia E, et al. Gene co-expression network analysis reveals coordinated regulation of three characteristic secondary biosynthetic pathways in tea plant ( Camellia sinensis ). BMC genomics . 2018;19(1):616. Yu Y, Shi J, Li X, Liu J, Geng Q, Shi H, Ke Y, Sun Q. Transcriptome analysis reveals the molecular mechanisms of the defense response to gray leaf spot disease in maize. BMC genomics . 2018;19(1):742. Abascal F, Tress ML, Valencia A. The evolutionary fate of alternatively spliced homologous exons after gene duplication. Genome biology and evolution. 2015;7(6):1392-1403. Chen H, Li Z, Xiong L. A plant microRNA regulates the adaptation of roots to drought stress. FEBS letters . 2012;586(12):1742-1747. Xu J, Li J, Cui L, Zhang T, Wu Z, Zhu PY, Meng YJ, Zhang KJ, Yu XQ, Lou QF, et al. New insights into the roles of cucumber TIR1 homologs and miR393 in regulating fruit/seed set development and leaf morphogenesis. BMC plant biology . 2017;17(1):130. Marin E, Jouannet V, Herz A, Lokerse AS, Weijers D, Vaucheret H, Nussaume L, Crespi MD, Maizel A. miR390, Arabidopsis TAS3 tasiRNAs, and their AUXIN RESPONSE FACTOR targets define an autoregulatory network quantitatively regulating lateral root growth. Plant cell . 2010;22(4):1104-1117. An W, Gong W, He S, Pan Z, Sun J, Du X. MicroRNA and mRNA expression profiling analysis revealed the regulation of plant height in Gossypium hirsutum. BMC genomics . 2015;16:886. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic acids research . 2012;40(1):37-52. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome biology . 2009;10(3):R25. Yang X, Li L. miRDeep-P: a computational tool for analyzing the microRNA transcriptome in plants. Bioinformatics . 2011;27(18):2614-2615. Fahlgren N, Howell MD, Kasschau KD, Chapman EJ, Sullivan CM, Cumbie JS, Givan SA, Law TF, Grant SR, Dangl JL, et al. High-throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of MIRNA genes. PloS one . 2007;2(2):e219. Ahsan MU, Hayward A, Irihimovitch V, Fletcher S, Tanurdzic M, Pocock A, Beveridge CA, Mitter N. Juvenility and Vegetative Phase Transition in Tropical/Subtropical Tree Crops. Frontiers in plant science . 2019;10:729. Parizotto EA, Dunoyer P, Rahm N, Himber C, Voinnet O. In vivo investigation of the transcription, processing, endonucleolytic activity, and functional relevance of the spatial distribution of a plant miRNA. Genes & development . 2004;18(18):2237-2242. Bo X, Wang S. TargetFinder: a software for antisense oligonucleotide target site selection based on MAST and secondary structures of target mRNA. Bioinformatics . 2005;21(8):1401-1402. Allen E, Xie Z, Gustafson AM, Carrington JC. microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell . 2005;121(2):207-221. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome biology . 2010;11(2):R14. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC bioinformatics . 2008;9:559. Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Statistical applications in genetics and molecular biology . 2005;4:Article17. Table Table.1 Summary of high-quality reads mapped to miRBase. Sample High-quality reads Reads with at least one alignment in miRBase Reads without alignment in miRBase Anther 6826041 11337 (0.17%) 6814704 (99.83%) Unpollinated carpel 7164301 10947(0.15% ） 7153354(99.85%) Pollinated carpel 8943060 13735 (0.15%) 8929325 (99.85%) Petiole 1472529 6118 (0.42%) 1466411 (99.58%) Leaf 1707006 7964 (0.47%) 1699042 (99.53%) Petal 1518064 5984 (0.39%) 1512080 (99.61%) Additional Files Additional file 1: Figure S1-S18 Additional file 2: Table S1-S8 Supplementary Files SupplementFigure513Final.pdf SupplementaryTables.xls Cite Share Download PDF Status: Published Journal Publication published 25 Jun, 2020 Read the published version in BMC Genomics → Version 3 posted Review # 1 received at journal 08 Jun, 2020 Review # 2 received at journal 01 Jun, 2020 Reviewer # 4 agreed at journal 20 May, 2020 Review # 4 received at journal 20 May, 2020 Review # 3 received at journal 20 May, 2020 Reviewer # 2 agreed at journal 19 May, 2020 Reviewer # 3 agreed at journal 19 May, 2020 Reviewers invited by journal 18 May, 2020 Reviewer # 1 agreed at journal 18 May, 2020 Editor assigned by journal 13 May, 2020 Submission checks completed at journal 12 May, 2020 Editor invited by journal 12 May, 2020 You are reading this latest preprint version Show more versions Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-9943","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research article","associatedPublications":[],"authors":[{"id":574585,"identity":"acc2451c-82e6-4de9-89f3-2ecad21e0acc","order_by":1,"name":"Yue Zhang","email":"","orcid":"","institution":"Wuhan Botanical Garden","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Yue","middleName":"","lastName":"Zhang","suffix":""},{"id":574586,"identity":"587d6c19-75a9-422f-94dc-8998b2fb0044","order_by":2,"name":"Razgar Seyed Rahmani","email":"","orcid":"","institution":"Universiteit Gent","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Razgar","middleName":"Seyed","lastName":"Rahmani","suffix":""},{"id":574587,"identity":"07d88908-5030-4ef0-8f39-7b04870a926a","order_by":3,"name":"Xingyu Yang","email":"","orcid":"","institution":"Wuhan Institute of Landscape Architecture","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Xingyu","middleName":"","lastName":"Yang","suffix":""},{"id":574588,"identity":"13395552-d858-48e9-8d5f-a2982ead9561","order_by":4,"name":"Jinming Chen","email":"","orcid":"","institution":"Wuhan Botanical Garden","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Jinming","middleName":"","lastName":"Chen","suffix":""},{"id":574589,"identity":"8b16d433-5a29-4d1a-b18b-fae0fd7d76b9","order_by":5,"name":"Tao Shi","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA6ElEQVRIiWNgGAWjYDACCcYGIFkjJ8/ewEa0lkagnmPGhj0HiNbCALKGObHhRgKRWuRnN7c/+FHBxtg4843Zwy81dvYM7GcP4NVicOdgY2PPGRlmdukcc2OZY8mJDTx5Cfi1SCQ2NvC2sbExzs7dJi3ZwJzAIMFjgN9hMxIbG/+2MfMw3DwL0lJvT1ALw43ExmbeNmYJhhu82yQ/NhxmbCCkxQCoZbbMmWMGhj35340Zjh1PbOPJIeSw9Acf31TU1M9nP5b28EdNtT0/+xkCDkMGQB8xMBAboRDA+IMk5aNgFIyCUTBSAAAcv0Y1TIJBDAAAAABJRU5ErkJggg==","orcid":"","institution":"Wuhan Botanical Garden","correspondingAuthor":true,"submittingAuthor":false,"prefix":"","firstName":"Tao","middleName":"","lastName":"Shi","suffix":""}],"badges":[],"createdAt":"2019-12-18 15:52:32","currentVersionCode":3,"declarations":"","doi":"10.21203/rs.2.19330/v3","doiUrl":"https://doi.org/10.21203/rs.2.19330/v3","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1186/s12864-020-06853-y","type":"published","date":"2020-06-25T12:00:00+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":1177294,"identity":"ff3013ee-257d-4848-96ae-eec602b3d428","added_by":"auto","created_at":"2020-05-26 16:07:10","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":713461,"visible":true,"origin":"","legend":"Summary of the miRNA expression. (a). A global view of the expression profile of all mature miRNAs in six tissues. (b). The UpSet plot summarizes the presence of mature miRNA in six tissues. The bottom left horizontal bar graph shows the total number of mature miRNA in per tissue. The circles in each panel’s matrix represent the unique and common parts in Venn diagram sections (unique and overlapping mature miRNAs). Connected circles indicate a certain intersection of mature miRNAs between tissues. The top bar graph in each panel summarizes the number of mature miRNAs for each unique or overlapping combination. (c). The bar plot of differentially expressed miRNAs between six samples. The red is the up-regulated miRNA and the blue is the down-regulated miRNA.","description":"","filename":"Figure.1.png","url":"https://assets-eu.researchsquare.com/files/rs-9943/v3/Figure.1.png"},{"id":1177297,"identity":"f9f85cbd-0f1e-44c7-b65a-3891f326918d","added_by":"auto","created_at":"2020-05-26 16:07:11","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":4764873,"visible":true,"origin":"","legend":"Relationship between miRNAs and their target gene in expression and location. (a). Impact of differentially expressed miRNAs (DEM) on expression of their target genes. Green (DEMTGs): the proportion of differentially expressed genes (DEGs) as targets of differentially expressed miRNAs; brown (DEGs): the proportion of DEGs in the genome background. (b) The distribution of the number of miRNA-target pairs according to Pearson’s correlation coefficient of target gene expression and miRNA expression. (c) The CIRCOS plot of distribution of pre-miRNAs and miRNA target genes in chromosome 1-8. Seven circles from the outside to the inside were the distribution of pre-miRNAs, miRNA target genes in anther, miRNA target genes in leaf, miRNA target genes in petal, miRNA target gene in petiole, miRNA target in unpollinated carpel and miRNA target gene in pollinated carpel, respectively.","description":"","filename":"Figure.2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9943/v3/Figure.2.jpg"},{"id":1177298,"identity":"2faed64c-5b78-449b-b343-33bfcb82b880","added_by":"auto","created_at":"2020-05-26 16:07:11","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":2833396,"visible":true,"origin":"","legend":"Expression profile of selected miRNAs and targeted genes. The expression of miRNAs is showed in line graph on the left column. The expression of targeted gene is showed in the histogram graph on the right column. The correlation coefficient between miRNA and targeted gene is showed. Gay: RT-PCR result; red: RNA-seq result. The a-d figure show “miR159-3p” and its corresponding targeted gene. The e-h figure show “N_miR171a_207a” and its corresponding targeted gene. The i-l figure show “nnu-miR293” and its corresponding targeted gene.","description":"","filename":"Figure.3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9943/v3/Figure.3.jpg"},{"id":1177299,"identity":"9e7d978a-8773-4332-9b3b-141a1ddd4b07","added_by":"auto","created_at":"2020-05-26 16:07:11","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":749293,"visible":true,"origin":"","legend":"The co-expression network of filtered isoforms. (a). Hierarchical cluster tree and color bands indicating 9 modules identified by weighted isoforms co-expression network analysis. (b). The analysis of module-trait correlation. Each row represents a module and each column represents a specific sample. Each cell at the row-column intersection is color-coded by correlation according to the color legend.","description":"","filename":"Figure.4.png","url":"https://assets-eu.researchsquare.com/files/rs-9943/v3/Figure.4.png"},{"id":1177300,"identity":"d19a7023-dd18-4a48-a081-5dee2a5bf060","added_by":"auto","created_at":"2020-05-26 16:07:11","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":1070855,"visible":true,"origin":"","legend":"Correlation analysis of expression between miRNAs and auxin-associated ortholog genes in lotus.","description":"","filename":"Figure.5.png","url":"https://assets-eu.researchsquare.com/files/rs-9943/v3/Figure.5.png"},{"id":13532463,"identity":"2c0a52c2-319f-4526-9567-0f82c32927fe","added_by":"auto","created_at":"2021-09-17 01:17:07","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1226225,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-9943/v3/3ab82f2c-b7b3-45eb-bce7-7464ebbeda04.pdf"},{"id":1177296,"identity":"50f27b05-2953-426f-8601-1d146a50c31e","added_by":"auto","created_at":"2020-05-26 16:07:11","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":19899158,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementFigure513Final.pdf","url":"https://assets-eu.researchsquare.com/files/rs-9943/v3/SupplementFigure513Final.pdf"},{"id":1177295,"identity":"c2087b13-56bd-4866-b7b0-24b93146e561","added_by":"auto","created_at":"2020-05-26 16:07:11","extension":"xls","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":12377088,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTables.xls","url":"https://assets-eu.researchsquare.com/files/rs-9943/v3/SupplementaryTables.xls"}],"financialInterests":"","formattedTitle":"\u003cp\u003eIntegrative expression network analysis of microRNA and gene isoforms in sacred lotus\u003c/p\u003e","fulltext":[{"header":"Background","content":"\u003cp\u003eThe genetic central dogma only illustrates a portion of gene regulation since gene expression regulation is a multi-layer mechanism involving more processes such as alternative splicing of pre-mRNAs, and non-coding RNA regulation. Among non-coding RNAs, microRNAs (miRNAs) are one of the most important groups that can interact with the gene at the RNA level. In plants, microRNAs (miRNAs) are a class of small endogenous single-stranded noncoding RNAs ranging from 18 to 24 nucleotides in length [1]. The primary miRNAs are derived from MIRNA genes transcribed by RNA polymerase II and further processed by dicer-like 1 (DCL1) to yield the precursor-miRNAs (pre-miRNAs) [2, 3]. The pre-miRNAs are later diced into short miRNA duplexes containing one or two mature miRNAs. Given that many miRNAs are tissue or species-specific, much research has been conducted to explore the function of plant miRNAs indicating that the plant miRNAs play key roles in response to plant development, abiotic and biotic stresses through regulating their target genes [4-6].\u003c/p\u003e\n\u003cp\u003eThe silencing or translational repression of genes containing miRNA binding sites is a post-transcriptional mechanism of gene regulation [7]. Several studies have suggested that a substantial amount of the miRNA targets are transcription factors or stress-response factors that are essential for biological processes. Lacking miRNA regulation, plants would face multiple developmental defects in many critical developmental stages [8-10]. High throughput small RNA sequencing is efficient and accurate to elucidate miRNA expression profiles and has been employed in many plant studies to uncover the roles of miRNAs in organ growth and response to the environmental stimuli [11-14]. Through differential expression analyses, studies found many differentially expressed miRNAs that participate in different processes and pathways such as auxin signal transduction during pollination of maize silks [15] and root development in Arabidopsis [16, 17].\u003c/p\u003e\n\u003cp\u003eRNA alternative splicing (AS) is another important post-transcriptional regulation mechanism, producing diverse transcript isoforms encoded by the same genes [18]. With the widespread application of full-length transcriptome sequencing technology, plenty of isoforms produced by alternative splicing events were identified in plants [19-21]. The structure variation in transcript isoforms can often result in proteins with altered physical characteristics and molecular functions [22]. In some cases, the presence or absence of the miRNA binding site in the isoform determines the possibility of its silencing by a complementary miRNA, allowing some isoforms to escape from being targeted due to lack of the miRNA binding site. This phenomenon of miRNA escaping through mRNA splicing has been identified in cotton and maize, indicating the gene regulation which can be interplayed by both miRNAs and AS [23, 24]. Nowadays, investigations on the regulated network of miRNA-mRNA interactions have been carried out in some model plants, such as Arabidopsis and rice, to identify the key genes related to abiotic stress [25, 26]. These studies focused on the regulation of miRNA on target gene expression, but the influence of miRNAs on the co-expression network of different splicing isoforms calls for further investigation in the plant. Besides, our understanding of expression and functional divergence of isoforms in response to different developmental and growth factors is impeded by the paucity of relevant case studies in plants [19-22].\u003c/p\u003e\n\u003cp\u003eLotus or sacred lotus (\u003cem\u003eNelumbo nucifera\u003c/em\u003e) is an important aquatic plant with utility in horticulture, landscape, and medicine, which is widely cultivated in Asia. Our previous deep-sequencing of miRNAs in six different tissue samples uncovers the evolution and diversity of miRNAs in lotus [27]. Meanwhile, by combining the full-length transcriptome sequencing and RNA-seq dataset of lotus, we also identified a large amount of AS events showing tissue-specific regulatory manner [28]. However, the interactions between miRNAs and targets at the isoform level, and the impact of miRNAs on target gene and isoform expression profiles are still unclear. In this study, comparative analyses of expression profiles between miRNAs and their target genes (and isoforms) were carried out, aiming to explore the spatial and temporal regulation of miRNAs in lotus. Combining the identified full-length isoforms and small RNA-seq data, we also comprehensively investigated the interactions between miRNAs and their target isoforms by WGCNA (weighted gene co-expression network analysis) to uncover the impact of miRNAs on the expression and function of their target isoforms.\u003c/p\u003e"},{"header":"Results","content":"\u003cp\u003e\u003cstrong\u003eIdentification of microRNAs in the new lotus reference genome\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo obtain a more comprehensive miRNA profile, we reanalyzed sRNA-seq datasets from six lotus tissues including leaf, petiole, petal, anther, unpollinated carpel and pollinated carpel, based on an updated miRbase database and an improved chromosome-level genome assembly of lotus. A total of 22.2 million filtered reads were mapped to the known miRNAs in miRBase (Table 1). The ratio of filtered high-quality reads mapped to the miRBase is 0.33%, i.e. a total of 50 866 reads aligned to the reference genome (nelumbo.biocloud.net) (Table 1) [29]. After merging with previous lotus miRNAs [27] and removing the redundant (overlapping) hairpin loci, a total of 1103 potential mature miRNA and 104 miRNA-star sequences were identified, and these miRNAs are produced by 1416 pre-miRNAs (hairpin loci) (Fig. 1a)(Additional file 2: Table S1 and S2). The number of detected mature miRNAs is less than pre-miRNAs because many pre-miRNAs from distinct duplicate MIRNA genes can produce identical (short) mature miRNA sequence, which was also reported in other species (http://mirbase.org). Comparing the origin of the pre-miRNAs with transposable elements (TE) region in genome, 623 (43.99%) pre-miRNAs were found to be TE-related, suggesting that a substantial number of the miRNAs originate from TEs [30, 31]. In addition, only 444 (40.25%) of those mature miRNAs were identified as miRNA in the previous analysis [27]. Furthermore, 235 (19.46%) of miRNAs were known sequences in miRBase database and 528 (43.74%) are novel miRNAs identified in this study. Among these currently identified novel miRNAs, 348 (65.9% of novel) are potentially produced by TE-related MIRNA-likes genes. By length, the 24bp miRNAs are the most abundant while 388 (58.43%) of which are TE-related, supporting that the emerging of novel miRNAs from TEs [32, 33](Additional file 1: Fig.S1). Furthermore, we observed that the frequency of each nucleobase (A, U, C and G) in the miRNAs was close to 25% (Additional file 1: Fig.S2). However, we also determined the frequency of the base of the mature miRNAs, the result showed that the 20bp, 21bp, and 22bp miRNAs preferentially start with \u0026lsquo;U\u0026rsquo; at the first base (46.96%, 55.37%, and 61.22%, respectively) (Additional file 1: Fig.S3), while 24bp miRNAs preferred \u0026lsquo;A\u0026rsquo; (58.5%). Comparing with miRNA\u0026rsquo;s first nucleotide bias analysis in other species, we found the bias tendency in 21bp, 22bp and 24bp miRNAs is similar to \u003cem\u003eCamellia japonica \u003c/em\u003e[34], pomegranate[35].\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eExpression dynamics of miRNAs and their target genes across different tissues\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThrough differential regulation in different tissues or developmental stages, miRNAs play pivotal roles in diverse biological processes including development [4, 5]. To gain insight into the miRNA expression pattern across different lotus tissues, we first performed hierarchical clustering on the expression data from our identified mature miRNAs (Fig.1a). Interestingly, we found that the majority of miRNAs are preferentially expressed in specific tissues. Only 110 miRNAs are commonly expressed in all tissues; carpel has the most specific miRNAs, followed by anther (Fig.1b). A total of 1,003 differentially expressed miRNAs were identified. We identified differentially expressed miRNA in other tissues relative to pollinated carpel, and the up-regulated miRNAs outnumber the down-regulated miRNA in the pollinated carpel, indicating that there could be intensive activation of miRNAs in carpel after pollination (Fig.1c).\u003c/p\u003e\n\u003cp\u003eThe Pearson correlation correlation coefiicients among gene expression profiles generated by the RNA-seq analysis of biological replicates suggested the high reproducibility between replicates (ave \u003cem\u003er\u003c/em\u003e \u0026gt; 0.859, all p-value \u0026lt; 0.0001)(Additional file 1: Fig.S4). To explore the expression pattern of miRNA target genes among different tissues, pairwise comparisons of these six samples were conducted to identify differentially expressed genes (DEGs). A total of 28 701 DEGs were identified by using the edgeR package. The comparison between anther and petiole shows the most DEGs, whereas the comparison between pollinated carpel and unpollinated carpel reveals the least DEGs (Fig.2a). To explore whether differentially expressed miRNAs might escalate the expression difference of their target genes between tissue samples, we calculated the proportion of DEGs in the target genes of those differentially expressed miRNAs (DEMTGs) and compared it to DEGs in the genome background. The comparison between anther and petiole also exhibits the highest percentage 49.26% (740) of DEMTGs, while the comparison in pollinated carpel and unpollinated carpel has the lowest percentage of 5.07% (18) (Fig.2a). The proportion of DEGs in DEMTGs is generally higher than that of DEGs in all genes for most between-tissue comparisons, especially in the comparison between carpel and leaf, between carpel and petiole ( \u0026chi;\u003csup\u003e2\u003c/sup\u003e test, all pvalue\u0026lt;0.01), except for the comparison between petiole and leaf. This indicates that the differentially expressed miRNAs among tissues might influence the expression of their targeted gene to some extent.\u003c/p\u003e\n\u003cp\u003eTo further explore how intensively the expression pattern of target genes was influenced by the miRNA, the expression correlation analyses between target genes and miRNAs across different tissue samples were carried out (Additional file 2: Table S3). In this study, the correlation coefficient (\u003cem\u003er\u003c/em\u003e) between miRNA and target gene is divided into six levels: strong negative correlation (-1 to -0.75), intermediate negative correlation (-0.75 to -0.25), weak negative correlation (-0.25 to 0), weak positive correlation (0 to 0.25), intermediate positive correlation (0.25 to 0.75) and strong positive correlation (0.75 to 1). The result showed a substantial bias toward negative correlations such that the negative correlations are about double comparing with positive correlations (Fig.2b). The intermediate negative correlations and weak negative correlations are prevalent, and the strong negative correlations are the least, suggesting that miRNAs still mainly repress their target genes (Fig.2b). We further investigated the expression level of targeted genes in different samples, which revealed that the expression of targeted genes is varied between samples possibly due to the expression difference of miRNAs between samples (Fig.2c). To validate the potential regulation of miRNA targets, we randomly selected 15 miRNA targeted genes to perform real-time qPCR experiments. We carried out correlation analyses between miRNAs expression\u0026nbsp;and RT-PCR result of target genes and compared with corresponding correlation obtained from RNA-seq expression. Among 15 pairs of correlation between miRNA and target genes, 12 pairs (80%) showed the negative correlation based on both results from RT-PCR and RNA-seq, further revealing the complex regulatory relationships between miRNAs and target genes. (Fig.3, Additional file 1: Fig.S5).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eDifferentially expressed miRNA and their target isoforms\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTaking advantage of transcript isoform analyses from our previous study [28], we further analyzed the miRNA-target isoforms instead of genes. A total of 10345 unique target isoforms were predicted (Additional file 2: Table S4). Most target isoforms (8850, 85.54%) contain only one miRNA target site; a small portion of isoforms (847, 8.18%) contain two miRNA target sites; the rest contain more than two miRNAs target sites (Additional file 1: Fig.S6a). Notably, the isoforms \u0026lsquo;Nn8g40904.1\u0026rsquo; and \u0026lsquo;Nn8g40902.1\u0026rsquo; can be bound by many miRNAs, with 38 and 31 homologous miRNAs from the family miR169, respectively. We also calculated the number of regulatory miRNA\u003cu\u003es per target \u003c/u\u003egene, and expectedly the distributions of the number of regulatory miRNAs for miRNA-targeted genes and miRNA-targeted isoforms are similar (Additional file 1: Fig.S6b). Not all miRNA-targeted genes have all their corresponding isoforms being targeted by miRNAs--there are only 1637 target genes having all of their isoforms targeted by the specific miRNAs, such as \u0026lsquo;Nn3g21300\u0026rsquo; (\u003cem\u003eAFB3\u003c/em\u003e) (Additional file 1: Fig.S7), whereas there are 2449 target genes with only a portion of their isoforms being targeted, such as \u0026lsquo;Nn3g21564\u0026rsquo; (Additional file 1: Fig.S7). We further compared the expression level of miRNA-targeted isoforms and non-miRNA-targeted isoforms from the same genes. Interestingly, we found that miRNA targeted isoforms tend to have significantly higher expression level in all investigated tissue samples, suggesting that the isoforms containing miRNA binding sites are under miRNA-mediated expression tuning and buffering likely because of their high expression level representing the functional importance (Additional file 1: Fig.S8). The most miRNA target sites in gene bodies are on coding regions (CDSs) (74.76%), whereas the 5\u0026rsquo;-UTRs (9.59%) and 3\u0026rsquo;-UTRs (15.65%) regions have fewer target sites by miRNAs. Given that a substantial number of TE-related miRNAs were found in this study, it is essential to know if they also have a regulatory role in gene expressions. We found that 43.57% of TE-related miRNAs have a target gene while 50.28% of non-TE-related miRNAs have a target gene, suggesting that the TE-related miRNAs also play an important role in regulating genes (Additional file 2: Table.S2, S4).\u003c/p\u003e\n\u003cp\u003eTo understand the biological functions of miRNAs, especially those tissue-specific miRNAs, functional annotation based on gene ontology (GO) was used. We found that only 1979 out of 4086 miRNA target genes were annotated by GO categories (Additional file 2: Table S5; Additional file 1: Fig.S9). Among the most significantly enriched GO terms of target genes are \u0026ldquo;endonuclease activity,\u0026rdquo; \u0026ldquo;regulation of transcription, DNA-templated\u0026rdquo; and \u0026ldquo;Cul4-RING ubiquitin ligase complex,\u0026rdquo; indicating that the genes targeted by miRNA can regulate numerous key processes and many belonging to transcription factors [36, 37]. The specific miRNA may regulate specific genes being crucial in the different developmental stages, and therefore GO functional enrichment analysis was conducted for six samples (Additional file 1: Fig.S10). In anther, the most enriched GO terms are related to plant reproductive processes such as \u0026ldquo;microtubule organizing center,\u0026rdquo; \u0026ldquo;auxin-activated signaling pathway\u0026rdquo; and \u0026ldquo;endonuclease activity.\u0026rdquo; In petiole, the miRNA target genes are enriched in \u0026ldquo;chloroplast stromal thylakoid\u0026rdquo; and \u0026ldquo;leaf development.\u0026rdquo; Both in the pollinated and unpollinated carpel, the most enriched GO terms are the same, i.e. \u0026ldquo;sepal development,\u0026rdquo; \u0026ldquo;regulation of anthocyanin biosynthetic process\u0026rdquo; and \u0026ldquo;miRNA binding.\u0026rdquo; These results collectively revealed that the functions of the miRNA target genes are closely related to the tissue-specification.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunctional differentiation of isoforms in the co-expression networks\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eIt is often assumed that the tightly connected genes in the co-expression network are likely participating in the same biological process, and therefore it provides a means to identify functional divergence between isoforms. Here, we performed WCGNA at the transcript isoform level. We found that some isoforms are exhibiting dramatic expression differences among different tissues. To explore the potential function of miRNA-targeted isoforms in different tissue, we first performed a hierarchical clustering analysis of total isoforms, and we found that a substantial portion of isoforms showed strong tissue-specificity (Additional file 1: Fig.S11). After filtering out the lowly expressed (FPKM \u0026lt;0.1) and universally expressed (C.V. of FPKMs across six tissue samples \u0026lt;2) isoforms, 56,583 isoforms were retained to construct a co-expression network by using WGCNA. A total of 10 modules were defined as clusters of major tree branches (Fig.4a), with the module size ranging from 766 to 13,309, and isoforms within the same cluster have high correlation coefficients among each other (Additional file 2: Table S6, Fig.4b). We further investigated correlations between the tissues and the 10 co-expression modules. Most modules are significantly (p \u0026lt;0.05) correlated with single tissue, except that the black module is significantly correlated with both pollinated carpel and unpollinated carpel. Basically, isoforms in each module are over-represented in the corresponding tissue, and the 150 candidate hub isoforms for each module were assigned (Additional file 1: Fig.S12). The correlation analysis between the modules revealed that black, cyan, green and pink module, which are significantly correlated with the three floral organs, also have high correlation among each other, proving the accuracy of the module clustering and the homology of differentiated floral organs (Additional file 1: Fig.S13). Because the leaf and petiole are both vegetative tissues, six modules are significantly correlated with leaf or petiole, respectively. To explore the influence of miRNAs on the co-expression network of isoforms, we calculated the content of miRNA-targeted isoforms and the number of hub isoforms in every module(Additional file 1: Fig.S14). Moreover, our further \u0026chi;\u003csup\u003e2\u003c/sup\u003e test analysis at module level revealed that only the proportion of isoforms in brown modules being targeted by miRNAs (184/2260, 8.14%) is significantly lower than the corresponding proportion of isoforms in hubs (51/150, 34%) (\u0026chi;\u003csup\u003e2\u003c/sup\u003e test, p \u0026lt;0.01) (Additional file 1: Fig.S14). This suggested that miRNAs preferentially target hub isoforms in the brown module, which is highly correlated with leaves.\u003c/p\u003e\n\u003cp\u003eThe isoforms from the same gene are often translated into protein variants with different structures and, hence, performing different functions [22]. To understand the scale of functional differentiation among isoforms from the same gene, we identified isoforms that were assigned to different modules in the co-expression network. Interestingly, among 11,302 genes with multiple isoforms being assigned to modules, 3029 genes have their isoforms being assigned into different modules (GIDDM). Moreover, 464 of these GIDDMs were targeted by miRNAs. This supports that substantial genes with multiple isoforms show functional divergence between isoforms. For example, \u0026ldquo;Nn5g29774\u0026rdquo;, annotated as \u0026lsquo;responding to salt stress\u0026rsquo;, produce a total of 41 isoforms, and 18 of them were clustered into five modules, including 12 in cyan, three in red, one in pink, one in black and one in brown (Additional file 1: Fig.S15). Among these 18 isoforms belonging to different modules, and five of them were regulated by two miRNAs, one by nnu-miR200 and one by miR-1655-3p.\u003c/p\u003e\n\u003cp\u003eIf the isoforms of the same gene are functionally divergent, we assume that these different isoforms might likely convert into different genes (duplicates) to play their independent functions during the long-term evolution. To validate this assumption, we searched the closest homologous isoform in rice and Arabidopsis, respectively, for each lotus isoform. After filtering out genes with only one isoform, the gene can be divided into three categories: different isoforms from the same lotus gene with their closest homologs being different genes in rice or Arabidopsis (I); different isoforms from the same lotus gene with their closest homologs being the same isoform from the same gene in rice or Arabidopsis (II); different isoforms from the same lotus gene with their closest homologs being different isoforms from the same gene in rice or Arabidopsis (III). The results show that the number of genes in category II is the most (Additional file 2: Table S7, Additional file 1: Fig.S16). However, interestingly, when only considered the GIDDM, the proportions of isoforms in category I, were largely increased by 10.3% and 9.8%, respectively in rice and Arabidopsis. The result further substantiates that different isoforms, belonging to different co-expression modules, from the same gene tend to evolve into more divergent sequence structures. Meanwhile, this somehow shows that these isoforms were more likely to convert into different duplicate copies during long-term evolution.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eMiRNA-targeted isoforms in plant hormone signaling pathways\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo further elucidate the functions of miRNAs and their target isoforms, we focused on the phytohormone pathways enrichment since they are essential in almost all biological processes in the plant. First, KEGG annotation found that a total of 397 miRNA target genes were assigned to 106 pathways. \u0026lsquo;Plant hormone signaling transduction\u0026rsquo; was the third most enriched pathways and represented by 20 genes. These 20 genes are in auxin-, cytokinin-, gibberellin-, abscisic acid-, ethylene, brassinosteroid- and jasmonic acid-associated signaling pathways, targeted by 24 miRNAs (Additional file 1: Fig.S17). Among the 20 signaling genes, 16 of them were assigned to different modules in the co-expression network at isoform level and four were not assigned to any module, suggesting that miRNA target genes in the hormone pathways are mostly functionally relevant in different lotus tissues (Additional file 1: Fig.S18).\u003c/p\u003e\n\u003cp\u003eIn auxin signaling pathways, the auxin receptor \u003cem\u003eTRANSPORT INHIBITOR RESPONSE1\u003c/em\u003e (\u003cem\u003eTIR1\u003c/em\u003e), the \u003cem\u003eauxin-responsive gene auxin/indole-3-acetic acid\u003c/em\u003e (\u003cem\u003eAUX/IAA\u003c/em\u003e), the \u003cem\u003eauxin response factor\u003c/em\u003e (\u003cem\u003eARF\u003c/em\u003e) the \u003cem\u003esmall auxin up RNA\u003c/em\u003e (\u003cem\u003eSAUR\u003c/em\u003e) are targeted by \u0026ldquo;nnu-miR393b-1s\u0026rdquo;, \u0026ldquo;nnu-miR393b-2b\u0026rdquo;, \u0026ldquo;NmiRNA#40_469\u0026rdquo;, \u0026ldquo;nnu-miR102\u0026rdquo;, \u0026rdquo;nnu-miR156c-1*\u0026rdquo;. Combining the expression of miRNA and their target isoforms, the result revealed that the altered expression levels of targeted isoforms were not always negatively correlated with their corresponding miRNAs (Fig.5). For example, the expression pattern of most isoforms from \u003cem\u003eTIR1\u003c/em\u003e was almost negatively regulated by miRNAs (Fig.5). However, the two isoforms of \u003cem\u003eTIR1\u003c/em\u003e, \u0026ldquo;Nn3g21300.4\u0026rdquo; and \u0026ldquo;Nn4g26020.5\u0026rdquo;, were highly expressed in leaf and petiole, same with the expressed pattern of the corresponding miRNAs \u0026ldquo;nnu-miR393b-1s\u0026rdquo; and \u0026ldquo;nnu-miR393b-2b\u0026rdquo;. In another example, we found that the low expression of \u0026ldquo;nnu-miR102\u0026rdquo; in petiole might be associated with high expression of most targeted isoforms \u0026ldquo;Nn1g04271\u0026rdquo; (\u003cem\u003eARF\u003c/em\u003e) in petiole except for \u0026ldquo;Nn1g04271.8\u0026rdquo; (Fig.5). Meanwhile, there are four isoforms, transcribed by Nn1g04271, clustered into red and magenta modules highly correlated with petiole, suggesting the important regulatory relationship of nnu-miR102-\u003cem\u003eARF\u003c/em\u003e in auxin signaling of the petiole. In addition to the auxin signaling genes, similar regulatory relationship, and expression patterns in abscisic acid signaling were also observed (Additional file 1: Fig.S19). For example, the \u003cem\u003ephosphatase 2C\u003c/em\u003e (\u003cem\u003ePP2C\u003c/em\u003e) is targeted by five miRNAs, which are lowly expressed in the unpollinated carpel. Isoforms from one of two miRNA-targeted genes, \u0026ldquo;Nn6g35319\u0026rdquo;, are highly expressed in the unpollinated carpel, whose five isoforms were clustered into the black module. Nevertheless, the miRNA \u0026ldquo;nnu-miR8\u0026rdquo; and its targeted genes, homologous to serine/threonine-protein kinase \u003cem\u003eCTR1\u003c/em\u003e, have high expression in anther, which, however, do not have a negative expression pattern between the miRNAs and target genes. Although we found many miRNAs and their target isoforms in plant hormone signaling pathways, their regulatory relationships are much more complex.\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eGene expression is complicated and regulated by multiple mechanisms, such as non-coding RNAs and transcript splicing. Increasing pieces of evidence indicate that miRNAs play a vital role in plant growth and development by regulating their target genes [4]. Previously, we have identified a great number of miRNAs in six tissues of lotus, and also unveiled the evolutionary patterns for miRNA families with different ages [27]. Further, our other study on lotus transcript isoforms has been facilitated by full-length transcript sequencing using a combination of PacBio and Illumina [28]. The strategy of the combination of miRNAs and mRNA by deep sequencing has been successfully applied in many plant species, such as soybean[38], peanut[39] and cotton[40]. Therefore, in this study, we focused on how the expression of miRNAs in different lotus tissues influences the expression dynamics of their target genes and especially targeted isoforms. Currently, based on an updated miRbase database and our new version of lotus genome assembly, we identified 1,207 unique mature miRNAs in these lotus tissues, which help us to discover some novel miRNAs which were missed \u0026shy;by the previous study. Most of the novel identified miRNAs in our current study are related to transposable elements. Several studies found that many miRNAs, including \u003cem\u003eTa\u003c/em\u003emiR1123 that function in vernalization in wheat, 10 miRNAs in Arabidopsis and 38 in rice are derived from miniature inverted-repeat transposable elements (MITEs) [30, 41, 42]. As we found a great number of TE-related miRNA that can also have target genes in lotus, we demonstrated that TEs are an important source to give birth to novel functional miRNAs. Meanwhile, combining miRNA sequencing with the corresponding RNA-seq transcriptome profiles, our study revealed the importance of miRNA-mediated regulation in the growth of different lotus tissues. More interestingly, we elucidated the interactions between miRNA and different transcript isoforms. By building a weight isoform co-expression network, we evaluated the impact of miRNA on isoform expression pattern and uncovered the functional divergence of many isoforms originating from the same gene by partitioning into different co-expression modules. These findings will facilitate our current understanding of gene regulation by miRNAs and splicing isoforms in the plant.\u003c/p\u003e\n\u003cp\u003eIn this study, most miRNAs are found to be preferentially accumulated in specific lotus tissues, suggesting that most miRNAs play more specialized roles during plant growth and development. In many case studies in plants, the miRNA sequences also appear to be tissue-specific, supporting that the specialized functions for most miRNAs [11, 43]. In our study, the tissue-specific miRNAs are preferentially enriched in the reproductive tissues, such as carpel and anther. Functional enrichment of miRNA-targeted genes showed that the tissue-specific miRNAs were important regulators in anther development, such as the formation of spindle fiber during mitosis. For example, the miR393 was identified to regulate \u003cem\u003eMPS1\u003c/em\u003e, which regulates cell cycle function during anther development in cotton[44], and the homologous miRNA \u0026ldquo;nnu-miR393a\u0026rdquo; was also found to be specifically expressed in lotus anther, reflecting its conserved biological function. Besides, the miR172 family members, which regulate AP2, are involved in forming the sepal and petal primordia in \u003cem\u003eArabidopsis\u003c/em\u003e[45]. Several miR172 family members, such as \u0026ldquo;nnu-miR172b,\u0026rdquo; were specifically expressed in lotus anther to probably keep it from being developed into stamen petaloid. Our further miRNA differential expression analyses showed that there are more miRNAs induced in the pollinated carpel than in other tissues, indicating that miRNAs might be more active in tissues that undergo dramatic physiological changes, such as pollen-pistil interaction.\u003c/p\u003e\n\u003cp\u003eOne novelty of our study is to investigate the relationship between miRNAs and their target isoforms because there are few studies focus on this aspect. Studies in human elucidated the interaction between miRNA and hub gene in a gene network [46, 47], yet not at isoforms levels. However, the intensive interactions between miRNAs and the brown module hub genes from the co-expression network in our current lotus study highlight the vital regulatory roles of miRNAs in the leaf. More interestingly, due to the different regions of isoforms from the same gene caused by alternative splicing or alternative transcriptional start site, we found some miRNAs can only target a portion of isoforms from the same gene while they miss the target for the other isoforms. Another novelty of our study is that we further explored the functional relevance of target genes between tissues and functional divergence of isoforms by building an isoform-based co-expression network. Meanwhile, the reason why we selected WGCNA for co-expression analysis is its diverse functions and maturity of data processing which has been successfully applied in many gene network and functional genomic studies [38, 48, 49]. Our findings revealed the case that some isoforms from the same gene appear to be divergent in the isoform co-expression network and selectively targeted by miRNAs due to their difference in sequence structure, suggesting the coordination of both structural and expression alteration during gene regulation. Our further homology-based analysis uncovered that those isoforms from the same gene, partitioned into different modules, also tend to have their homologs matching different genes in another plant species. We speculated that gene duplication and splicing isoforms might be interchangeable during the long-term evolution of plants. This phenomenon of \u0026lsquo;isoform-duplicate conversion\u0026rsquo; is also found during vertebrate evolution including humans [50].\u003c/p\u003e\n\u003cp\u003eGiven that several studies suggest that miRNA has emerged as key regulators of plant hormone response pathways by affecting their metabolism, distribution, and perception[16], we also focused on miRNA-targeted-isoform relationships in these pathways. It was shown that the miR393 regulates the lateral root development in \u003cem\u003eArabidopsis\u003c/em\u003e and leaf morphogenesis in cucumber by down-regulating the receptor \u003cem\u003eTIR1\u003c/em\u003e during auxin signaling [51, 52], and the similar miRNA-target relationship was found in lotus, suggesting the conserved regulatory relationship during plant evolution. We also found that the \u003cem\u003eSAUR\u003c/em\u003e gene related to auxin signaling pathway was targeted by miR156 (nnu-miR156c-1*), as this is also found in \u003cem\u003eArabidopsis \u003c/em\u003e[53]. Additionally, the \u003cem\u003ePPC2\u003c/em\u003e, abscisic acid signal pathway genes, were found to be targeted by miR166 and involved in regulating the plant height in \u003cem\u003eGossypium hirsutum\u003c/em\u003e[54]. Our study found that the \u003cem\u003ePPC2\u003c/em\u003e orthologous gene in lotus is targeted by five miRNAs, and the target isoforms were lowly expressed in the petiole, implying that the miR166-\u003cem\u003ePPC2\u003c/em\u003e might participate in the regulation of plant height in lotus. As plant hormones have pivotal roles in plant development, these hormone-related miRNAs found in our study might provide an important genetic basis for future molecular genetic studies and manipulation.\u003c/p\u003e"},{"header":"Conclusion","content":"\u003cp\u003eIn the present study, we conducted systematic miRNA and mRNA transcriptome analyses in six lotus tissues and discovered more novel miRNAs. We showed that most miRNAs have tissue-specificity in expression and a negative correlation with the expression of their targeted genes. Meanwhile, the genes regulated by miRNA are involved in multiple biological processes, especially plant hormone response pathways. Our co-expression network at isoform level by using the WGCNA highlights the core regulatory role of miRNAs as the intensive interactions between miRNAs and hub isoforms. We also found that the isoforms from the same gene can be selectively targeted by miRNAs, and we further explored the functional divergence of these isoforms in both structure and expression. Collectively, this study on the interactions between the miRNAs and isoforms, and functional divergence of isoforms can facilitate our current understanding of the complexity of gene regulation in plants.\u003c/p\u003e"},{"header":"Methods","content":"\u003cp\u003e\u003cstrong\u003ePlant material and RNA extraction\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe tissues including leaf, petiole, petal and mature anther were harvested from\u003cem\u003e N. nucifera\u003c/em\u003e var. \u0026lsquo;China Antique\u0026rsquo; in Wuhan Botanical Garden, CAS (114\u003csup\u003eo\u003c/sup\u003e30\u0026rsquo;E, 30\u003csup\u003eo\u003c/sup\u003e60\u0026rsquo;N). The unpollinated and pollinated carpel were collected on a blooming day according to our previous study [27]. The samples were immediately frozen in liquid nitrogen and the total RNAs were extracted from each sample using RNAprep pure Plant Kit (TIANGEN). And the RNA was only used for RT-PCR validation experiments in this study. Meanwhile, the corresponding RNA-seq datasets from our previous study were downloaded from the NCBI Short Read Archive (SRA) accession number PRJNA503979 and PRJNA492157 [59].\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eSmall RNA-seq analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe small RNA-seq dataset was downloaded from the SRA accession number SRX1591010 [27]. The latest plant miRNA sequences were downloaded from miRbase 22.0 for miRNA annotation. For each sample, the high-quality reads without adaptors were filtered to retain sRNAs reads length from 18bp to 29bp. Before miRNA identification, the retained reads were further processed using the miRDeep2 [55]. After quality filtering, the Bowtie was used to map the high-quality reads to the reference genome (nelumbo.biocloud.net) with zero mismatches [29, 56]. The mapped small RNAs were guided by the known precursor miRNA dataset (miRbase 22.0) to identify potential miRNA precursors in the lotus genome, allowing the duplicate (hairpin) loci less than five. The hairpin structures and the aligned small RNAs were processed as described in the miRDeep-P package [57]. The overlapping (redundant) precursors were removed. Subsequently, we identified TE-related miRNAs as those pre-miRNAs overlapping with transposable elements in the lotus genome. To quantify the expression level of miRNA, the read counts were converted into TPM [58]. Differential expression analysis was performed using edgeR. Differentially expressed miRNAs between two samples were defined as those with Benjam\u0026iacute;n-Hochberg false discovery rate (FDR) less than 0.05 (FDR \u0026lt; 0.05) and log2 of fold_change more than one.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eIsoform identification and quantification\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo identify the isoforms, the PacBio full-length transcripts were downloaded and mapped to the new lotus reference genome as previously described [59]. Illumina RNA-Seq from petiole, leaf, petal, anther, unpollinated carpel and pollinated carpel were also downloaded to quantify isoform expression. In brief, the high-quality reads were mapped to the reference by HISTA2 v2.1.0, and then the FPKM (fragments per kilobase per million) value of the genes and isoforms were calculated by StringTie v1.3.5 using a combination of Illumina and full-length transcript based annotations. The differentially expressed genes and isoforms were identified by the edgeR package using the same threshold used for miRNA differential expression analysis.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePrediction of miRNA target genes and isoforms\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eComparing with the animal miRNA, most miRNAs and their target mRNAs in the plant have near-perfect sequence complementarity [60]. The miRNAs target isoforms were predicted by TargetFinder [61], with a strict screening criterion. Briefly, FASTA search was applied to find the potential targets, while a position-dependent mispair penalty system was used to score the target sequence after assessment of penalties for mismatches, bulges, gaps (+1 per position) and G:U pairs (+0.5 per position)[62]. Penalties were doubled if the mismatch, bulge, gap, or G:U pairs occurred at position 2 to 13 relative to the 5\u0026rsquo; end of the miRNA. Only one single nucleotide bulge or singe nucleotide gap was allowed while other parameters were set at default. Finally, only predicted targets with scores of three or less were retained.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eGene Ontology, KEGG pathway analyses, and ortholog identification\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe GOseq package based on Wallenius non-central hypergeometric distribution [63] was used to identify the significantly enriched GO terms and KEGG pathways for each gene set, including miRNA target genes and the differentially expressed genes. The KEGG pathway annotation was performed using the target gene as queries in GhostKOALA (https://www.kegg.jp/). The orthologs among Arabidopsis (TAIR11), rice (http://rice.plantbiology.msu.edu/, version 7.0) and lotus were identified using the OrthoMCL with an e-value \u0026lt; 1e-15 and an inflation parameter of 2.0.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCo-expression network analysis of transcript isoforms by WCGNA\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo determine whether isoform expression is correlated with the tissues, the co-expression networks were constructed based on 15 RNA-seq samples by using WGCNA (v1.0) in R[64]. These 15 samples include five leaf samples, two petioles, two mature anthers, two pollinated carpels, and unpollinated carpels (PRJNA503979 and PRJNA492157). To filter out the silent or constantly expressed isoforms, those with an averaged FPKM \u0026gt; 0.1 and coefficient of variation (C.V.) of FPKM \u0026gt; 2 were retained for the subsequent pipeline of WGCNA. According to Topological Overlap Matrix [65], the transcripts were first clustered hierarchically. The transcripts were assigned to nine co-expression modules using a bottom-up algorithm known as the dynamic hybrid cut method, which was named after randomly assigned colors. Modules were identified under the minimum module size of 600, and the threshold of the module similarity was set at 0.5. The kME of isoform, which measures the correlation between isoforms and modules, was calculated. The top 150 isoforms, which were most correlated with each corresponding module, were defined as hub isoforms according to the WGCNA protocol.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eRT-qPCR analysis of sampled target genes\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe cDNA libraries were constructed using 1 \u0026micro;g total RNA from five samples and they were diluted to 100 \u0026micro;L before performing RT-PCR. The miRNA targeted gene-specific primers were designed by using Primer Premier (v.5.0) (Additional file 2: Table S8). The qRT-PCR reactions were performed on QuantStudio 6 Flex (Life Technologies, USA) in a final volume of 10 \u0026micro;g containing 1 \u0026micro;g cDNA. The reaction procedure was initiated at 95\u0026deg;C for 10 min, followed by 42 cycles of 95\u0026deg;C for 15 s, 60\u0026deg;C for 30 s, and 72\u0026deg;C for 30 s with fluorescence detection. The miRNA targeted gene expressions were analyzed using the 2\u003csup\u003e-\u003c/sup\u003e\u003csup\u003e△△\u003c/sup\u003e\u003csup\u003eCT \u003c/sup\u003emethod, with the lotus \u003cem\u003e\u0026beta;-Actin\u003c/em\u003e gene used as an internal standard.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsent for pulication\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAvailability of data and materials\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll data that support the findings of this study are included in the manuscript and supplementary information.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no competing interests.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis work was supported by grants from the Strategic Priority Research Program CAS (No. XDB31000000), the National Natural Science Foundation of China (No. 31570220, No. 31870208 and No. 31700197), Youth Innovation Promotion Association of Chinese Academy of Sciences (No. 2019335) and Bureau of Landscaping and Forestry of Wuhan Municipality (No.2018-28). The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthors\u0026rsquo; Contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTS and JC conceived and designed the research. YZ and XY participated in the preparation of the reagents and materials. YZ conducted the experiments. TS and YZ analyzed the data. YZ and RSR wrote the manuscript. All the authors read and approved the manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgements\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe thank all individuals who participated in this study. We thank our friend Owen Edward (USA) for proofreading the manuscript.\u003c/p\u003e"},{"header":"Abbreviations","content":"\u003cp\u003eAS: Alternative Splicing\u003c/p\u003e\n\u003cp\u003eDEG: Differentially Expressed Gene\u003c/p\u003e\n\u003cp\u003eDEMTG: Differentially Expressed miRNA Targeted Gene\u003c/p\u003e\n\u003cp\u003eFDR: False Discovery Rate\u003c/p\u003e\n\u003cp\u003eFPKM: Fragments Per Kilobase of transcript per Million fragments mapped\u003c/p\u003e\n\u003cp\u003eGIDDM: Genes have their Isoforms Divided into Different Modules\u003c/p\u003e\n\u003cp\u003eGO: Gene Ontology\u003c/p\u003e\n\u003cp\u003eKEEG: Kyoto Encyclopedia of Gene and Genomes\u003c/p\u003e\n\u003cp\u003eMIRNAs: miRNA genes\u003c/p\u003e\n\u003cp\u003eTPM: Transcripts per million reads\u003c/p\u003e\n\u003cp\u003eWGCNA: Weight Gene Co-expression Network Analysis\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eReinhart BJ, Weinstein EG, Rhoades MW, Bartel B, Bartel DP. MicroRNAs in plants. Genes \u0026amp; development. 2002;16(13):1616-1626.\u003c/li\u003e\n\u003cli\u003eLee Y, Kim M, Han J, Yeom KH, Lee S, Baek SH, Kim VN. MicroRNA genes are transcribed by RNA polymerase II. EMBO journal. 2004;23(20):4051-4060.\u003c/li\u003e\n\u003cli\u003eCui J, You C, Chen X. The evolution of microRNAs in plants. Current opinion in plant biology\u003cem\u003e. \u003c/em\u003e2017;35:61-67.\u003c/li\u003e\n\u003cli\u003eLi C, Zhang B. MicroRNAs in Control of Plant Development. Journal of cellular physiology\u003cem\u003e. \u003c/em\u003e2016;231(2):303-313.\u003c/li\u003e\n\u003cli\u003eFerdous J, Hussain SS, Shi BJ. Role of microRNAs in plant drought tolerance. Plant biotechnology journal. 2015;13(3):293-305.\u003c/li\u003e\n\u003cli\u003eYao S, Yang Z, Yang R, Huang Y, Guo G, Kong X, Lan Y, Zhou T, Wang H, Wang W, et al. Transcriptional Regulation of miR528 by OsSPL9 Orchestrates Antiviral Response in Rice. Molecular plant. 2019;12(8):1114-1122.\u003c/li\u003e\n\u003cli\u003eZhang S, Yan S, Zhao J, Xiong H, An P, Wang J, Zhang H, Zhang L. Identification of miRNAs and their target genes in Larix olgensis and verified of differential expression miRNAs. BMC plant biology\u003cem\u003e. \u003c/em\u003e2019;19(1):247.\u003c/li\u003e\n\u003cli\u003eSong X, Li Y, Cao X, Qi Y. MicroRNAs and Their Regulatory Roles in Plant-Environment Interactions. Annual review of plant biology\u003cem\u003e. \u003c/em\u003e2019;70:489-525.\u003c/li\u003e\n\u003cli\u003eWang JW. Regulation of flowering time by the miR156-mediated age pathway. Journal of experimental botany\u003cem\u003e. \u003c/em\u003e2014;65(17):4723-4730.\u003c/li\u003e\n\u003cli\u003eChen X. A microRNA as a translational repressor of APETALA2 in Arabidopsis flower development. Science\u003cem\u003e. \u003c/em\u003e2004;303(5666):2022-2025.\u003c/li\u003e\n\u003cli\u003ePokoo R, Ren S, Wang Q, Motes CM, Hernandez TD, Ahmadi S, Monteros MJ, Zheng Y, Sunkar R. Genotype- and tissue-specific miRNA profiles and their targets in three alfalfa (\u003cem\u003eMedicago sativa L\u003c/em\u003e) genotypes. BMC genomics\u003cem\u003e. \u003c/em\u003e2018;19(10):913.\u003c/li\u003e\n\u003cli\u003eCelik O, Akdas EY. Tissue-specific transcriptional regulation of seven heavy metal stress-responsive miRNAs and their putative targets in nickel indicator castor bean (\u003cem\u003eR. communis L.\u003c/em\u003e) plants. Ecotoxicology and environmental safety\u003cem\u003e. \u003c/em\u003e2019;170:682-690.\u003c/li\u003e\n\u003cli\u003eTaller D, Balint J, Gyula P, Nagy T, Barta E, Baksa I, Szittya G, Taller J, Havelda Z. Expansion of Capsicum annuum fruit is linked to dynamic tissue-specific differential expression of miRNA and siRNA profiles. PloS one\u003cem\u003e. \u003c/em\u003e2018;13(7):e0200207.\u003c/li\u003e\n\u003cli\u003eChen Q, Deng B, Gao J, Zhao Z, Chen Z, Song S, Wang L, Zhao L, Xu W, Zhang C, et al. Comparative Analysis of miRNA Abundance Revealed the Function of Vvi-miR828 in Fruit Coloring in Root Restriction Cultivation Grapevine (\u003cem\u003eVitis vinifera L.\u003c/em\u003e). International journal of molecular sciences\u003cem\u003e. \u003c/em\u003e2019;20(16).\u003c/li\u003e\n\u003cli\u003eLi XM, Sang YL, Zhao XY, Zhang XS. High-throughput sequencing of small RNAs from pollen and silk and characterization of miRNAs as candidate factors involved in pollen-silk interactions in maize. PloS one\u003cem\u003e. \u003c/em\u003e2013;8(8):e72852.\u003c/li\u003e\n\u003cli\u003eCuraba J, Singh MB, Bhalla PL. miRNAs in the crosstalk between phytohormone signalling pathways. Journal of experimental botany\u003cem\u003e. \u003c/em\u003e2014;65(6):1425-1438.\u003c/li\u003e\n\u003cli\u003eJin D, Wang Y, Zhao Y, Chen M. MicroRNAs and their cross-talks in plant development. Journal of genetics and genomics\u003cem\u003e. \u003c/em\u003e2013;40(4):161-170.\u003c/li\u003e\n\u003cli\u003eFilichkin S, Priest HD, Megraw M, Mockler TC. Alternative splicing in plants: directing traffic at the crossroads of adaptation and environmental stress. Current opinion in plant biology\u003cem\u003e. \u003c/em\u003e2015;24:125-135.\u003c/li\u003e\n\u003cli\u003eQiao D, Yang C, Chen J, Guo Y, Li Y, Niu S, Cao K, Chen Z. Comprehensive identification of the full-length transcripts and alternative splicing related to the secondary metabolism pathways in the tea plant (\u003cem\u003eCamellia sinensis\u003c/em\u003e). Scientific reports\u003cem\u003e. \u003c/em\u003e2019;9(1):2709.\u003c/li\u003e\n\u003cli\u003eWang L, Jiang X, Wang L, Wang W, Fu C, Yan X, Geng X. A survey of transcriptome complexity using PacBio single-molecule real-time analysis combined with Illumina RNA sequencing for a better understanding of ricinoleic acid biosynthesis in Ricinus communis. BMC genomics\u003cem\u003e. \u003c/em\u003e2019;20(1):456.\u003c/li\u003e\n\u003cli\u003eWang J, Deng Y, Zhou Y, Liu D, Yu H, Zhou Y, Lv J, Ou L, Li X, Ma Y, et al. Full-length mRNA sequencing and gene expression profiling reveal broad involvement of natural antisense transcript gene pairs in pepper development and response to stresses. Plant journal. 2019;99(4):763-783.\u003c/li\u003e\n\u003cli\u003eSatyawan D, Kim MY, Lee SH. Stochastic alternative splicing is prevalent in mungbean (\u003cem\u003eVigna radiata\u003c/em\u003e). Plant biotechnology journal. 2017;15(2):174-182.\u003c/li\u003e\n\u003cli\u003eWang M, Wang P, Liang F, Ye Z, Li J, Shen C, Pei L, Wang F, Hu J, Tu L, et al. A global survey of alternative splicing in allopolyploid cotton: landscape, complexity and regulation. New phytologist\u003cem\u003e. \u003c/em\u003e2018;217(1):163-178.\u003c/li\u003e\n\u003cli\u003eThatcher SR, Danilevskaya ON, Meng X, Beatty M, Zastrow-Hayes G, Harris C, Van Allen B, Habben J, Li B. Genome-Wide Analysis of Alternative Splicing during Development and Drought Stress in Maize. Plant physiology\u003cem\u003e. \u003c/em\u003e2016;170(1):586-599.\u003c/li\u003e\n\u003cli\u003eChen J, Li L. Multiple Regression Analysis Reveals MicroRNA Regulatory Networks in Oryza sativa under Drought Stress. International journal of genomics\u003cem\u003e. \u003c/em\u003e2018; 9395261.\u003c/li\u003e\n\u003cli\u003eWu C, Li X, Guo S, Wong SM. Analyses of RNA-Seq and sRNA-Seq data reveal a complex network of anti-viral defense in TCV-infected Arabidopsis thaliana. Scientific reports\u003cem\u003e. \u003c/em\u003e2016;6:36007.\u003c/li\u003e\n\u003cli\u003eShi T, Wang K, Yang P. The evolution of plant microRNAs: insights from a basal eudicot sacred lotus. Plant journal. 2017;89(3):442-457.\u003c/li\u003e\n\u003cli\u003eZhang Y, Nyong AT, Shi T, Yang P. The complexity of alternative splicing and landscape of tissue-specific expression in lotus (\u003cem\u003eNelumbo nucifera\u003c/em\u003e) unveiled by Illumina- and single-molecule real-time-based RNA-sequencing. DNA research. 2019;26(4):301-311.\u003c/li\u003e\n\u003cli\u003eShi T, Rahmani RS, Gugger PF, Wang M, Li H, Zhang Y, Li Z, Wang Q, Van de Peer Y, Marchal K et al. Distinct expression and methylation patterns for genes with different fates following a single whole-genome duplication in flowering plant. Molecular biology and evolution 2020.\u003c/li\u003e\n\u003cli\u003eLi Y, Li C, Xia J, Jin Y. Domestication of transposable elements into MicroRNA genes in plants. PloS one\u003cem\u003e. \u003c/em\u003e2011;6(5):e19212.\u003c/li\u003e\n\u003cli\u003eLu S. De novo origination of MIRNAs through generation of short inverted repeats in target genes. RNA biology\u003cem\u003e. \u003c/em\u003e2019;16(6):846-859.\u003c/li\u003e\n\u003cli\u003eBarrera-Figueroa BE, Gao L, Wu Z, Zhou X, Zhu J, Jin H, Liu R, Zhu JK. High throughput sequencing reveals novel and abiotic stress-regulated microRNAs in the inflorescences of rice. BMC plant biology\u003cem\u003e. \u003c/em\u003e2012;12:132.\u003c/li\u003e\n\u003cli\u003eRoberts JT, Cardin SE, Borchert GM. Burgeoning evidence indicates that microRNAs were initially formed from transposable element sequences. Mobile genetic elements\u003cem\u003e. \u003c/em\u003e2014;4:e29255.\u003c/li\u003e\n\u003cli\u003eYin H, Fan Z, Li X, Wang J, Liu W, Wu B, Ying Z, Liu L, Liu Z, Li J. Phylogenetic tree-informed microRNAome analysis uncovers conserved and lineage-specific miRNAs in Camellia during floral organ development. Journal of experimental botany\u003cem\u003e. \u003c/em\u003e2016;67(9):2641-2653.\u003c/li\u003e\n\u003cli\u003eSaminathan T, Bodunrin A, Singh NV, Devarajan R, Nimmakayala P, Jeff M, Aradhya M, Reddy UK. Genome-wide identification of microRNAs in pomegranate (\u003cem\u003ePunica granatum L.\u003c/em\u003e) by high-throughput sequencing. BMC plant biology\u003cem\u003e. \u003c/em\u003e2016;16(1):122.\u003c/li\u003e\n\u003cli\u003eSamad AFA, Sajad M, Nazaruddin N, Fauzi IA, Murad AMA, Zainal Z, Ismail I. MicroRNA and Transcription Factor: Key Players in Plant Regulatory Network. Frontiers in plant science\u003cem\u003e. \u003c/em\u003e2017;8:565.\u003c/li\u003e\n\u003cli\u003eKoyama T, Sato F, Ohme-Takagi M. Roles of miR319 and TCP Transcription Factors in Leaf Development. Plant physiology\u003cem\u003e. \u003c/em\u003e2017; 175(2):874-885.\u003c/li\u003e\n\u003cli\u003eYu JY, Zhang ZG, Huang SY, Han X, Wang XY, Pan WJ, Qin HT, Qi HD, Yin ZG, Qu KX, et al. Analysis of miRNAs Targeted Storage Regulatory Genes during Soybean Seed Development Based on Transcriptome Sequencing. Genes\u003cem\u003e. \u003c/em\u003e2019;10(6).\u003c/li\u003e\n\u003cli\u003eChen H, Yang Q, Chen K, Zhao S, Zhang C, Pan R, Cai T, Deng Y, Wang X, Chen Y, et al. Integrated microRNA and transcriptome profiling reveals a miRNA-mediated regulatory network of embryo abortion under calcium deficiency in peanut (\u003cem\u003eArachis hypogaea L.\u003c/em\u003e). BMC genomics\u003cem\u003e. \u003c/em\u003e2019;20(1):392.\u003c/li\u003e\n\u003cli\u003eXie F, Zhang B. microRNA evolution and expression analysis in polyploidized cotton genome. Plant biotechnology journal. 2015;13(3):421-434.\u003c/li\u003e\n\u003cli\u003ePiriyapongsa J, Jordan IK. Dual coding of siRNAs and miRNAs by plant transposable elements. RNA\u003cem\u003e. \u003c/em\u003e2008;14(5):814-821.\u003c/li\u003e\n\u003cli\u003eYu M, Carver BF, Yan L. TamiR1123 originated from a family of miniature inverted-repeat transposable elements (MITE) including one inserted in the Vrn-A1a promoter in wheat. Plant science\u003cem\u003e. \u003c/em\u003e2014;215-216:117-123.\u003c/li\u003e\n\u003cli\u003eWang R, Yang Z, Fei Y, Feng J, Zhu H, Huang F, Zhang H, Huang J. Construction and analysis of degradome-dependent microRNA regulatory networks in soybean. BMC genomics\u003cem\u003e. \u003c/em\u003e2019;20(1):534.\u003c/li\u003e\n\u003cli\u003eChen J, Su P, Chen P, Li Q, Yuan X, Liu Z. Insights into the cotton anther development through association analysis of transcriptomic and small RNA sequencing. BMC plant biology\u003cem\u003e. \u003c/em\u003e2018;18(1):154.\u003c/li\u003e\n\u003cli\u003eLuo Y, Guo Z, Li L. Evolutionary conservation of microRNA regulatory programs in plant flower development. Developmental biology\u003cem\u003e. \u003c/em\u003e2013;380(2):133-144.\u003c/li\u003e\n\u003cli\u003eJu Q, Zhao YJ, Dong Y, Cheng C, Zhang S, Yang Y, Li P, Ge D, Sun B. Identification of a miRNA-mRNA network associated with lymph node metastasis in colorectal cancer. Oncology letters\u003cem\u003e. \u003c/em\u003e2019;18(2):1179-1188.\u003c/li\u003e\n\u003cli\u003eZhou XG, Huang XL, Liang SY, Tang SM, Wu SK, Huang TT, Mo ZN, Wang QY. Identifying miRNA and gene modules of colon cancer associated with pathological stage by weighted gene co-expression network analysis. OncoTargets and therapy. 2018;11:2815-2830.\u003c/li\u003e\n\u003cli\u003eTai Y, Liu C, Yu S, Yang H, Sun J, Guo C, Huang B, Liu Z, Yuan Y, Xia E, et al. Gene co-expression network analysis reveals coordinated regulation of three characteristic secondary biosynthetic pathways in tea plant (\u003cem\u003eCamellia sinensis\u003c/em\u003e). BMC genomics\u003cem\u003e. \u003c/em\u003e2018;19(1):616.\u003c/li\u003e\n\u003cli\u003eYu Y, Shi J, Li X, Liu J, Geng Q, Shi H, Ke Y, Sun Q. Transcriptome analysis reveals the molecular mechanisms of the defense response to gray leaf spot disease in maize. BMC genomics\u003cem\u003e. \u003c/em\u003e2018;19(1):742.\u003c/li\u003e\n\u003cli\u003eAbascal F, Tress ML, Valencia A. The evolutionary fate of alternatively spliced homologous exons after gene duplication. Genome biology and evolution. 2015;7(6):1392-1403.\u003c/li\u003e\n\u003cli\u003eChen H, Li Z, Xiong L. A plant microRNA regulates the adaptation of roots to drought stress. FEBS letters\u003cem\u003e. \u003c/em\u003e2012;586(12):1742-1747.\u003c/li\u003e\n\u003cli\u003eXu J, Li J, Cui L, Zhang T, Wu Z, Zhu PY, Meng YJ, Zhang KJ, Yu XQ, Lou QF, et al. New insights into the roles of cucumber TIR1 homologs and miR393 in regulating fruit/seed set development and leaf morphogenesis. BMC plant biology\u003cem\u003e. \u003c/em\u003e2017;17(1):130.\u003c/li\u003e\n\u003cli\u003eMarin E, Jouannet V, Herz A, Lokerse AS, Weijers D, Vaucheret H, Nussaume L, Crespi MD, Maizel A. miR390, Arabidopsis TAS3 tasiRNAs, and their AUXIN RESPONSE FACTOR targets define an autoregulatory network quantitatively regulating lateral root growth. Plant cell\u003cem\u003e. \u003c/em\u003e2010;22(4):1104-1117.\u003c/li\u003e\n\u003cli\u003eAn W, Gong W, He S, Pan Z, Sun J, Du X. MicroRNA and mRNA expression profiling analysis revealed the regulation of plant height in Gossypium hirsutum. BMC genomics\u003cem\u003e. \u003c/em\u003e2015;16:886.\u003c/li\u003e\n\u003cli\u003eFriedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic acids research\u003cem\u003e. \u003c/em\u003e2012;40(1):37-52.\u003c/li\u003e\n\u003cli\u003eLangmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome biology\u003cem\u003e. \u003c/em\u003e2009;10(3):R25.\u003c/li\u003e\n\u003cli\u003eYang X, Li L. miRDeep-P: a computational tool for analyzing the microRNA transcriptome in plants. Bioinformatics\u003cem\u003e. \u003c/em\u003e2011;27(18):2614-2615.\u003c/li\u003e\n\u003cli\u003eFahlgren N, Howell MD, Kasschau KD, Chapman EJ, Sullivan CM, Cumbie JS, Givan SA, Law TF, Grant SR, Dangl JL, et al. High-throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of MIRNA genes. PloS one\u003cem\u003e. \u003c/em\u003e2007;2(2):e219.\u003c/li\u003e\n\u003cli\u003eAhsan MU, Hayward A, Irihimovitch V, Fletcher S, Tanurdzic M, Pocock A, Beveridge CA, Mitter N. Juvenility and Vegetative Phase Transition in Tropical/Subtropical Tree Crops. Frontiers in plant science\u003cem\u003e. \u003c/em\u003e2019;10:729.\u003c/li\u003e\n\u003cli\u003eParizotto EA, Dunoyer P, Rahm N, Himber C, Voinnet O. In vivo investigation of the transcription, processing, endonucleolytic activity, and functional relevance of the spatial distribution of a plant miRNA. Genes \u0026amp; development\u003cem\u003e. \u003c/em\u003e2004;18(18):2237-2242.\u003c/li\u003e\n\u003cli\u003eBo X, Wang S. TargetFinder: a software for antisense oligonucleotide target site selection based on MAST and secondary structures of target mRNA. Bioinformatics\u003cem\u003e. \u003c/em\u003e2005;21(8):1401-1402.\u003c/li\u003e\n\u003cli\u003eAllen E, Xie Z, Gustafson AM, Carrington JC. microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell\u003cem\u003e. \u003c/em\u003e2005;121(2):207-221.\u003c/li\u003e\n\u003cli\u003eYoung MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome biology\u003cem\u003e. \u003c/em\u003e2010;11(2):R14.\u003c/li\u003e\n\u003cli\u003eLangfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC bioinformatics\u003cem\u003e. \u003c/em\u003e2008;9:559.\u003c/li\u003e\n\u003cli\u003eZhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Statistical applications in genetics and molecular biology\u003cem\u003e. \u003c/em\u003e2005;4:Article17.\u003c/li\u003e\n\u003c/ol\u003e"},{"header":"Table","content":"\u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:justify;font-size:14px;font-family:DengXian;line-height:200%;\"\u003e\u003cspan style='font-size: 15px; line-height: 200%; font-family: \"Times New Roman\", serif; color: rgb(0, 0, 0);'\u003eTable.1 Summary of high-quality reads mapped to miRBase.\u003c/span\u003e\u003c/p\u003e\n\u003ctable style=\"width: 4.4e+2pt;border: none;border-collapse:collapse;\"\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd style=\"width:68.25pt;border-top:windowtext;border-left:white;border-bottom:windowtext;border-right:white;border-style:solid;border-width: 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:53.75pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size:15px;font-family:\"Times New Roman\",serif;'\u003eSample\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:87.2pt;border-top:solid windowtext 1.0pt;border-left:none;border-bottom:solid windowtext 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:53.75pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003eHigh-quality reads\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:120.5pt;border-top:solid windowtext 1.0pt;border-left:none;border-bottom:solid windowtext 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:53.75pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003eReads with at least one alignment in miRBase\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:168.05pt;border-top:solid windowtext 1.0pt;border-left:none;border-bottom:solid windowtext 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:53.75pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003eReads without alignment in miRBase\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width:68.25pt;border:solid white 1.0pt;border-top:none;background:white;padding:.75pt .75pt 0in .75pt;height:32.25pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003eAnther\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:87.2pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:32.25pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e6826041\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:120.5pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:32.25pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e11337 (0.17%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:168.05pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:32.25pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e6814704 (99.83%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width:68.25pt;border:solid white 1.0pt;border-top:none;background:white;padding:.75pt .75pt 0in .75pt;height:36.45pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003eUnpollinated carpel\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:87.2pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:36.45pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e7164301\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:120.5pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:36.45pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e10947(0.15%\u003c/span\u003e\u003cspan style=\"font-size: 15px; font-family: SimSun;\"\u003e）\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:168.05pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:36.45pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e7153354(99.85%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width:68.25pt;border:solid white 1.0pt;border-top:none;background:white;padding:.75pt .75pt 0in .75pt;height:36.6pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003ePollinated carpel\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:87.2pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:36.6pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e8943060\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:120.5pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:36.6pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e13735 (0.15%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:168.05pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:36.6pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e8929325 (99.85%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width:68.25pt;border:solid white 1.0pt;border-top:none;background:white;padding:.75pt .75pt 0in .75pt;height:39.8pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003ePetiole\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:87.2pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:39.8pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e1472529\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:120.5pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:39.8pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e6118 (0.42%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:168.05pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:39.8pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e1466411 (99.58%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width:68.25pt;border:solid white 1.0pt;border-top:none;background:white;padding:.75pt .75pt 0in .75pt;height:31.2pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003eLeaf\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:87.2pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:31.2pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e1707006\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:120.5pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:31.2pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e7964 (0.47%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:168.05pt;border-top:none;border-left:none;border-bottom:solid white 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:31.2pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e1699042 (99.53%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width:68.25pt;border-top:none;border-left:solid white 1.0pt;border-bottom:solid windowtext 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:33.75pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003ePetal\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:87.2pt;border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:33.75pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e1518064\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:120.5pt;border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:33.75pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e5984 (0.39%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width:168.05pt;border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;border-right:solid white 1.0pt;background:white;padding:.75pt .75pt 0in .75pt;height:33.75pt;\"\u003e\n \u003cp style=\"margin:0in;margin-bottom:.0001pt;text-align:center;font-size:14px;font-family:DengXian;\"\u003e\u003cspan style=\"color: rgb(0, 0, 0);\"\u003e\u003cspan style='font-size: 15px; font-family: \"Times New Roman\", serif;'\u003e1512080 (99.61%)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n\u003c/table\u003e"},{"header":"Additional Files","content":"\u003cp\u003eAdditional file 1: Figure S1-S18\u003c/p\u003e\n\u003cp\u003eAdditional file 2: Table S1-S8\u003c/p\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"bmc-genomics","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"gics","sideBox":"Learn more about [BMC Genomics](http://bmcgenomics.biomedcentral.com/)","snPcode":"","submissionUrl":"https://www.editorialmanager.com/gics","title":"BMC Genomics","twitterHandle":"#BMCGenomics","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"em","reportingPortfolio":"BMC Series","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"microRNA, transcript isoforms, co-expression network, sacred lotus","lastPublishedDoi":"10.21203/rs.2.19330/v3","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.2.19330/v3","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003e\u003cstrong\u003eBackground: \u003c/strong\u003eGene expression is complex and regulated by multiple molecular mechanisms, such as miRNA-mediated gene inhibition and alternative-splicing of pre-mRNAs. However, the coordination of interaction between miRNAs with different splicing isoforms, and the change of splicing isoform in response to different cellular environments are largely unexplored in plants. In this study, we analyzed the miRNA and mRNA transcriptome from lotus (\u003cem\u003eNelumbo nucifera\u003c/em\u003e), an economically important flowering plant.\u003c/p\u003e\u003cp\u003e\u003cstrong\u003eResults: \u003c/strong\u003eThrough RNA-seq analyses on miRNAs and their target genes (isoforms) among six lotus tissues, expression of most miRNAs seem to be negatively correlated with their targets and tend to be tissue-specific. Further, our results showed that preferential interactions between miRNAs and hub gene isoforms in one coexpression module which is highly correlated with leaf. Intriguingly, for many genes, their corresponding isoforms were assigned to different co-expressed modules, and they exhibited more divergent mRNA structures including presence and absence of miRNA binding sites, suggesting functional divergence for many isoforms is escalated by both structural and expression divergence. Further detailed functional enrichment analysis of miRNA targets revealed that miRNAs are involved in the regulation of lotus growth and development by regulating plant hormone-related pathway genes.\u003c/p\u003e\u003cp\u003e\u003cstrong\u003eConclusions: \u003c/strong\u003eTaken together, our comprehensive analyses of miRNA and mRNA transcriptome elucidate the coordination of interaction between miRNAs and different splicing isoforms, and highlight the functional divergence of many transcript isoforms from the same locus in lotus.\u003c/p\u003e","manuscriptTitle":"Integrative expression network analysis of microRNA and gene isoforms in sacred lotus","msid":"","msnumber":"","nonDraftVersions":[{"code":3,"date":"2020-05-26 16:07:10","doi":"10.21203/rs.2.19330/v3","editorialEvents":[{"type":"communityComments","content":0},{"type":"editorInvitedReview","content":"","date":"2020-06-08T12:00:00+00:00","index":1,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"editorInvitedReview","content":"","date":"2020-06-01T12:00:00+00:00","index":2,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"reviewerAgreed","content":"","date":"2020-05-20T12:00:00+00:00","index":4,"fulltext":""},{"type":"editorInvitedReview","content":"","date":"2020-05-20T12:00:00+00:00","index":4,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"editorInvitedReview","content":"","date":"2020-05-20T12:00:00+00:00","index":3,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"reviewerAgreed","content":"","date":"2020-05-19T12:00:00+00:00","index":2,"fulltext":""},{"type":"reviewerAgreed","content":"","date":"2020-05-19T12:00:00+00:00","index":3,"fulltext":""},{"type":"reviewersInvited","content":"","date":"2020-05-18T12:00:00+00:00","index":"","fulltext":""},{"type":"reviewerAgreed","content":"","date":"2020-05-18T12:00:00+00:00","index":1,"fulltext":""},{"type":"editorAssigned","content":"","date":"2020-05-13T12:00:00+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2020-05-12T12:00:00+00:00","index":"","fulltext":""},{"type":"editorInvited","content":"","date":"2020-05-12T12:00:00+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"bmc-genomics","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"gics","sideBox":"Learn more about [BMC Genomics](http://bmcgenomics.biomedcentral.com/)","snPcode":"","submissionUrl":"https://www.editorialmanager.com/gics","title":"BMC Genomics","twitterHandle":"#BMCGenomics","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"em","reportingPortfolio":"BMC Series","inReviewEnabled":true,"inReviewRevisionsEnabled":true}},{"code":2,"date":"2020-03-18 17:08:10","doi":"10.21203/rs.2.19330/v2","editorialEvents":[{"type":"communityComments","content":0},{"type":"editorInvitedReview","content":"","date":"2020-04-14T12:00:00+00:00","index":4,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"decision","content":"Major revision","date":"2020-04-14T12:00:00+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2020-04-08T12:00:00+00:00","index":2,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"editorInvitedReview","content":"","date":"2020-04-05T12:00:00+00:00","index":3,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"reviewerAgreed","content":"","date":"2020-03-28T12:00:00+00:00","index":4,"fulltext":""},{"type":"reviewerAgreed","content":"","date":"2020-03-26T12:00:00+00:00","index":3,"fulltext":""},{"type":"reviewerAgreed","content":"","date":"2020-03-25T12:00:00+00:00","index":2,"fulltext":""},{"type":"reviewersInvited","content":"","date":"2020-03-23T12:00:00+00:00","index":"","fulltext":""},{"type":"reviewerAgreed","content":"","date":"2020-03-23T12:00:00+00:00","index":1,"fulltext":""},{"type":"editorInvitedReview","content":"","date":"2020-03-23T12:00:00+00:00","index":1,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"editorAssigned","content":"","date":"2020-03-17T12:00:00+00:00","index":"","fulltext":""},{"type":"editorInvited","content":"","date":"2020-03-16T12:00:00+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2020-03-12T12:00:00+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"bmc-genomics","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"gics","sideBox":"Learn more about [BMC Genomics](http://bmcgenomics.biomedcentral.com/)","snPcode":"","submissionUrl":"https://www.editorialmanager.com/gics","title":"BMC Genomics","twitterHandle":"#BMCGenomics","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"em","reportingPortfolio":"BMC Series","inReviewEnabled":true,"inReviewRevisionsEnabled":true}},{"code":1,"date":"2019-12-19 20:21:48","doi":"10.21203/rs.2.19330/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Major revision","date":"2020-01-25T12:00:00+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2020-01-14T12:00:00+00:00","index":2,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"editorInvitedReview","content":"","date":"2020-01-14T12:00:00+00:00","index":1,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"editorInvitedReview","content":"","date":"2020-01-14T12:00:00+00:00","index":3,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"editorInvitedReview","content":"","date":"2020-01-13T12:00:00+00:00","index":5,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"editorInvitedReview","content":"","date":"2020-01-13T12:00:00+00:00","index":4,"fulltext":"Recommendation: Reviewer's comments unavailable due to the journal's policy.\n"},{"type":"reviewerAgreed","content":"","date":"2019-12-30T12:00:00+00:00","index":5,"fulltext":""},{"type":"reviewerAgreed","content":"","date":"2019-12-29T12:00:00+00:00","index":2,"fulltext":""},{"type":"reviewerAgreed","content":"","date":"2019-12-29T12:00:00+00:00","index":3,"fulltext":""},{"type":"reviewerAgreed","content":"","date":"2019-12-29T12:00:00+00:00","index":4,"fulltext":""},{"type":"reviewerAgreed","content":"","date":"2019-12-23T12:00:00+00:00","index":1,"fulltext":""},{"type":"reviewersInvited","content":"","date":"2019-12-18T12:00:00+00:00","index":"","fulltext":""},{"type":"editorInvited","content":"","date":"2019-12-16T12:00:00+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2019-12-12T12:00:00+00:00","index":"","fulltext":""},{"type":"submitted","content":"","date":"2019-12-11T12:00:00+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2019-12-11T12:00:00+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"bmc-genomics","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"gics","sideBox":"Learn more about [BMC Genomics](http://bmcgenomics.biomedcentral.com/)","snPcode":"","submissionUrl":"https://www.editorialmanager.com/gics","title":"BMC Genomics","twitterHandle":"#BMCGenomics","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"em","reportingPortfolio":"BMC Series","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"1e137d90-2345-4b7c-83c7-b760b79f63d9","owner":[],"postedDate":"May 26th, 2020","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[{"id":71856,"name":"Epigenetics \u0026 Genomics"}],"tags":[],"updatedAt":"2020-06-29T19:16:16+00:00","versionOfRecord":{"articleIdentity":"rs-9943","link":"https://doi.org/10.1186/s12864-020-06853-y","journal":{"identity":"bmc-genomics","isVorOnly":false,"title":"BMC Genomics"},"publishedOn":"2020-06-25 12:00:00","publishedOnDateReadable":"June 25th, 2020"},"versionCreatedAt":"2020-05-26 16:07:10","video":"","vorDoi":"10.1186/s12864-020-06853-y","vorDoiUrl":"https://doi.org/10.1186/s12864-020-06853-y","workflowStages":[]},"version":"v3","identity":"rs-9943","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"identity":"rs-9943","version":["v3"]},"buildId":"7rjqhiLT3MXkJMwkYKINL","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}