The complexity and robustness of gene network inside the story between miRNAs and splicing isoforms in sacred lotus

BACKGROUND. Gene expression is complex and regulated by multiple molecular mechanisms, such as miRNA-mediated gene inhibition and alternative-splicing of pre-mRNAs. However, coordination of interaction between miRNAs with different splicing isoforms, and the role 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 expression analyses among six lotus tissues, the negative regulatory roles of most miRNAs are reflected by their tissue-biased expression and the negative correlation with their targets in expression. Further, the central roles of miRNAs in the gene network was unveiled as there are more frequent interactions between miRNAs and hub isoforms than between miRNAs and non-hub isoforms. Surprisingly, 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. The gene function enrichment analysis of miRNA targeted reveals that miRNAs are involved in regulation of lotus growth and development by regulating plant hormone-related pathway genes. CONCLUSION. Taken together, we carry out a comprehensive and deep analysis between miRNA and mRNA transcriptome to study coordination of interaction between miRNAs with different splicing isoforms. Our study on lotus highlights not only the complicate interactions between the miRNAs and transcript isoforms but also functional divergence of many transcript isoforms from the same locus in plant.

there are more frequent interactions between miRNAs and hub isoforms than between miRNAs and non-hub isoforms. Surprisingly, 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. The gene function enrichment analysis of miRNA targeted reveals that miRNAs are involved in regulation of lotus growth and development by regulating plant hormone-related pathway genes.
CONCLUSION. Taken together, we carry out a comprehensive and deep analysis between miRNA and mRNA transcriptome to study coordination of interaction between miRNAs with different splicing isoforms. Our study on lotus highlights not only the complicate interactions between the miRNAs and transcript isoforms but also functional divergence of many transcript isoforms from the same locus in plant.

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) is one of the most important group that can interact with gene in the RNA level. In plant, microRNAs (miRNAs) are 3 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 transcripted 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, a lot of researches were carried out to explore the function of plant miRNAs. Intriguingly, a growing number of studies indicate that the plant miRNAs play key roles in response to plant development, abiotic and biotic stresses through regulating their target genes [4][5][6].
The silencing or translational repression of miRNA target genes has emerged as a post-transcriptional mechanism to regulate the expression of genes [7]. Studies have confirmed that a substantial amount of the miRNAs targets are transcription factors or stress-response factors, involving in important biological processes [8][9][10]. Lacking miRNA regulation, plants would face multiple developmental defects in many critical periods, indicating that miRNAs are a group of the key regulators in plant development. High through-put small RNA sequencing has been proven to be efficient and unbiased to elucidate miRNAs expression profiles and has been employed by many plant studies to uncover miRNAs roles in tissue growth and responding to the environmental factors [11][12][13][14]. Through differential expression analyses, studies found many differentially expressed miRNAs which participate in different processes and pathways such as auxin signal transduction during pollination of maize silks [15], auxin signaling during root development in Arabidopsis [16,17].
RNA alternative splicing (AS) is another important post-transcriptional regulation mechanism, resulting in diverse transcript isoforms encoded by the same gene [18]. With the widespread application of full-length transcriptome sequencing technology, comprehensive isoforms produced by alternative splicing events were identified in many plants [19][20][21]. The structure variation in isoforms can 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 determine the possibility of its silencing by a complementary miRNA, allowing some isoforms to escape from being targeted due to lacing the miRNA binding site. This phenomenon of miRNA escaping through mRNA splicing has 4 been identified in cotton and maize, indicating the gene regulation which can be interplayed by both miRNAs and AS [23,24]. Nowadays, the investigations on the regulated network of miRNA-mRNA showing tissue-specific regulatory manner are identified [28]. However, the interactions between miRNAs and targets at the isoform level and the effects of miRNAs on target gene and isoforms expression profiles are still unclear. In this study, comparative analyses of expression profiles between miRNAs and their target genes were carried out aiming to explore the spatial and temporal regulation of miRNAs in lotus. Combining the identified full-length isoforms dataset and small RNA-seq dataset, we also comprehensively analyzed the interactions between miRNAs and their target isoforms by WGCNA (weighted gene co-expression network analysis) to understand the impact of miRNAs on the expression and function of their target isoforms.

