Quantitative Amplification of Cleaved Ends (qACE) to assay miRNA-directed target cleavage

Plant material

Arabidopsis thaliana wild type Ler and hen1-1 mutant seeds (obtained from ABRC, Columbus, OH: Stock#CS6583) were surface-sterilized and grown in Sunshine mix #1 (Tessman Company, Sioux Falls, SD) under 16h light at 25°C. Leaf, stem and flower tissues were harvested at the same developmental stages, immediately frozen in liquid N₂ and stored at -70°C. Soybean seeds Glycine max cv. Williams-82, the genotype used for genome sequencing project²⁰ were surface-sterilized²¹ and germinated in 4” pots filled with a mixture of vermiculite: perlite (Hummert International, MO) in the ratio of 1:3 and watered with nitrogen free plant nutrient solution²². The plants were grown in a controlled environment vertical growth chamber (Conviron Growth chamber, Manitoba, Canada). Growth conditions used were: 16h light and 8h dark and 50% relative humidity with a day and night temperature of 25°C and 20°C respectively. For nodulation assays 5 days old germinated plants were inoculated with Bradyrhizobium japonicum (USDA110) grown in Vincent’s rich medium prepared and supplemented with Chloramphenicol (antibiotic selection marker) at 30°C²³. The plants were inoculated with B. japonicum at a concentration of OD₆₀₀=0.08 and 14 days post inoculation mature nodules and the root segment adjacent to the nodule were harvested separately.

Over expression of gma-miR164a

To over-express miR164, the precursor of gma-miR164 (miRBaseID MI0007209) was PCR amplified using G. max genomic DNA as template. The PCR product was initially cloned in to an entry vector PCR8/GW/TOPO (Invitrogen, Carlsbad, CA) and subsequently in to the binary vector pCAMGFP-CsVMV: GW²⁴ using a Gateway LR clonase (Invitrogen, Carlsbad, CA) reaction. The vector was transformed in to Agrobacterium rhizogenes K599 cells, using electroporation and composite plant transformation was performed as described²⁵. Transgenic GFP roots were collected on dry-ice and stored at -70°C until RNA isolation and cDNA synthesis.

RNA isolation and cDNA synthesis

Total RNA was isolated from the different tissues using TRI reagent (Product#T9424, Sigma Aldrich, St. Louis, MO) as described previously²¹. The RNA was quantified using Nanodrop spectrophotometer (ND1000, Thermo Scientific, Wilmington, DE) and RNA integrity was verified using agarose gel electrophoresis. Adapter ligated cDNA was prepared using GeneRacer kit (Product#L1502-01, Invitrogen, Carlsbad, CA) as per the manufacturer’s instructions except that calf intestinal phosphatase treatment to remove 5’ phosphate and 5’ Cap removal using tobacco acid pyrophosphatase steps were omitted. For Arabidopsis gene expression analysis 2µg of total RNA each of leaf, stem, and flower was pooled together and subject to adapter-ligated cDNA synthesis. For gene expression analysis in soybean nodules and adjacent root tissues, 7µg of total RNA was used in adapter-ligated cDNA synthesis. For miRNA stem-loop cDNA synthesis, 1µg of total RNA was used and a multiplexed cDNA was performed for all miRNAs examined using M-MuLV reverse transcriptase^26,27. miR1515 & U6 was used as normalization/house-keeping control for miRNAs²⁸.

qPCR and qACE assays

qPCR assays were performed in MX3000P thermocycler (Stratagene/Agilent technologies, Santa Clara, CA) using SYBR Advantage qPCR premix (Product#639676, Clontech, Mountain View, CA). The data was analyzed using the MxPro software. Relative expression values for full length transcripts were performed using the dCt method²⁹ using Actin for soybean and U6 for Arabidopsis as normalization control. Relative levels of cleaved transcripts (described as RCT in text) were calculated using the dCt method comparing the Ct values of full-length vs. cleaved transcripts in the same sample. For statistical analyses, we calculated the range of possible expression values in each sample based on the deviation between Ct values of replicates. Error bars indicate this range in each sample.

Principle of qACE

qACE is an extension of the modified RLM-5’RACE method^15,30 used to qualitatively validate target cleavage and the subsequent semi-quantitative method¹⁴ used to examine the levels of cleavage remnants. An RNA adapter is ligated to cleaved 5’-ends of polyA or total RNA preparations (“adapter-ligated RNA”). The adapter gets ligated to any available 5’-phosphate of ribose molecule of RNA arising due to miRNA-directed cleavage or other means (e.g. mRNA degradation), but not to full-length mRNAs, because of the 5-methyl Guanosine/CAP. Oligo-dT primed cDNAs generated from adapter-ligated RNA are subject to real-time PCR where a primer that spans the adapter junction is used as a forward primer (subsequently referred to as qACE forward primer) and a gene-specific reverse primer (Figure 1).

Figure 1. Quantitative Amplification of cDNA Ends (qACE).

Total RNA containing full length transcripts with 7-methyl Guanosine/CAP and cleaved transcripts arising due to miRNA activity is used as starting material. An RNA oligo adapter of known sequence is ligated to the cleaved transcripts (with a 5’-Phosphate group) giving rise to adapter ligated mRNA. Oligo dT-primed cDNA is used to assay the levels of full length mRNA (qPCR) and cleaved transcripts (qACE). The qPCR forward and reverse primers (qPCR-F & qPCR-R) were designed across the miRNA binding site whereas the qACE forward primer (qACE-F) is designed complementary to adapter sequence but with six nucleotides at the 3’end specific to the sequence at 5’end of cleaved transcript, and a gene specific reverse primer (qACE-R).

The forward primer had to be designed such that it does not amplify full-length cDNA molecules or bind non-specifically to other adapter-ligated molecules. We achieved this by careful design such that this primer corresponds to the adapter-cleavage product junction with 6 nucleotides at the 3’end of the primer corresponding to the nucleotides downstream of the cleavage site in the target gene. A shorter design (a 4nt sequence for example), had an increased probability of non-specific binding. On the other hand, longer than 6nt design had a higher thermodynamic stability making it possible to bind to full-length cDNA molecules (depending on the GC content). In addition, for target genes with identical miRNA cleavage signatures (and therefore identical qACE forward primers), specificity was achieved by carefully designing gene-specific reverse primers.

We expected that the qACE forward primer combined with a gene-specific reverse primer would help specifically amplify ligation products where the adapter is ligated to the predicted/validated cleavage site. Therefore, a qPCR assay using these primers will specifically quantify cleaved mRNAs resulting from miRNA-directed cleavage. We also expected that the qACE forward primer will not amplify full-length cDNA molecules and neither would it efficiently amplify ligation products where the cleaved end does not correspond to the miRNA-directed cleavage site (Figure S4a1). Finally, we also expected that the use of a gene-specific reverse primer would amplify linearly and help distinguish different targets of the same miRNA even if the miRNA-binding sites are identical/highly conserved (Figure S4a–Figure S4d). In addition, the same cDNA preparation can be used to quantify full-length molecules of target mRNA as well as using appropriate primers that span across the miRNA binding site (subsequently referred to as full-length qPCR primers). Ratio of cleaved transcripts vs. full-length transcripts will serve as a quantitative indicator of the extent of miRNA-directed cleavage of the target mRNA. The level of cleaved target transcripts also depends on the levels of transcriptional activity of the target gene. Therefore, comparing absolute levels of cleaved transcripts (or levels normalized to house-keeping genes) between two conditions or tissue types might not reveal differences in miRNA regulation. We decided to normalize the levels of cleaved transcripts to full-length transcripts using the dCt method to obtain the ratio of cleaved transcripts (RCT). The formula for calculating RCT is, 2^{(C_{t(Full length)-}C_{t(cleaved transcript)})}, where C_t stands for threshold cycle. It should be noted that qACE relies on prior knowledge about the validity of the cleavage site e.g. obtained through 5’-RACE or degradome/PARE analyses. It cannot be used to validate a predicted cleavage site. In addition, one has to assume that the stability of cleavage products are identical among multiple target genes that are compared. However, such an assumption has to be made even for alternate methods such as degradome/PARE analysis.

qACE demonstrates increased GmNAC1b transcript cleavage in soybean roots over-expressing miR164

To obtain proof of concept for the method, we examined the abundance of cleavage remnants in miR164 over-expressing soybean roots. We isolated total RNA from vector control and miR164 over-expressing roots and examined mature miR164 levels by stem-loop qPCR. We observed an 6-fold increase in miR164 over-expressing roots compared to control roots as expected (Figure 2a)³¹. We also examined the expression of GmNAC1b, a 5’-RACE-validated target of miR164 in soybean (Figure S1). GmNAC1b is a close ortholog of AtNAC1 (Figure S2). A potential cleavage product of AtNAC1 accumulates in the roots of Arabidopsis plants over-expressing miR164³³. To quantify the levels of GmNAC1b, we generated oligo-dT-primed cDNA and performed qPCR assays using primers designed across the miR164 binding site of GmNAC1b. Results from qPCR assays indicated that there was a 4.2-fold reduction in the levels of full-length GmNAC1b transcripts in miR164 over-expressing roots compared to control roots (Figure 2b). These results indicated that indeed over-expression of miR164 resulted in a reduction in full-length transcripts of its target.

Figure 2. Increased cleavage of GmNAC1b transcripts by over expression of miR164 in soybean roots.

(a) The relative expression level of miR164 (normalized to that of miR1515) in control soybean roots and those over-expressing miR164 (miR164ox) determined using stem-loop qPCR. (b) The relative expression levels of full-length GmNAC1b (normalized to that of GmActin) in control and miR164ox roots assayed using qPCR. The data shown in (a) and (b) are average of three replicate assays and error bars indicate the range of possible values based on deviation between Ct values of replicate assays. (c) The Ratio of Cleaved Transcripts (RCT) over full-length transcripts of GmNAC1b in control and miR164ox roots assayed using qACE. The data shown are average of three replicate assays and error bars indicate the range of possible values based on deviation between Ct values of three replicate assays.

Next, we used qACE to detect and quantify target cleavage in these roots. We designed qACE primers for GmNAC1b as outlined in Figure 1. First we examined if the qACE primer pair does not amplify full-length molecules of GmNAC1b using the above cDNA from total RNA without any adapter ligation. qPCR assays detected no amplification (Figure S4a1) indicating that indeed qACE primers did not amplify/detect full-length GmNAC1b molecules. Next, we prepared adapter-ligated cDNA by ligating a known RNA adapter (see methods) to total RNA and reverse-transcribing these molecules using an oligo-dT primer. We used this cDNA to assay full length and cleaved molecules of GmNAC1b using the appropriate primer pairs (Figure 1) and calculated RCT values for each NAC1b in control and miR164 over-expressing roots. There was an 11.6-fold increase in the ratio of cleaved GmNAC1b transcripts in miR164 over-expressing roots compared to the control roots (Figure 2c) indicating that indeed miR164 over-expressing roots had increased cleavage of GmNAC1b. This experiment demonstrated that qACE could detect and provide a quantitative indicator of the levels of target cleavage.

qACE demonstrates reduced target cleavage in miRNA-deficient Arabidopsis hen1-1 mutant plants

Next, we examined the levels of target cleavage in Arabidopsis wild-type (Ler) and miRNA-deficient hen1-1 mutant plants. It was previously demonstrated using Northern blots that hen1-1 mutants accumulated less amounts of miRNAs and that there was an increase in the levels of miRNA targets in these plants³⁴. We used stem-loop qPCR to determine the levels of three different conserved miRNAs: ath-miR160, ath-miR164 and ath-miR159 in Arabidopsis Ler and hen1-1 plants²⁷. As reported earlier, a clear reduction in the levels of mature miR160, miR164 and miR159 was observed in hen1-1 plants in comparison to Ler (control) (Figure 3a). However, it should be noted that the levels of reduction were not uniform for all three miRNAs; fold change of -3.2, -21, and -18 respectively.

Figure 3. Reduced cleavage of target transcripts in miRNA-deficient Arabidopsis hen1-1 mutants.

(a) Expression of ath-miR160, ath-miR164 and ath-miR159 in wild-type Ler and hen1-1 tissues assayed by stem-loop qPCR. The expression levels of the miRNAs were normalized to U6. The data shown are average and error bars indicate the range of possible values based on deviation between Ct values of three replicate assays. (b) Levels of full-length target transcripts (AtARF17 target of miR160, AtNAC1 & AtCUC2 target of miR164 and AtMYB65 target of miR159) in Ler and hen1-1 tissues analyzed using qPCR. Gene expression levels were normalized to U6 and further confirmed using two additional housekeeping genes, AtEF-α and AtGAPDH (data not shown). The data shown are average and error bars indicate the range of possible values based on deviation between Ct values of two replicate assays. (c) The Ratio of Cleaved Transcripts (RCT) over full-length transcripts of AtARF17, AtNAC1, AtCUC2, AtMYB65 in Ler and hen1-1 tissues assayed by qACE. Data shown are average and error bars indicate the range of possible values based on deviation between Ct values of two replicate assays.

Next, we used qPCR and qACE respectively to determine the levels of full length and cleaved transcripts of selected targets of the above miRNAs. AtARF17 is an auxin response factor post-transcriptionally regulated by ath-miR160³². As reported previously, we also observed an increase in the expression of AtARF17 in hen1-1 plants (~1.4-fold increase; Figure 3b). Similarly, AtNAC1 and AtCUC2, targets of ath-miR164 showed ~3.5 fold increase in gene expression levels in hen1-1 compared to Ler (Figure 3b). AtMYB65, known to be regulated by ath-miR159 had a 4.6-fold increase in transcript levels in hen1-1 compared to Ler (Figure 3b). qACE assays indicated a significant reduction in the abundance of miRNA-directed cleavage products for all these target genes in hen1-1. For example, the levels of cleaved AtARF17 levels were lower in hen1-1 plants resulting in a ~3.5-fold reduction in RCT (Figure 3c). Similarly, AtNAC1 had a 6.3-fold reduction in RCT in hen1-1 compared to Ler. AtCUC2 another target of miR164 had only a 2.5-fold decrease in RCT in hen1-1 vs. Ler. Note that both AtNAC1 and AtCUC2 had a 3.5-fold increase in the levels of full-length transcripts (hen1-1 vs. Ler) but the reduction in RCTs was very different. This data indicated that the extent of gene expression regulated by miR164 is perhaps different between these two targets. Possible reasons include differences in cleavage efficiency as well as differences in the extent of overlap in tissue domains where the miRNA and the targets are expressed. qACE helped identify this difference where as simple comparison of full-length transcript levels by qPCR would not have identified it. Finally, AtMYB65 showed a ~3.5-fold reduction in RCT in hen1-1 vs. Ler (Figure 3c).

As a comparison, we determined the ratio of cleaved vs. full length target transcripts using previously published northern blot data using Arabidopsis plants with increased or decreased levels of miR164³³. We performed image analysis to obtain band intensities of full length and cleaved products and calculated the ratio of cleaved transcripts to full length. In Arabidopsis plants over-expressing miR164 in a chemically inducible manner, a 48 hour induction led to a 1.5 fold increase in RCT in one of the transgenic lines (Table S2). In the other transgenic line, the cleaved product was below the limit of detection at this time point and hence RCT could not be calculated. Nevertheless, increased RCT in miR164 over-expressing plants was consistent with results from our soybean miR164 over-expression experiments (Figure 2). Similarly, in mir164 mutant lines, we observed a significant decrease in RCT (Table S3) consistent with results from our hen1-1 experiments (Figure 3). Evidence from these experiments clearly demonstrated the ability of qACE to provide a quantitative indicator of miRNA-directed cleavage of targets resulting from changes in cognate miRNA levels.

Nodule specific expression of NAC1 transcription factors in soybean

Having demonstrated the ability of qACE to detect and quantify cleavage of specific miRNA targets, we used the technique to examine the role of miRNAs in governing nodule-specific/enriched gene expression in soybean. We compiled a total of 326 miRNA genes classified in to 134 families from soybean and predicted a total of 596 genes targeted by these miRNAs in soybean³⁵. Among these, a set of miRNA-target pairs had inverse expression pattern between mature symbiotic nodules (MN) and adjacent root tissues (ABMN; data not shown). We identified two NAC transcription factors, GmNAC1a and GmNAC1b (5’-RACE validated targets of miR164) that were expressed in a nodule-enriched manner. While GmNAC1a had ~6.3-fold higher expression in MN vs. ABMN tissues, GmNAC1b had a 3.2-fold higher expression (Figure 4a). Interestingly, miR164 had a nodule-excluded expression pattern i.e. ~40-fold lower expression in MN vs. ABMN (Figure 4b).

Figure 4. Nodule-enriched gene expression directed by miRNA activity in soybean.

(a) Levels of full-length GmNAC1a & b transcript levels in mature nodules (MN) and adjacent root tissues (ABMN) assayed by qPCR and expression levels were normalized to GmActin. Data shown are average and error bars indicate the range of possible values based on deviation between Ct values of two replicate assays. (b) Level of gma-miR164 in MN and ABMN tissues assayed by stem-loop qPCR. miRNA expression levels were normalized to that of miR1515. Data shown are average and error bars indicate the range of possible values based on deviation between Ct values of three replicate assays. (c) The Ratio of Cleaved Transcripts (RCT) over full-length transcripts of GmNAC1a & b assayed by qACE. Data shown are average from two biological replicates and error bars indicate ± SD.

Based on such inverse expression of these miRNA-target pairs we hypothesized that nodule-specific expression of these target genes might be at least in part regulated by nodule-excluded expression of miR164. In other words, miR164 might actively cleave these target genes in adjacent root tissues to restrict their expression to the nodules. Such a mechanism combined with transcriptional regulation can be used to generate expression gradients as observed between miR166 and its HD-ZIP III proteins during xylem development³⁶. However, mere inverse expression of miRNAs and their targets is not conclusive evidence for our hypothesis. If indeed, these target genes were preferentially cleaved in ABMN tissues to ensure MN-specific expression, we would expect to observe a larger proportion of cleaved transcripts in ABMN tissues vs. MN tissues.

We used qACE to determine the cleavage levels of GmNAC1a and GmNAC1b in ABMN and MN tissues and calculated RCT values. It is worth noting that these targets share the same cleavage signature in degradome/PARE analyses (Table S1). We designed gene specific reverse primers that distinguished these genes (Table S4 and Figure S3a & b). Both GmNAC1a and GmNAC1b had increased cleavage in ABMN tissues compared to MN tissues. Interestingly, consistent with the higher enrichment of GmNAC1a in MN tissues, this gene had a higher cleavage in ABMN tissues (~23.5-fold higher than MN tissues). GmNAC1b also had higher cleavage in ABMN tissues (~15.8-fold higher than MN tissues; Figure 4c). Increased cleavage of GmNAC1a in ABMN compared to GmNAC1b was consistent with increased enrichment of full length GmNAC1a transcripts in MN compared to that of GmNAC1b. qACE assays strongly support the hypothesis that increased cleavage of targets in ABMN tissues by miR164 contributes at least in part to determine their nodule-enriched expression. The high sequence identity between the transcripts of GmNAC1a and b (~92%) precluded the use of Northern to perform such comparisons between these genes.

In summary, we have demonstrated that qACE can detect and quantify miRNA-directed cleavage of specific target genes using miRNA over-expression and miRNA-deficient tissues. We successfully used the technique to identify differential cleavage of specific target genes between nodule and adjacent root tissues revealing an additional miRNA-directed mechanism that results in nodule-specific/enriched gene expression.