Comprehensive analysis of small RNA-seq data reveals that combination of miRNA with its isomiRs increase the accuracy of target prediction in Arabidopsis thaliana

Identification of targets of miRNAs and their corresponding isomiRs

We became particularly interested in how miRNA terminal heterogeneity impacts its binding to target transcripts. Imprecise cleavage by DCL generates heterogeneous termini that can alter the miRNA seed region (2–8 nt from the 5′-end) and lead to different mRNA targets. However, to the best of our knowledge, no systematic study has been conducted yet to determine the effect of miRNA terminal heterogeneity on target mRNA selection. To understand the biological role of isomiRs, we predicted isomiR targets using psRNATarget.³⁸ miRNAs, and thus their isomiRs, were selected for this study if (1) experimental validation of target mRNAs was published, and (2) at least 8 known isomiRs are produced from the miRNA (Tables 1A and 1B). The isomiRs were mapped to the 5′-end of the canonical miRNA, with at most a 5-nucleotide shift in the upstream or downstream direction.

All predicted targets were further analyzed to understand the varied patterns that arise due to terminal heterogeneity of the miRNA. To accomplish this, a binary matrix was created in which each row indicated the predicted mRNA transcript (TAIR10, http://www.arabidopsis.org/) and each column specified an isomiR sequence. Values of “1” or “0” signified the presence or absence, respectively, of an isomiR target site within the transcript. The isomiRs were arranged based on their 5′-end position and length. Gitools (version 1.6.1, http://www.gitools.org/) was then used to generate a heat map of the matrix for easier visualization. Moreover, we clustered mRNAs based on the number of predicted isomiRs that may target them.

Intriguingly, we observed that most of the validated mRNAs could be targeted by the majority of isomiRs and the miRNA, indicating that “authentic” targets have binding sites for multiple isomiRs and not only by the “canonical” miRNA (Fig. 2, S1, and S2). We also identified several isomiRs with specificity toward distinct mRNAs. Figure S1 shows that ath-miR156a has a total of 17 isomiRs (including the miRNA) targeting 47 mRNA transcripts. However, all of the experimentally validated mRNAs are targeted by 15 or more (88–100%) isomiRs. In the case of ath-miR164a, 11 isomiRs were found to target 69 mRNAs. Each of the experimentally validated mRNAs is targeted by all of the isomiRs (Figure S2). Similarly, analysis of ath-miR172a produced 20 isomiRs that could target 90 mRNAs (Fig. 2). Ten of 14 experimentally validated target mRNAs were predicted to be recognized by more than 18 isomiRs (90–100%). Our analysis predicted that one mRNA (i.e., AT3G14770.1) could be targeted by 10 isomiRs, while 3 mRNAs (i.e., AT1G53780.1, AT1G53780.2, and AT1G53780.3) are targeted by one isomiR.

To understand the significance of isomiRs in target selection, we calculated the percentage of isomiRs that can target experimentally validated mRNAs. Here we only considered targets predicted by psRNATarget (highlighted in green in Tables 1A and 1B). This analysis revealed 31 miRNAs (and their isomiRs) capable of targeting 2,163 mRNAs, of which 2,075 are unique and 88 are targeted by multiple miRNAs. Table 2 demonstrates that if we increase the percentage of isomiRs targeting an mRNA, there is an increased likelihood of the mRNA being an experimentally validated target. The mRNA targeted by miRNA and its all isomiRs can cover ∼57% of validated targets (Table 2). This result indicates that “authentic” target mRNAs are recognized by multiple isomiRs, and not only by the canonical miRNA. Furthermore, we analyzed the recognition sites within a few transcripts and found that isomiRs can target the same site as the canonical miRNA with slight variations in a few nucleotides (Fig. 3A, 3B, and S3). Therefore, the possibility exists that isomiRs can sort into different AGO proteins to target the same transcript through different RNA interference (RNAi) mechanisms. Our computational models support this hypothesis (Tables 3 and S2). In addition, we checked the pattern of miRNA:target binding in humans. Like Arabidopsis, we observed that experimentally verified target mRNAs can be recognized by several human miRNA/isomiRs (see Supplemental Results and Discussion).

Gene ontology analysis of targets

Transcripts targeted by Arabidopsis isomiRs were further subjected to Gene Ontology (GO) analysis using the agriGO tool (http://bioinfo.cau.edu.cn/agriGO/).³⁹ We hypothesized that mRNAs targeted by a single miRNA and its isomiRs are natural targets involved in a common, specific biological function. Our results show that mRNAs targeted by more than 50% of isomiRs derived from ath-miR156a are significantly enriched for 23 GO terms, while another group of test transcripts targeted by up to 50% of isomiRs had no significant terms in common (Table 4 and Table S3). We also observed similar patterns of GO term enrichment for other test sets. For example, 5 and 4 GO terms were uncovered for ath-miR164a and ath-miR172a (Table 4 and Table S3), respectively. Moreover, we searched the significance of GO terms in all predicted targets of 31 miRNAs and their isomiRs. This analysis can further help us to understand the significance of isomiRs in target finding, and overall functional characteristics of miRNA target genes. Among 2,163 mRNAs, 1,792 mRNA are targeted by up to 50% of isomiRs, while only 371 mRNA are targeted by more than 50% of isomiRs. Interestingly, we found significantly enriched 96 GO terms in a group of mRNA targeted by more than 50% of isomiRs, while other genes targeted by up to 50% of isomiRs does not show significant GO terms (Table 4 and Table S3). Among them, the terms “transcription factor activity” and “transcription regulator activity” appeared regularly. This result further supports our observation that authentic miRNA target genes are also targeted by various isomiRs instead of by only the canonical miRNA, and that they most likely participate in transcription factor activity. Figure 4 shows interactive GO graphs according to biological process. Supplementary Figure S4 depicts results from the GO analysis according to cellular component and molecular function.

In addition, we studied how isomiRs and canonical miRNA influenced the performance of target predictions. Our result showed that targets predicted with only canonical miRNA have significant GO terms, while targets predicted with isomiRs without canonical miRNA does not have significant GO term (Table 4 and Table S3). Our result showed that 31 miRNAs targeting 569 mRNAs in the class of canonical miRNA, and among them 108 mRNAs are experimentally known (Table 4). On the other hand 371 mRNAs targeted by class of more than 50% of isomiRs, and among them 102 mRNAs are experimentally known. Therefore, sensitivity of target prediction increases from 18.98% to 27.49%, and enrichment of GO terms become stronger if isomiRs were combined with the canonical miRNA (Table 4). Though, complete targets of miRNAs were not experimentally known, but if we consider the available list of experimentally known targets, the true positive and true negative targets would be 108 and 1578, respectively for class of canonical miRNA. While true positive and true negative targets would be 102 and 1770, respectively for class of more than 50% of isomiRs. Therefore accuracy of target prediction increase from 77.9% to 86% if isomiRs were combined with the canonical miRNA (Table 4). This further supports that combination of isomiRs with miRNA can predict the target with high accuracy.

Furthermore, we assessed the significance of isomiRs in target prediction using computationally generated all possible isomiRs of ath-miR172a from its precursor sequence with length varies from 18 to 28 nt (see Length diversity of sib-miRs of Supplemental Results and Discussion). Two sets of isomiRs were generated having heterogeneous termini of: (I) -3/+3 nt, and (II) -5/+5 nt, relative to 5′-end of canonical miRNA. These isomiRs along with canonical miRNA were used for target prediction by psRNATarget, and compared with experimentally validated targets (Table 4 and Table S4). Compared to natural isomRs, randomly generated isomiRs can increase the sensitivity from 77% (10/13) to 91% (10/11) and give strong GO enrichment (Table 4). We observed that mRNA targeted by up to 50% of isomiRs does not show any GO terms, while mRNA targeted by more than 50% of isomiRs showed the GO enrichments which is similar to the target of natural isomiRs (Table 4). It was interesting that we observed that targets predicted with only canonical miRNA had prediction sensitivity of 50% (11/22) and targets predicted by more than 50% of isomiRs (including canonical miRNA) had significantly increased prediction sensitivity, which was 91% (10/11) (Table 4).

Evidence of target cleavage by isomiRs

Concentration level of isomiRs are considered to be lower compared to canonical miRNA, but still can cleave the target gene and have a biological impact on gene regulation. In order to prove this hypothesis, we investigated the isomiRs directed target cleavage through publically available degradome data of polyadenylated transcripts using a computational pipeline CleaveLand.⁴⁰ Experimental studies showed that miRNA cleaves the target mRNA precisely between the 10th and 11th nucleotides relative to the 5′-end of miRNA.^40-42 The resulting downstream fragment of the cleaved mRNA is uncapped and polyadenylated, called as degradome,⁴² and has been exploited to identifying miRNA guided cleavage site accurately in the mRNA transcripts ^40,43-45 Through degradome data analysis, we have detected the cleaved targets by both isomiRs and canonical miRNAs (Figure S5A–L). Our analysis also showed that some isomiRs are more effective in target cleavage than their canonical miRNA (Figure S5P–Q). Interestingly, we also detected the targets which are cleaved by isomiRs but not by its canonical miRNA (Figure S5M-N), and thus support the finding of other study showing significance of isomiRs in gene regulation.⁴⁶

Prediction of novel targets using miRNA and isomiRs

We observed that a “true” target is suppressed by various isomiRs and the canonical miRNA in our genome-wide bioinformatics analysis. Furthermore, GO analysis showed that target mRNAs of these isomiRs are enriched for a particular biological function, which indicates that natural miRNA targets may be recognized by several isomiRs. Based on these observations, we extended and applied our new strategy to predict novel miRNA targets. In this analysis, a transcript was considered a probable target if it contained a target site for 80% or more isomiRs (including the miRNA). We also lowered the threshold if targets were not found using this criterion. Therefore, the highest number of isomiRs associated with a target indicates a real target gene. Tables 5A and 5B list the miRNA targets predicted with this method.

Biological validation of predicted targets

The predicted targets of ath-miR158a were validated experimentally to confirm the accuracy of our method. Our analysis predicted 4 targets, namely, AT1G49910, AT1G64100, AT3G03580, and AT1G62860 (Table 5B), for ath-miR158a and its isomiRs. Therefore, we hypothesized that, in Arabidopsis, the absence of ath-miR158a and its isomiRs leads to increased expression of these 4 target genes. To test this hypothesis, we obtained 2 T-DNA knockout mutants for the ath-miR158a gene (Figure S6), namely SALK_025691 and SALK_083227, from ABRC (https://abrc.osu.edu/). The expression levels of miR158 and its isomiRs were quantified by RT-qPCR in plants that were homozygous for the respective T-DNA insertions. RT-qPCR was performed using forward primers against miR158 and an isomiR predicted to be generated toward the 5′-end of miR158. Wild-type (Col-0) plants showed amplification of both miR158 and isomiR158, indicating the expression of both sRNAs in Arabidopsis (Fig. 5). Our data also indicated that the expressions of the mature miR158a and one of its isomiRs were significantly reduced in both mutant plants (Fig. 5).

Endogenous levels of AT1G49910, AT1G64100, AT3G03580, and AT1G62860 mRNA were also assessed by RT-qPCR in the leaf tissues of mutant and wild-type plants. The SALK_025691 mutant plant exhibited an approximately 4-fold increase in AT3G03580 mRNA compared to wild-type (Fig. 6A). This mutant also showed significantly higher levels of AT1G64100, and AT1G62860 mRNA. Degradome data also support our finding that AT1G62860 mRNA is the target of miR158a (Figure S5O). Similarly, SALK_083227 mutant plants showed a 3-fold increase in AT3G03580 mRNA level (Fig. 6A). This increase in expression of 3 out of 4 predicted targets in ath-miR158a-deficient plants demonstrated that these mRNAs are regulated by ath-miR158a in wild-type plants. These mutants also showed significantly decreased level of expression of miR158-derived isomiRs (Fig. 5). However, not all predicted targets were confirmed. For example, AT1G49910 mRNA was expressed similarly in wild-type and mutant plant leaves (Fig. 6A). Also, unlike SALK_025691 mutants, SALK_083227 mutant plant leaves did not display an increase in AT1G49910, AT1G64100, and AT1G62860 mRNA levels.

To understand the reason for these variations, we investigated the expression pattern of these 4 mRNAs in different plant parts across different growth stages using Arabidopsis eFP Browser, a publicly available database (http://bar.utoronto.ca/efp/cgi-bin/efpWeb.cgi).⁴⁷ Data from the Arabidopsis eFP Browser indicates that all 4 mRNAs are highly expressed in reproductive tissues or meristems, especially flower tissues (Figure S7), but not in leaves. Thus, we further analyzed the endogenous expression level of all 4 mRNAs in the flowers of mutant and wild-type plants by (Col-0) RT-qPCR. Results from this experiment demonstrated that all 4 transcripts are significantly up-regulated in the mutant flower tissues compared to wild-type (Fig. 6B). Convincingly, AT1G62860 mRNA was highly expressed in the flower tissues of both mutant plants.

Through bioinformatics analysis, we also predicted that ath-miR869 can target AT4G08210 mRNA (Table 5B). To verify this, we analyzed the expression of AT4G08210 mRNA in an Arabidopsis mutant, SALK_052117, which contains a T-DNA insertion in the ath-miR869 gene (i.e., AT5G39693) by RT-qPCR. Our results showed that this mRNA was expressed more than 10-fold higher in mutant plant leaves compared to wild-type (Figure S8A). Similarly, expression of the ath-miR5028 target, namely AT1G68540, was analyzed (Table 5A) in SALK_025566 and SALK_108642 mutants, which possess a mutation in the ath-miR5028 gene. Consistent with the previous results, AT1G68540 mRNA was significantly increased in both the SALK mutants when compared to wild-type plant leaves (Figure S8B). These experimental results validated the predicted target genes, and thus demonstrated that our bioinformatics approach of combining the sequences of canonical miRNA and its isomiRs can accurately predict the target of miRNA in Arabidopsis.

Data sets

We retrieved sRNAs of Arabidopsis from Plant MPSS databases (http://mpss.udel.edu/) released on November 25, 2011.³⁷ The database contains information about sRNAs that are perfectly mapped to 231 Arabidopsis pre-miRNAs taken from 51 sRNA high-throughput Sequencing-by-Synthesis (SBS) sequencing libraries (104,576,308 sequences in total) covering a diversity of tissues and cell lines. We retrieved sequences mapped to each miRNA precursor with their normalized read number (read per million; RPM) from various libraries along with their genomic coordinates (TAIR10).

According to miRBase 18 (release November 2011), Arabidopsis possesses 291 miRNA precursors that generate 328 mature sequences. Among them, 200 were placed into 97 different families based upon their sequence similarity (miRFam.dat of miRBase). To achieve more accurate sequence analysis, a clean and non-redundant dataset was prepared according to 2 criteria: (A) taking one precursor from each family, and (B) taking all precursors not assigned into any family. Among the 231 precursors obtained from the Massively Parallel Signature Sequencing (MPSS) database, 90 belonged to category (A) while 55 were from category (B). Moreover, we removed sRNAs whose total read number was less than 10 RPM in all 51 sRNA libraries, and mapped to opposite strands of the genome relative to annotated pre-miRNA. Thus, the final data set (termed pre-miR_data) contained 132 precursor sequences and their perfectly matched sRNAs, which we called sib-miRs. The number of sib-miRs and their expression level across 51 sRNA libraries was given in the Table S1.

Furthermore, we constructed a dataset of 132 precursors and the number of species in which these precursors are conserved. We also included the number of unique sib-miRs and its read counts (RPM). The data set (called as conserve_pre-miR.txt) was used to study the relationship between miRNA gene conservation and their expression.

To study the cleavage accuracy of DCL on pre-miRNA and the effect of terminal heterogeneity of mature mRNA on targeting genes, a dataset of isomiRs, isomiR_data, was created. For this, sib-miRs mapping within 5 nucleotides upstream and downstream relative to end (5′-end or 3′-end) of canonical miRNAs were taken. Here we considered canonical miRNAs as annotated in miRBase 18.

We also compiled data of experimentally validated targets of miRNAs from miRTarBase, SeqTar, and the published literature into a data set called exp_valid_miR_target.txt.^12,43,53,72 We used an in-house-developed Perl script, along with careful manual inspection, to select and analyze the final data.

Data normalization

Datasets were normalized using 2 different methods. The first method involved normalizing the counts as RPM. This method involved calculating the read counts of an sRNA divided by the total number of read counts in that library and multiplying by 10⁶. The data set contains only “genome-matched reads" excluding t/r/sn/snoRNAs. Throughout this manuscript, we have used “read counts” to indicate “normalized read counts” unless otherwise mentioned. The second method involved normalizing within the pre-miRNA. Here, the read count of a sRNA was divided by the total number of sRNA read counts generated from the same pre-miRNA and multiplying by 10³. This calculation normalized the sib-miRs at the single pre-miRNA level. We considered read counts of sib-miRs in RPM, derived from 51 libraries, and mapped to a pre-miRNA.

Sorting into AGO proteins

All unique sib-miRs were first exactly matched with experimentally known sequences with AGO (Gene Expression Omnibus (GEO):GSE10036) and assigned a specificity for each AGO protein.⁷³ The remaining sequences were predicted for AGO specificity using SVM models developed through standard protocols.^9,74,75 These models were developed using features of nucleotide frequency of sRNA and binary patterns of 6 nt from both terminals of sRNA⁹ of AGO-specific experimental sequences (GEO:GSE10036). The binary models of AGO1, AGO2, AGO4, and AGO5 achieved an accuracy of 85%, 78%, 88%, and 93%, respectively, on threshold “0” using 5-fold cross validation techniques (Table S5). The performance of these models as assessed on independent datasets showed high accuracy (Table S5).

Analysis of DCL cutting accuracy

Terminal heterogeneity showed absolute distances between the observed 5′- or 3′-ends of sRNA and the reference miRNA (annotated in miRBase 18). We adopted a previously described method to calculate the terminal heterogeneity of each miRNA.^27,28 In brief, to calculate 5′-end heterogeneity, sequences mapped within 5 nucleotides upstream or downstream of the 5′-end of a mature miRNA was used. A heterogeneity value was calculated for each canonical miRNA using the total shift in cleavage site divided by the total copy of sequences.

For example, the 6 sRNA sequences mapped to the 5′-arm of ath-miR164a as follows:

1. ath-miR164a|0|21nt|50

2. ath-miR164a|-1|22nt|25

3. ath-miR164a|1|21nt|20

4. ath-miR164a|2|19nt|13

5. ath-miR164a|-3|23nt|7

6. ath-miR164a|3|22nt|5

where the format of ath-miR164a|-1|22nt|25: miRNA, -1 is (1 nt upstream) the shifted nucleotide position in the sequence compared to the 5′-end of the canonical ath-miR164a, 22 nt is the sequence length, and 25 is the read count in RPM. Total copy of sequences = 50+25+20+13+7+5 = 120 Total shift in cleavage site = 50*abs(0)+25*abs(-1)+20*abs(1)+13*abs(2)+7*abs(3) +5*abs(3) = 107 Heterogeneity value for ath-miR164a is : Total shift in cleavage site / Total copy of sequences =107/120= 0.89.

Finally, the mean heterogeneity of the 5′-end terminal was calculated as the average heterogeneity of all studied canonical mature miRNA.

miRNA/isomiRs and their target identification

The computational tool, psRNATarget (http://plantgrn.noble.org/psRNATarget/), was utilized to identify isomiR target mRNAs.³⁸ We used TAIR10 as a genomic library with default parameters except the hpsize, which we set to 15. We named each sequence of isomiRs in a format which gives complete information about the sequence relative to canonical miRNA. For instance, “ath-miR164a_-1_21_21_0.11359” indicates ath-miR164a as the canonical miRNA, -1 represents shift in one nucleotide upstream compared to the 5′-end of canonical ath-miR164a, 21 is the length of the isomiR, 21 is the length of canonical the ath-miR164a, 0.11359 is the normalized value of the isomiR within the pre-miRNA. Therefore, “ath-miR164a_0_21_21_985.78476” would be the name of canonical miRNA expressed in great abundance.

Gene Ontology (GO) analysis

The Singular Enrichment Analysis (SEA) of agriGO (v 1.2) tool was used for GO term analysis.³⁹ To accomplish this, we divided the predicted target transcripts of a miRNA and its isomiRs into 2 groups: (1) mRNA targeted by up to 50% of all isomiRs, and (2) mRNA targeted by more than 50% of all isomiRs. We then analyzed both groups for Gene Ontology enrichment by comparing all targets as “Customized reference." The Fisher statistical test, Bonferroni correction, a significance level of 0.01, and a minimum number of mapping entries equal to 5 were selected as additional parameters.

Degradome data analysis

We downloaded 7 libraries of degradome data of Arabidopsis from Plant MPSS databases (https://mpss.udel.edu/at_pare/): (1) AxIDT [GSM278334], (2) AxIRP [GSM278335], (3) AxSRP [GSM278370], (4) Col [GSM284751], (5) ein5 [GSM284752], (6) TWF [GSM280226], and (7) Tx4F [GSM2802267]. These files were already normalized, where information of each sequence was given in the form of sequence and its normalized copy number value. In order to use in CleaveLand, each sequence was represented with multiple copy of sequence according to its copy number, and all files were merged into a single file. Sliced targets of isomiRs and miRNA were identified and classified using CleaveLand 4.3 pipeline. For analysis, we used 3 different input files: (1) degradome data, (2) Our's list of isomiRs and miRNA sequences, and (3) Transcript database of Arabidopsis previously used in psRNATarget analysis (http://plantgrn.noble.org/psRNATarget/).

Plant material and growth conditions

Arabidopsis T-DNA insertion mutants of the ath-miR158a, ath-miR869, and ath-miR5028 genes were obtained from the Arabidopsis Biological Resource Center in Columbus, Ohio, USA. Arabidopsis ecotype Col-0 was used as the wild-type control. All mutant plants used in this study were confirmed to be homozygous for their respective T-DNA insertion (described below). Detailed information about the mutant plants and the respective gene IDs used in this study is described in Supplementary Table S6. Wild-type Col-0 and mutant seeds were cold-treated for 3 d at 4°C and germinated in Metro Mix Professional Growing Mix (SUNGRO) in a growth chamber maintained at 20° to 22°C. The plants were grown at this temperature under short day conditions (8 h of light).

Zygosity test for mutant plants

Arabidopsis wild-type Col-0 and respective mutant plants were grown for 3 weeks and leaf samples were collected to test for zygosity using the method recommended by SALK (http://signal.salk.edu/tdnaprimers.2.html). Genomic DNA was isolated by using the REDExtract-N-Amp Plant PCR Kit (Sigma-Aldrich Co. St. Louis MO) according to the manufacturer's recommended protocol. Two paired PCR reactions using primer combinations, namely LP+RP and LB+RP (primer details are provided in Supplementary Table S7), were performed. Reactions that showed amplification of the expected PCR product with the LB+RP primers and no amplification in the LP+RP reaction were considered homozygous for the respective mutation.

Quantification of endogenous miRNAs

The homozygous mutant plants for miR158, namely SALK_025691 and SALK_083227, were confirmed to be deficient in the endogenous mature miRNA transcript using primers as described in Supplementary Table S8. Total RNA was extracted from 5-week-old Arabidopsis plant leaves using the RNeasy plant mini kit (Qiagen Valencia, CA, USA) according to the manufacturer's instructions. The High-specificity miRNA RT-qPCR detection kit and protocol (Agilent Technologies Inc.. Clara, CA USA) were used to quantify miRNAs and isomiRs. This kit has been optimized for the detection and quantification of miRNA from 15–250 ng of total RNA. RT-qPCR (BioRad-CFX connect real-time system, Hercules, CA, USA) was performed using a miRNA-specific forward primer (Supplementary Table S8) and universal reverse primer provided in the kit. The forward primer for isomiR158 was designed using the miRNA primer design tool at http://www.astridbio.com/ and synthesized by Integrated DNA Technologies (IDT, Inc.., Coralville, IA, USA). Forward primer design and RT-qPCR were performed as described previously⁷⁶ with modifications. Relative expression as fold increase was calculated using the ΔCP method⁷⁷ with comparison to a reference (i.e., AtUbiquitin 5 and miR159).

Real-time PCR for quantification of endogenous mRNA expression

Total RNA was extracted from the leaf or floral tissues as described above. RNA samples were treated with RNAse-free DNAse prior to reverse transcription using the Superscript^TM First Strand Synthesis System for RT-PCR (Invitrogen, Life Technologies, Carlsbad, CA, USA). First strand cDNA was synthesized from 2 μg of total RNA. Quantitative real-time PCR (RT-qPCR) was conducted with a SYBR green (Sigma Aldrich Co. St. Louis MO) reaction mix using a 1:20 dilution of Cdna.⁷⁸ Using the Ct values, the relative expression of target transcripts was calculated by the ΔCP method⁷⁷ with comparison to the reference genes AtEF1α, and AtUbiquitin 5. Primers for the respective endogenous transcripts were designed using Primer Quest software (IDT, Inc..) and the primer sequences provided in Supplementary Table S9.