Identification of microRNAs on 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 5 filtered clean reads mapped to the miRBase is 0.33%, i.e. a total of 50 866 reads are aligned to the reference genome (Table 1). After merging with previous lotus miRNAs[27] and removing the redundant (overlapping) MIRNAs (miRNA genes), a total of 1103 potential mature miRNA and 104 miRNA-star sequences were identified, and these miRNAs are produced by 1416 pre-miRNAs ( Fig. 1a) (Additional file 2: Table S1 Fig.S1).
Furthermore, the summary of miRNA nucleotide bias at each position showed there is no nucleotide bias in lotus miRNA nucleotide content (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 first nucleotide bias analysis in other species, we found the bias tendency in 21bp, 22bp and 24bp miRNAs is similar to Camellia japonica [33], pomegranate [34].

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 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 1003 differentially expressed miRNAs were identified in the differential expression analysis. When comparing the miRNA expression level of pollinated carpel with all the other samples, the up-regulated miRNAs outnumbers the down-regulated miRNA, indicating that there is intensive activation of miRNAs in carpel after pollination (Fig.1c).
To explore the miRNA target genes expression pattern among different tissues, pairwise comparisons of these six samples was conducted to identify differentially expressed genes (DEGs). A total of 28 701 DEGs were identified using 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. 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 5.07% (18) (Fig.2a). The proportion of DEGs in DEMTGs is generally higher than that of DEGs in all genes in most between-tissue comparisons, especially in the comparison between carpel and leaf, carpel and petiole (Χ2 test, pvalue<0.01), except for the comparison between petiole and leaf. This, indicates that the differentially expressed miRNAs among tissues have influence on targeted gene expression in some extent.
To explore how intensive the expression pattern of target genes was influenced by the miRNA, the correlation analyses between expression level of the target gene and miRNA across different tissue samples were carried out (Additional file 2: Table S3). In this study, the correlation coefficient (Cor) is

Differentially expressed miRNA and their target isoforms
Taking advantage of transcript isoform analyses from our previous study [28], we analyzed the miRNAtargeted isoforms instead of genes. A total of 10345 unique target isoforms were found (Additional file 2: Table S4). Most target isoforms (10586, 86.42%) are only regulated by one miRNA, while 904 (7.38%) isoforms are regulated by two miRNAs and the rest are affected by more than two miRNAs (Additional file 1: Fig.S5a). Intriguingly the isoforms 'Nn8g40904.1' and 'Nn8g40902.1' can be bound by the many miRNAs, with 38 and 31 homologous miRNAs from the family miR169, respectively. We also calculated the number of regulatory miRNAs 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.S5b). 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.S6), whereas there are 2449 target with only a portion of their isoforms being targeted, such as 'Nn3g21564' gene (Additional file 1: Fig.S6). 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 (Additional file 1: Fig.S7).
The most miRNA target regions in gene body are on coding regions (CDSs) (74.76%), whereas the 5'-8 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 genes expressions. We found that 43.57% of TE-related miRNAs have target gene, and 50.28% of non-TE-related miRNAs have target gene, suggesting that the TE-related miRNAs also play an important role in regulating genes.
To understand the biological functions of miRNAs, especially tissue-specific miRNAs, gene ontology (GO) functional annotation was used. We found that only 1979 out of 4086 miRNA target genes were annotated which their function based on different GO categories (biological processes, molecular functions and cellular components) are provided in supplementary files (Additional file 2: Table S5; Additional file 1: Fig.S8). The most significantly enriched GO terms were "endonuclease activity", "regulation of transcription, DNA-templated" and "Cul4-RING ubiquitin ligase complex", indicating the genes targeted by miRNA regulated numerous key processes and many belonging to transcription factors [35,36]. The specific miRNA may regulate specific gene in different developmental stage, therefore, GO functional enrichment analysis was conducted for six samples (Additional file 1: Fig.S9).
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 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 revealed that the functions of the miRNA target genes are closely related to the tissuespecification.

Functional differentiation of isoforms in the co-expression networks
It is assumed that the closely located gene in the co-expression network, are participating in the same biological process. To investigate lotus isoforms divergence in functions, especially those are miRNA regulated targets, we performed WCGNA at transcript isoform level. We found that some isoforms are exhibiting dramatic expression difference among different tissues. To explore the 9 potential function of miRNA-targeted isoforms in different tissue, we first performed a hierarchical clustering analysis of total isoforms, and found that a substantial portion of isoforms showed strong tissue specificity (Additional file 1: Fig.S10). After filtering out the low 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 the WGCNA. A total 9 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 coexpression 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 identified (Additional file 1: Fig.S11). 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 homology of differentiated floral organs (Additional file 1: Fig.S12). 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 isoforms module allocation, we calculated the content of miRNAs targeted isoforms and number of hub isoform in every module. The proportion of isoforms in all modules being targeted by miRNAs (1785/29391, 6.07%) is significantly lower than the corresponding proportion of isoforms in hubs (121/1500, 8.07%) (Χ 2 test, p <0.01) (Additional file 1: Fig.S13). This suggests that miRNAs preferentially target hub isoforms and play a central role in a gene network.
The isoforms from the same gene are often translated into a protein with different structures and, hence, different functions [22]. To understand the scale of functional differentiation among isoforms from the same gene, we identified isoforms which were assigned to different modules in the coexpression network. Interestingly, among 11,302 genes with multiple isoforms being assigned to modules, 3029 genes have their isoforms divided into different modules (GIDDM). Moreover, 464 of GIDDM were targeted by miRNAs. This supports that a substantial of 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.S14). 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 isoforms of the same gene have functionally diverged, we assume that these different isoforms might convert into different genes (duplicates) to play their independent functions during 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, 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). These results showed that category II is prevalent in this isoform-homology analysis (Additional file 2: Table S7, Additional file 1: Fig.S15). However, interestingly, when we only included the GIDDM, and the proportion of isoforms in category III and especially category I, which represent a stronger degree of sequence structural differentiation between different isoforms from the same lotus gene, were largely increased. This further substantiates that different isoforms from the same gene with functional divergence in the gene network also tend to have higher structural differentiation between each other and more likely to evolve into different duplicate copies.

MiRNA targeted isoforms in plant hormone signaling pathways
To further elucidate the functions of miRNAs and their target isoforms, we investigate the hormones pathways enrichment since they are essential in almost all biological processes in plant. First, KEGG annotation was carried out for miRNA target genes in lotus. 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-, cytokinine-, gibberellin-, abscisic acid-, ethylene, brassinosteroid-and jasmonic acid-associated signaling pathways targeted by 24 miRNAs (Additional file 1: Fig.S16). 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.S17).
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 petiole. In addition to the auxin signaling genes, similar regulatory relationship and expression pattern in abscisic acid signaling were also observed (Additional file 1: Fig.S18). For example, the phosphatase 2C (PP2C) is targeted by five miRNAs, which are lowly expressed in unpollinated carpel. Isoforms from one of two miRNA-targeted genes, "Nn6g35319", are highly expressed in unpollinated carpel, whose five isoforms were clustered into 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 relationship is much complex, suggesting miRNA and gene splicing might play important roles in lotus hormone signaling.

Discussion
The gene regulation model is complicated by our increasing knowledge about non-coding RNAs and transcript splicing. Increasing pieces of evidence indicate that miRNAs play a vital role in plant growth and development through regulating their target genes [4]. Previously, a great number of miRNAs in six tissues have been identified in lotus, and the evolutionary patterns for miRNA families with different age are also unveiled in that study [27]. Also, a great number of transcript isoforms are found in lotus thanks to full length transcript sequencing by combination of PacBio and Illumina [28].
The strategy of combination of miRNAs and mRNA by deep sequencing has been successfully applied in many plant species, such as soybean [37], peanut [38] and cotton [39]. Yet, how the expression of miRNAs in different lotus tissues influences the expression dynamics of their target genes and especially target isoforms across these tissues is largely unknown. Currently, based on an updated miRbase database and new lotus genome assembly, we identified 1 207 non-duplicated mature miRNAs in these lotus tissues, which led us to discover some novel miRNAs which were missed by the previous study. This seems that most of novel identified miRNAs in our current study, are related to transposable elements. Several studies found that many miRNAs, including TamiR1123 that functions in vernalization in wheat, 10 miRNAs in Arabidopsis and 38 in rice are derived from miniature inverted-repeat transposable elements (MITEs) [29, 40,41]. 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 differences in interactions between the miRNA and different transcript isoforms. By building a weight isoform coexpression network, we evaluated the impact of miRNA on isoform expression pattern and uncovered 13 the functional divergence of many isoforms originating from the same gene by partitioning into different functional modules. These results can facilitate our current understanding on gene regulation robustness and complexity in plant.
In this study, most miRNAs are found to be preferentially accumulated in specific lotus tissues, suggesting that most miRNAs have more specialized roles (functions) 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 miRNAs [11,42]. In our study, the tissue-specific miRNAs are preferentially enriched in the reproductive tissues (organs), such as carpel and anther. The target genes functional enrichment showed that the tissue-specific miRNAs were important regulators in anther development, such as formation of spindle fiber during mitosis. For example, the miR393 was identified to regulate the MPS1 gene which regulated cell cycle function during anther development in cotton [43], and the homologous miRNA "nnu-miR393a" was also found to be specifically expressed in lotus anther. Besides, the miR172 family members, which regulate the AP2, are involved in forming the sepal and petal primordia in Arabidopsis [44]. Several miR172 family members, such as "nnu-miR172b", were specifically expressed in lotus anther to probably keep it from being 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 were more active in tissues that undergo dramatic physiological changes, i.e. the pollen-pistil interaction. In addition, the generally higher enrichment of DEGs in DEMTGs than the genome-wide background suggested the differentially expressed miRNAs can influence the expression pattern their target isoforms across different tissue samples.
One novelty of our study is to investigate the relationship between miRNAs and their target isoforms because there is few study focus on this relationship. Studies in human found the interaction between miRNA and hub gene in a gene network [45,46], yet not at isoforms levels. Therefore, the intensive interactions between miRNAs and hub genes from co-expression network in our current lotus study highlights the core regulatory roles of miRNAs. More interestingly, due to the different region 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 to further explore the functional relevance of target genes in different tissues and functional divergence of isoforms using the WGCNA given the increased popularity and success of this approach in many gene networks and functional genomic studies [37,47,48]. As we found that some isoforms from the same gene not only diverge in function by partitioning into different co-expression modules but also exhibit difference of presence and absence in their upstream regulatory miRNAs which is with few reports in plant. Our further homology-based analysis on functional divergence of those isoforms from the same gene also supported that isoforms from the same gene with different expression preference also tend to have their homologs being different genes in another plant species, indicating that gene duplication and splicing isoform might be interchangeable during evolution. This phenomenon of 'isoform-duplicate conversion' is also found during vertebrate evolution including human [49].
Given that a number of studies suggest that miRNA have emerged as key regulators of plant hormone response pathways by affecting their metabolism, distribution, and perception [16], we also focuses on miRNA-targeted isoform relationships in plant hormones signaling pathway. It was well identified that the miR393 regulates the lateral root development in Arabidopsis and leaf morphogenesis in cucumber by down-regulating the receptor TIR1 during auxin signaling [50,51], and the similar miRNA-target relationship was found in our result, 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 [52].
Additionally, the PPC2, abscisic acid signal pathway genes, was found to be targeted by miR166 and involved in regulating the plant height in Gossypium hirsutum [53]. Our result found that PPC2 ortholog gene was targeted by five miRNAs, and the target isoforms were low expressed in petiole, revealing that the miRNA-PPC2 might participate in regulation of lotus plant height. As plant hormones have pivotal role in plant development, these hormone-related miRNAs found in our study might provide the important genetic basis for future genetic studies and manipulation. 15

Conclusion
In the present study, we conducted systematic miRNA and mRNA transcriptome analyses in six lotus tissues and discover novel miRNAs missed by pervious study. The results showed that most miRNAs have tissue-specifically expressed characterization and regulated negatively with their targeted gene.
We constructed a co-expression network in isoforms level by using the WGCNA, highlighting the core regulatory role of miRNAs because of the intensive interactions between miRNAs and hub isoforms by partitioning into different co-expression modules but also exhibit difference of presence and absence in their upstream regulatory miRNAs. Finally, the complicate interactions between the miRNAs and isoforms and functional divergence of isoforms from the same gene can facilitate our current understanding on gene regulation robustness and complexity in plant. sample, the clean reads without adaptors were filtered to retain sRNAs reads length from 18bp to 30bp. Prior to miRNA identification, the retained reads were further processed using miRDeep2

Methods
package [54]. After quality filtering, Bowtie tool was used to map the clean reads to the reference genome (from nelumbo.biocloud.net) with zero mismatch [55]. The mapping result was guided by the known precursor miRNA dataset (miRbase 22.0) in order to identify potential miRNA precursors in the lotus genome while keeping the duplicate loci less than five The hairpin structures and the mapping information was processed as described in the miRDeep-P package [56]. The overlapping precursors were removed and manually checked. To quantify the expression level of miRNA, the count data were normalized to TPM values [57]. Differential expression analysis was performed using edgeR package by controlling the Benjamín-Hochberg false discovery rate (FDR) at 0.05 level (FDR < 0.05) and miRNAs with log fold_change above 1 considered as differentially expressed between two samples.

Isoform identification and quantification
To identify the isoforms, the PacBio full-length transcripts dataset was downloaded and was mapped to the reference as previous described [58]. The petiole, leaf, petal, anther, unpollinated carpel and pollinated carpel Illumina RNA-Seq datasets were also downloaded for quantification of isoform expression. In brief, the clean reads were mapped to the reference by HISTA2 v2.1.0, and then the FPKM (fragments per kilobase per million) value in the genes and isoforms level were calculated by StringTie v1.3.5 using the combination of Illumina and full-length transcripts based annotations. The differentially expressed genes and isoforms are identified by edgeR package using the same threshold used for miRNA differential expression analysis.

Prediction of miRNA target genes
Comparing with the animal miRNA, the most miRNAs and their target mRNAs in plant have nearperfect sequence complementarity [59]. The miRNAs target isoforms were identified by TargetFinder [60], with a high percentage of true positive results as compared to other algorithms [61].
Briefly, FASTA searches was used 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, bugles, 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 used at default. Finally, only predicted targets with scores of three or less were considered.

Gene Ontology enrichment, KEGG pathway analyses and ortholog identification
The GOseq package based on Wallenius non-central hyper-geometric distribution [63] was used to identify the significantly enriched GO terms and KEGG pathway for each gene set, including miRNA target genes and the differentially expressed genes. The KEGG pathway analyswas performed using target gene IDs as queries. The orthologs among Arabidopsis, rice and lotus were identified using the OrthoMCL with an e-value < 1e-15 and an inflation parameter of 2.0.

WCGNA for gene isoforms
To determine whether isoform expression is correlated with the tissues, the co-expression networks were constructed using the WGCNA (v1.0) package in R [64]. To filter non-expressed and invariant isoforms in expression, those with an averaged FPKM > 0.1 and coefficient of variation (C.V.) of FPKM > 2 were retained for the WGCNA. According to Topological Overlap Matrix measure [65], the transcripts were clustered first hierarchically. The transcripts were assigned to 9 co-expression modules using a bottom-up algorithm, known as the dynamic hybrid cut method, named after randomly assigned colors. Modules were made under the parameter that min module size was 600, and the threshold value of the module similarity was 0.5. The isoform kME, a measure of the correlation between isoforms and module, was calculated for each isoform in order to find the relationship between isoforms and modules.

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

Availability of data and materials
All data that support the findings of this study are included in the manuscript and supplementary information.

Author information
to the secondary metabolism pathways in the tea plant (Camellia sinensis).   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).

Scientific reports
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.
33 Figure 5 Correlation analysis of expression between miRNAs and auxin-associated ortholog genes in lotus.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.
Supplementary Tables.xls Supplement Figure.pdf