In vivo mRNA structure regulates miRNA cleavage in Arabidopsis

: MicroRNA(miRNAs)-mediated cleavage is involved in numerous essential cellular pathways. miRNAs recognize target RNAs via sequence complementarity. In addition to complementarity, in vitro and in silico studies have suggested that RNA structure may influence the accessibility of mRNAs to miRNA-Induced Silencing Complexes (miRISCs), thereby affecting RNA silencing. However, the regulatory mechanism of mRNA structure in miRNA 20 cleavage remains elusive. Here, we investigated the role of in vivo RNA secondary structure in miRNA cleavage by developing the new CAP-STRUCTURE-seq method to capture the intact mRNA structurome in Arabidopsis thaliana . This approach revealed that miRNA target sites were not structurally accessible for miRISC binding prior to cleavage in vivo . Instead, the unfolding of the target site structure is the primary determinant for miRISC activity in vivo . Notably, we found 25 that the single-strandedness of the two nucleotides


Main Text:
MicroRNAs (MiRNAs) are ~21 nucleotide RNAs that impact various aspects of development and stress responses by post-transcriptionally regulating gene expression (1). MiRNAs are loaded onto ARGONAUTE proteins (AGO) to form functional post-transcriptional gene silencing effector complexes, miRNA-Induced Silencing Complexes (miRISCs) (2). miRISC 5 is guided by the miRNA to target RNAs through sequence complementarity and cleave the target RNAs (3,4). However, previous studies found that sequence complementarity is not the sole factor dictating miRNA cleavage (2). The structure of an RNA has been suggested to influence the silencing efficiency (5-7). However, these studies were unable to reveal native RNA structure features adopted through evolution because they indirectly assessed RNA structure by inserting a 10 long sequence that was predicted to form a strong structure, such as a hairpin (5-7). Additionally, these structure assessments examined the target site together with its long flanking regions (5)(6)(7). This led to difficulties in dissecting the individual contributions from the different regions and confounded the identification of a specific RNA structure motif that regulated miRNA cleavage. Furthermore, these in vitro and in silico studies could not reflect the RNA structure folding status 15 in living cells (8-10). Recently, several transcriptome-wide structure probing methods for RNA in vivo have been established (8-10), which provide powerful tools to understand RNA structures under physiological conditions. However, these methods have not assessed the causal relationship between RNA structure and miRNA cleavage because the methods cannot discern whether the obtained RNA structure information is from endogenous degraded RNAs or from intact RNAs. 20 Deciphering the in vivo relationship between mRNA structure and miRNA cleavage requires in vivo structures of target mRNAs before cleavage and the outcome after miRNAmediated cleavage. To obtain the RNA structure of intact mRNA, we performed in vivo selective 2′-hydroxyl acylation analyzed by primer extension (SHAPE) chemical probing on Arabidopsis thaliana (A. thaliana) with optimized conditions (Fig. 1A and Fig. S1A). Next, we enriched the 25 intact mRNAs through the terminator exonuclease treatment (Fig. S1B) (11), and then polyA+ purification to remove the degraded mRNAs. We generated two independent biological replicates of (+)SHAPE (samples with SHAPE treatment) and (-)SHAPE (control samples without SHAPE treatment) libraries according to the protocol described previously (8,12). Between 90-97% of 340~380 million reads were mapped onto mRNAs ( Fig. S1C and D) with the reproducibility of 30 the CAP-STRUCTURE-seq library confirmed by comparing the two biological replicates ( Fig.  S2A and B). Nucleotide occurrence was consistent in both (-)SHAPE and (+)SHAPE libraries (Fig. S1E). To validate CAP-STRUCTURE-seq, we compared the SHAPE reactivity of the 18S rRNA with the corresponding phylogenetic covariance structure (Fig. S3A) and the 3D structure (Fig. S3B). We found that CAP-STRUCTURE-seq can accurately probe RNA structure in vivo, 35 and in addition it outperforms the previous dimethyl sulphate (DMS)-based method, STRUCTURE-seq (8) (Fig. S1F).
To further validate CAP-STRUCTURE-seq we performed meta-property analyses with over 16,576 transcripts of sufficient mRNA structure information (Fig. S4A). We found that the region immediately upstream of the start codon showed particularly high SHAPE reactivity (Fig. 40 S4B). This result further supports the notion that less structured regions near the start codon may facilitate translation (13,14). We also found a periodic reactivity trend in the coding sequence (CDS) but it was absent in the untranslated regions (UTRs) (Fig. S4C), which is consistent with previous studies (8,15). Similar to a RNase-based structure study in human (15), we also observed an asymmetric mRNA structure signature at the exon-exon junction in A. thaliana (Fig. S4D). These conserved mRNA structure features suggest that CAP-STRUCTURE-seq successfully provides a global mRNA structure model in plants. Moreover, existing RNA structure methods (9, 10, [16][17][18][19] are not able to discern whether the RNA structure information belongs to the 5 endogenous degraded mRNAs or to the intact mRNAs (Fig. S4E). Additionally, degraded mRNAs are capable of introducing false positive signals in the previous methods ( Fig. S4E) (9, 10, 20, 21). For instance, the miRNA cleavage products led to a skewed DMS reactivity profile at the miRNA cleavage sites in the previous DMS-based method (8) (Fig. S4F), while no such false positive signal occurs in our CAP-STRUCTURE-seq (Fig. S4G). Thus, our CAP-STRUCTURE-seq 10 provides RNA structure information for intact mRNAs while excluding degradation signals ( Fig.  S4E and G), thereby overcoming the limitations of previous methods that include degradation products ( Fig. S4E and F) (8, 10, 18,22,23).
Having developed a reliable method to accurately probe the intact mRNA structure, we then focused on estimating the in vivo miRNA cleavage efficiency (CE) at a genome-wide scale. 15 Previously, qualitative measurement of miRNA cleavage events was achieved from degradome library analysis (24)(25)(26). Inspired by the definition of enzymatic activity (27), we quantitatively estimated CE by measuring how many degradation products were generated from one unit of substrate mRNA by one unit of miRNA (Methods). We generated degradome libraries to estimate the degradation products ( Fig. S5A and Methods) and miRNA-seq libraries to estimate miRNA 20 abundance ( Fig. S5A and Methods), with library reproducibility confirmed by comparing the biological replicates (Fig. S2C, S2D and Methods). We then combined the degradome, (-)SHAPE and miRNA-seq libraries to estimate CE ( Fig. S5A and Methods). We verified the consistency of our CE against previously reported targets (Table S1). For example, the CE of AP2 targeted by miRNA172, which has been shown to act through translational repression rather than mRNA 25 cleavage (28)(29)(30), was zero as expected (Table S1). SNZ (Table S1) is another target of miRNA172, and showed no evidence of miRNA cleavage, consistent with the previous result (31). In contrast, TOE2, which is cleaved by miRNA172, had relatively high CE (Table S1). Additionally, TAS1a and TAS2, which must be cleaved by miRNA173 to then serve as templates for trans-acting siRNA (tasi-RNA) (32), had high CE (Table S1). These observations were 30 consistent with their previous reported biological functions (28)(29)(30)(31)(32). Then, we systematically examined the relationship between sequence complementarity and CE. Globally, we found that sequence complementarity and CE were uncorrelated (Spearman correlation -0.015, Fig. S5E). In addition, targets with mismatches and/or GU wobble pairs (sequence complementarity penalty score, SCPS>0) were sometimes more effectively cleaved than targets with perfect 35 complementarity (SCPS=0) (Fig. S5E). Our results indicate that besides sequence complementarity between miRNA and mRNA there may be other factors affecting CE, for example, mRNA structure. In summary, both the RNA structure of the intact mRNAs and miRNA cleavage in vivo can be quantitatively measured.
With CAP-STRUCTURE-seq elucidating the RNA structure, we could begin to answer the 40 elusive question about whether miRNA target sites were structurally accessible in vivo. Since our CAP-STRUCTURE-seq directly measures the in vivo structural accessibility via SHAPE reactivity (33), we assessed the SHAPE reactivity profiles across the miRNA target sites on the intact mRNAs. SHAPE reactivities of the target sites showed no significant difference from the upstream region (two-sided Mann-Whitney-U test, P value is 0.31, Fig. S4G), and were lower than the downstream region (one-sided Mann-Whitney-U test, P value is 0.0038, Fig. S4G). These features indicate that under physiological conditions the target sites are not highly accessible, which may provide a protective mechanism for target sites, mitigating against other cellular ribonucleases. These relatively inaccessible target sites prompted us to ask whether the target site structure affects miRNA cleavage in vivo. To address this question, we examined two alternative energetic 5 landscapes associated with the miRISC cleavage process in vivo: an enzyme-limiting scenario and a structure-limiting scenario (Fig. 1B). In the enzyme-limiting scenario, the energy barrier (ΔG ⧧ open) between the inaccessible and accessible structural states (i.e., the unfolding of the target site) is lower than the barrier for catalytic cleavage (black line in Fig. 1B). Thus, the target sites equilibrate quickly between inaccessible and accessible structural states during the binding step 10 prior to the catalytic step of miRNA cleavage. In this scenario, the CE would vary with the free energy required to surmount the AGO catalytic barrier, ΔG ⧧ cutting (Methods), and would be less affected by the RNA structure of the target site. In the structure-limiting scenario, the energy barrier (ΔG ⧧ open) between the inaccessible and accessible structural states is higher than the barrier for cleavage (red line in Fig. 1B). Therefore, the target sites cannot achieve equilibrium binding 15 with miRISC before catalytic cleavage. In this scenario, CE would vary with the free energy of opening the target site structure, ΔG ⧧ open, rather than ΔG ⧧ cutting. We used our in vivo structures to computationally approximate these two scenarios and explored a range of flanking lengths upstream and downstream of the target site (Fig. 1C, D and Fig. S6A, B). Analysis of our SHAPE reactivity-informed structures revealed that, for most flank sizes, CE anti-correlated with 20 ΔG ⧧ open with a broad maximum centered around flanks of 50 nucleotides upstream and downstream (Spearman correlation of -0.23, P = 6.3e-9) ( Fig. 1C and Fig. S6A). However, for most flank sizes, CE had no correlation with ΔG ⧧ cutting ( Fig. 1D and Fig. S6B), contrary to the reaction kinetics where the energy barrier is anti-correlated with reaction processivity. These results indicate that the unfolding of target site structure is the rate-limiting step that determines miRISC activity in vivo. 25 Furthermore, this structure-limiting scenario reveals that the ribonuclease AGO undergoes "sticky regime" activation (34), where substrate mRNAs associate and dissociate with AGO more slowly than they are being cleaved. This scenario further explains previous AGO RIP-seq results, where target transcripts were only captured in the slicing-defective AGO1 mutant (AGO1 DAH ) but not in wild type AGO (35). We found that ΔG ⧧ open anti-correlated with the enrichment ratio of target 30 RNAs associated with AGO1 protein from previous AGO1 DAH -RIP-seq results (35) (Spearman correlation -0.21, P = 0.05). In contrast, the free energy of binding of the miRNA-target duplex (ΔGduplex) and ΔG ⧧ cutting show no correlations (Spearman correlation 0.06 with P = 0.32 and -0.11 with P = 0.16, respectively). These observations suggest that the target sites are not structurally accessible in vivo, but rather the unfolding of the target site structure is the primary determinant 35 for target RNA processing by AGO.
Having revealed that the target site structure affects cleavage in vivo, we then investigated whether the structure of the target site flanking regions is involved with miRNA cleavage. We assessed the RNA secondary structure by separating the miRNA targets into non-cleaved (zero CE) and cleaved (positive CE) groups. We found higher SHAPE reactivity at the +1 and +2 40 nucleotides immediately downstream of target sites in the cleaved group relative to the noncleaved group ( Fig. 2A), suggesting that these two nucleotides are more single-stranded than their neighbors. To confirm this observation, we used the SHAPE reactivity with the ViennaRNA RNAfold utility (36) to calculate the base pairing probabilities (BPPs) for each nucleotide. We found that the BPPs of the +1 and +2 nucleotides were much lower than their neighboring 45 author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.21.885699 doi: bioRxiv preprint nucleotides (Fig. S7), indicating an increased likelihood of single-strandedness in the cleaved group compared to the non-cleaved group. Furthermore, the single-strandedness of the two nucleotides was unlikely to be due to sequence composition (Fig. 2B) or AT content (Fig. 2C) because there was no difference between the non-cleaved and cleaved groups. Our results reveal that a secondary structure feature, specifically single-strandedness of the two nucleotides adjacent 5 to the 3' end of the miRNA target site, generally exists in vivo in intact mRNAs that will undergo cleavage. We named this structure feature 'Target Adjacent structure Motif' (TAM).
To explore the functional role of TAM in miRNA cleavage, we designed a structure assay (Methods) by concatenating miRNA156 target sites (20 nt) with a designed stable structure module (either a G-quadruplex structure or a stem-loop) to mimic the base-pairing state of the two 10 nucleotides immediately downstream of the target site ( Fig. 3A and Fig. S8A). To maintain the single-strandedness of the TAM we inserted two Adenines (AA), or "slippery sequence", immediately between the target site and the designed stable structure module ( Fig. 3A and Fig.  S8A). We confirmed the formation of TAM in vivo by using capillary electrophoresis (37) to resolve the in vivo RNA structure for each design (Methods, Fig. 3B and Fig. S8B). We then 15 assessed the miRNA cleavage in vivo by measuring the levels of non-cleaved substrate mRNA. We found that the non-cleaved mRNA level in the substrates with TAM was significantly lower than those without TAM ( Fig. 3C and Fig. S8C), which suggested increased cleavage when TAM was present. To further confirm whether the presence of TAM exclusively determines target cleavage, we performed an in vitro AGO cleavage assay using HA immuno-affinity-purified 20 wildtype AGO protein and we found that the target RNA was cleaved only when TAM was present (Fig. 3D). Our results reveal that TAM is essential for miRISC nuclease activity. One might expect TAM in the target mRNA to facilitate AGO binding instead of directly triggering the nuclease activity of AGO proteins. To test the possibility that TAM affects target binding, we conducted an in vivo binding assay (Methods) by using the slicing-defective AGO1 mutant (AGO1 DAH ). We 25 found that AGO1 was able to bind the target RNAs with the same binding affinity irrespective of whether the TAM was present or absent ( Fig. 3E and F). Therefore, our data reveal that TAM regulates miRISC cleavage activity rather than affecting target binding.
Following our investigation into the functional role of TAM, we further examined the biological impact of TAM by using a native target gene. An important flower development gene, 30 APETALA2 (AP2), is targeted by miRNA172 and the interaction between AP2 and miRNA172 leads to translational inhibition rather than cleavage (29). In our CE results, AP2 cleavage efficiency was zero, representing non-cleavage. Interestingly, the in vivo RNA structure of AP2 in our CAP-STRUCTURE-seq shows that the two nucleotides downstream of the target site (+1 and +2 positions) were base-paired with the 879 th and 880 th nucleotides in AP2, indicating that TAM 35 was absent (Fig. 4A). To assess whether the non-cleavage of AP2 is due to the absence of TAM, we introduced TAM by mutating the A at the 879 th position with G, U, and C (labelled A879G, A879U, A879C) to disrupt the base-pairing states while preserving the codon. We then transformed the wild-type construct without TAM (non-TAM) and the synonymous mutation constructs with TAM (TAM) into protoplasts of the stable MIR172 over-expression line. We 40 measured the miRNA cleavage of AP2 by primers designed to span the cleavage site to detect the non-cleaved mRNAs (Fig. 4B). We found increased cleavage of AP2 (measured by decreased levels of non-cleaved AP2) when the TAM was present compared to when TAM was absent (Fig.  4B), indicating the TAM can switch on the miRNA172 cleavage of AP2. Therefore, our results suggest that the TAM is able to change the cleavage fate of target genes. 45 author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.21.885699 doi: bioRxiv preprint To further assess the genetic effects of TAM, we performed stable transformations of these non-TAM and TAM constructs (A879G, A879U, A879C) into the loss-of-function ap2-5 mutant. As a floral homeotic gene, AP2 determines floral organ identity and the ap2-5 mutant results in the replacement of the perianth organs by the reproductive organs ( Fig. 4C and D). In the non-TAM transgenic lines, floral organ defects were restored by complementing ap2-5 mutants with 5 the wild-type AP2, in which no TAM formed ( Fig. 4E and F). In contrast to the non-TAM transgenic lines, the severe floral organ defects of ap2-5 could not be restored when ap2-5 mutants were transformed with the TAM constructs (A879G, A879U, A879C) ( Fig. 4G-4L). We further assessed the miRNA cleavage efficacy of AP2 where the phenotypic restorations are different between non-TAM and TAM transgenic lines. We measured the levels of non-cleaved AP2 in 10 these transgenic lines. We found increased cleavage of AP2 (measured by decreased levels of noncleaved AP2) in the TAM transgenic plants compared to the non-TAM transformants (Fig. 4M), while the miRNA172 levels were largely unchanged between non-TAM and TAM transgenic plants (Fig. 4N). Taken together, our results indicate that switching on TAM is sufficient to turn on the miRNA cleavage of AP2, and so significantly affects plant development. 15 In summary, by deciphering intact mRNA structures in vivo through CAP-STRUCTUREseq, we found that miRNA target sites are not structurally accessible in vivo and we demonstrated that the unfolding of the miRNA target site structure predominantly affects miRISC activity in vivo. Furthermore, we discovered that the native RNA structure motif, TAM, is sufficient to regulate miRNA cleavage in vivo. The mechanism that we found here provides evidence of mRNA 20 structure-dependent regulation of biological processes in vivo. Our study reveals that in vivo mRNA structure serves as an additional regulator of miRISC activity, which will facilitate the biotechnological engineering of gene silencing, and possibly provide an additional avenue towards crop improvement. 25 author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.21.885699 doi: bioRxiv preprint Fig. 1. Measurement of in vivo RNA structure of intact mRNA and in vivo miRNA cleavage efficiency. A, Schematic of (+)SHAPE sequencing library generation showing NAI treatment, nucleotide modification and purification of intact mRNA steps. A. thaliana etiolated seedlings were treated with NAI. After extraction of total RNA, degraded mRNAs (dark yellow) were removed leaving intact mRNAs characterized by 5'CAP and 3' poly(A) tails (dark blue). cDNAs 5 were obtained and subjected to an established library construction. B, The miRISC cleavage reactions include target binding and cleavage catalysis. Two alternative scenarios demonstrate the energetic landscape of miRNA cleavage (black and red). In the enzyme-limiting scenario (black), target site structure equilibrates quickly between inaccessible (closed) and accessible (open) states in the binding step compared to the catalytic step of miRNA cleavage. In this scenario, the apparent 10 activation energy is ΔG ⧧ cutting, which measures the energy required to raise the initial substrate target RNA to the transition catalysis-compatible state. Alternatively, in the structure-limiting scenario (red), the target site cannot achieve equilibrium binding before cleavage. In this scenario, the energy barrier between the target site and the transient state is higher than the barrier for cleavage. And the apparent activation energy is equal to ΔG ⧧ open, which measures the energy 15 required to open the target site structure. C, Spearman correlation between ΔG ⧧ open and cleavage efficiency (647 target sites with the upstream and downstream flank lengths of 50 nucleotides, P = 6.3e-9 ***). D, A similar analysis to C, but for ΔG ⧧ cutting and cleavage efficiency (P = 0.46). author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.21.885699 doi: bioRxiv preprint . The profiles show the per nucleotide mean +/-SEM across transcripts, aligned by target site start (left panels) and end position (right panels). Two nucleotides in the cleaved group (positive CE group), immediately downstream of target sites (TAM region), show significantly higher SHAPE reactivities compared to their neighbors (by Mann-Whitney-U tests). Compared to the upstream region (target sites), P = 4.8e-7*** for both 1 st and 2 nd nucleotides; Compared to the downstream region, P = 3.9e-3** for both 1 st and 2 nd nucleotides. The two individual nucleotides of the TAM region in the non-cleaved group (zero CE group) are not significantly higher than their neighbors by Mann-Whitney-U tests. B, Sequence composition around the target sites for the total cleaved (675 target sites) and non-cleaved groups (571 target sites). C, AT content around the target sites for the total cleaved (675 target sites) and 15 non-cleaved (571 target sites) groups. author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.21.885699 doi: bioRxiv preprint Fig. 3. Validation of TAM functionality by a designed structure assay. A, Cartoon representation of the protoplast transformation assay to validate the TAM functionality using a designed structure assay. GQS refers to a G-quadruplex. The miRNA156 target sites (blue comb) followed by 0 or 2 Adenines (As) and ending with a GQS. The prefixes, "0A" and "2A", indicate 5 the number of Adenines. B, In vivo RNA structures of 0A_GQS and 2A_GQS. C, The non-cleaved mRNA abundance for the substrates in B was measured by qRT-PCR (dark yellow bars). The antisense target sites were used as controls (teal bars). Data are mean +/-SEM from three independent biological replicates. D, In vitro AGO1 cleavage assay shows that target RNA is cleaved when the TAM is present. The target RNAs were incubated for 1h or 2h, where two 10 cleavage products were present in the target with TAM on the X-ray film (as indicated). The asterisk indicates the background bands present in both control and experiment groups. E, In vivo AGO1 binding assay shows no difference between target RNAs with TAM and without TAM. RENILLA and ACTIN were used as the control. miRNA156 levels were measured in all the samples. F, The RNA abundance enrichment in E was quantified by amplicon intensities and normalized 15 by input. Data are mean +/-SEM from three independent biological replicates.
author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.21.885699 doi: bioRxiv preprint Acknowledgments: We thank Dame Prof. Caroline Dean (John Innes Centre), Prof. Giles Oldroyd (SLCU, Cambridge), Dr. Desmond Bradley (John Innes Centre) and the group members in Ding lab for discussions with this work. We thank Mr. Mirko Ledda (UC. Davis) for discussions on data analysis. We thank Jim Carrington's lab for providing us with the AGO construct for performing the binding assay.  The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.21.885699 doi: bioRxiv preprint

Plant materials and growth conditions
A. thaliana seeds of the Columbia (Col-0) and the xrn4 mutant accession (38, 39) were sterilized with 70% (v/v) ethanol and plated on half-strength Murashige and Skoog medium (1/2 MS). The plates were wrapped in foil and stratified at 4℃ for 3-4 days and then grown in a 22-5 24℃ growth chamber for 5 days.

Degradome library construction
Five-day-old A. thaliana etiolated seedlings were grown as described above. They were collected and immediately frozen in liquid nitrogen and stored at −80°C. The seedlings were ground into powder. Total RNA was extracted using RNeasy Plant Mini Kit (Qiagen) according 10 to the manufacturer's instructions. On-column DNAaseI treatment was carried out according to RNase-Free DNase Set (Qiagen). To construct the Illumina library for degradome analysis, polyA+ selection was carried out using the Poly(A)Purist Kit (Ambion TM ). Selectively captured polyadenylated RNAs (1μg) were ligated directly to an DNA/RNA hybrid adapter (5'-CTACAC GACGCTrCrUrUrCrCrGrArUrCrUrNrNrN-3') using T4 RNA ligase (NEB) at 37°C for 30 15 minutes. The ligated RNAs were subjected to RT by SuperScript III First-Strand Synthesis System (Invitrogen) with random hexamers fused with Illumina TruSeq adapters (5'-CAGACGTGTGCTCTTCCGATCTNNNNNN-3'). PCR amplification was performed on the ligated cDNA using Illumina TruSeq Primers. Two different barcode indices were used for two degradome biological replicates. The final dsDNA degradome libraries were subjected to next-20 generation sequencing on Illumina HiSeq 4000.

miRNA library construction
The same seedling samples stored at −80°C, as described above, were ground into powder using liquid nitrogen. Total RNA was extracted using mirVana miRNA Isolation Kit (Ambion TM , Austin, TX, USA) following the manufacturer's instructions. The integrity analysis was performed 25 on a Bioanalyzer by the Beijing Genomics Institute (BGI), Shenzhen, China, which also performed the library construction according to standard protocols.

Gel-based 18S rRNA structure probing
The gel-based method of structure probing used the same in vivo total RNA pools as for CAP-STRUCTURE-seq. To accomplish gel-based structure probing, reverse transcription was 30 performed using 18S gene-specific DNA primers with 5' end labelled Cy5 (TAGAATTACTACGGTTATCCGAGTA). The whole procedure was performed according to Ding,et.al(8). Each gel was detected by Typhoon FLA 9500 (GE Healthcare).

(+)SHAPE and (-)SHAPE CAP-STRUCTURE-seq library construction
We modified the in vivo chemical probing protocol (8) by changing the reagent from 35 dimethyl sulphate (DMS) to the SHAPE reagent, 2-methylnicotinic acid (NAI). NAI was prepared author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.21.885699 doi: bioRxiv preprint as reported previously (40). Briefly, five-day-old A. thaliana etiolated seedlings were suspended and completely covered in 20 ml 1X SHAPE reaction buffer (100mM KCl, 40mM HEPES (pH7.5) and 0.5mM MgCl2) in a 50 ml Falcon tube. NAI was added to a final concentration of 144mM and the tube swirled on a shaker (1,000rpm) for 15min at room temperature (22°C). This NAI concentration and reaction time had been optimized to allow NAI to penetrate plant cells and 5 modify the RNA in vivo under single-hit kinetics conditions (Fig. S1A). After quenching the reaction with freshly prepared dithiothreitol (DTT), the seedlings were washed with deionized water and immediately frozen with liquid nitrogen and ground into powder. Total RNA was extracted using RNeasy Plant Mini Kit (Qiagen) according to the manufacturer's instructions, followed by on-column DNaseI treatment in accordance with the manufacturer's protocol. The 10 control group was prepared using DMSO (labelled as (-)SHAPE), following the same procedure as described above.
To capture the structure information around the cleavage site of miRNA target genes, we adopted the feature of 5PSeq (11). The whole CAP-STRUCTURE-seq procedure is illustrated in Fig. 1A. In our method, the (+)SHAPE and (-)SHAPE RNA samples were treated with 15 Terminator™ 5´-Phosphate-Dependent Exonuclease (TER51020, EPICENTRE co.), which processively digests RNA with 5´-monophosphate ends, thereby leaving mRNAs with 5'cap structures. Following the 5'cap enrichment, polyA+ selection was carried out using the PolyA purist Kit (AmbionTM) leaving intact (pre-cleaved) mRNAs with enriched 5'cap and 3'poly(A) tails. The resultant mRNAs were subjected to library construction following the STRUCTURE-20 seq procedure. The name of CAP-STRUCTURE-seq refers to 5' cap-enriched and 3' poly(A)enriched RNA structure sequencing.

Degradome analysis
Raw reads were processed to remove 5'and 3' adapter sequences. Degradome reads were mapped to the TAIR10 transcript reference and a degradome density file was generated. The 25 degradation level of target genes was normalized by reads per kilobase per million mapped reads (RPKM).

miRNA-seq analysis
The small RNA sequences were processed by BGI to filter out the 5′ adapter sequences, 3′ adapter sequences and low-quality reads. We mapped two biological replicates against 253 30 miRNA sequences confidently annotated as A. thaliana mature miRNAs (41). We used Bowtie (42) for the mapping using the command 'bowtie -f -a -S --best --strata -v 1'. pysam (42) was used to count the mapped reads.

CAP-STRUCTURE-seq analysis
We merged the biological replicates of the transcript-level reverse transcription (RT) stop 35 counts to obtain a single (-)SHAPE library and a single (+)SHAPE library. We calculated the SHAPE reactivity using a slightly modified version of the formula in Ding et al. where Pi is the (+)SHAPE RT count and Mi is the (-)SHAPE RT count at nucleotide i. The factor, α (= min (1, ∑ log(1 + 1 ) 1 / ∑ log (1 + 1 ) 1 ) is a simple library size correction factor. Setting α=1 recovers the reactivity formula in Ding et al.(8). The reactivities were then normalised using the box-plot method (43). For the SHAPE reactivity profiles, we extracted values in the 50 nts 5 upstream and downstream of target sites and calculated a per nucleotide mean and SEM.
The cleavage efficiency (CE) calculation miRNAs regulate mRNAs through translational repression, mRNA de-stabilization and 10 mRNA cleavage. In animals, it has been suggested that translational repression is prevalent, which can then be followed by de-stabilization, such as shortening of the poly(A) tails at the 3' end and removal of the cap at the 5' end (44). In plants, cleavage is the dominant pathway for miRNAs to regulate their target mRNAs (45). In order to compare the ability of different miRNAs to cleave their targets on a global scale, we need to quantify the cleavage efficiency (CE) of miRNAs at their 15 target sites.
Our CE calculation is based on two underlying facts (41,46): miRNA-mediated cleavage is the major mRNA turnover pathway for target genes; the 5' cleaved products are located within binding sites, which are transiently stable. Therefore, the degradation signal within target sites reflects the cleavage products from miRISC cleavage. These two facts were also confirmed by our 20 analysis below.
Firstly, to confirm that the degradation signal within target sites is mainly from miRNAmediated cleavage, we mapped the 5'end of our wildtype (WT) A. thaliana degradome reads to previously validated cleavage sites (41). We found that most of the read ends were mapped at the tenth nucleotide of the miRNA complementary sites (Fig. S5C), which provides strong evidence 25 that miRNA-mediated cleavage is the major mRNA turnover pathway for target genes. Secondly, to confirm that the 5' cleaved products are transiently stable, we mapped the read ends of the cleavage products in the xrn4 mutant to previously validated cleavage sites (41). The cleavage site distribution in the xrn4 mutant exhibited the same pattern as WT (Fig. S5D), which is consistent with the notion that miRNA cleavage products are temporally stable intermediates, resistant to 30 cellular XRN4 exonuclease in A. thaliana, although the precise mechanism is currently unknown (46,47). Thus, we counted the degradation reads within target sites as the outcomes of cleavage products from miRISC cleavage.
Based on these two observations, we can derive the miRNA CE in vivo. Since AGO and miRNA form an enzyme complex, we defined the (CE) in a similar way to enzyme activity. In 35 details, the catalytic ability of an enzyme can be defined as the amount of product generated by one unit of enzyme from one unit of substrate, which led us to define: author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.21.885699 doi: bioRxiv preprint Cleaved transcripts can be characterized by the 5' phosphate featured on 3' cleavage products. We constructed a degradome library to capture cleavage products. RNA abundance of degradation products can be measured by the Reads Per Kilobase of transcript per Million mapped reads (RPKM). RPKM is a measure of relative RNA concentration in the whole transcriptome. In the degradome, the RPKM of each mRNA means relative degraded mRNA fragment number 5 compared to all degradation products, i.e., [ ] ≈ .

10
To quantify the miRNA-mediated degradation products, we used TargetFinder (48) to predict the miRNA target sites on the expressed transcripts in our samples and removed any duplicated target sites from the same miRNA family. TargetFinder predicts target sites with high specificity in A. thaliana by assigning a sequence complementarity penalty score (SCPS) (49) (Fig.  S5B). We then counted and designated the reads mapped within each target site as the products of 15  The total degradation products, total RNAs and total miRNAs should be constant and be reflected by the library sequencing depth. Therefore, α is a constant.
The population of each target mRNA is constant over time due to the dynamic equilibrium of an intact mRNA and its degraded products (50), i.e., 40 author/funder. All rights reserved. No reuse allowed without permission. In short, we can combine miRNA-seq, (-)SHAPE and degradome libraries to quantify miRNA CE in vivo. The CE pipeline is illustrated in Fig. S5A.
Multiple sequence alignment and Visualization of RNA structure Multiple sequence alignment was conducted with Jalview (51) and visualized by WebLogo 25 (52). The secondary structures and the corresponding SHAPE values were visualized using the VARNA Java applet (53). ΔG ⧧ open was computed as the difference between the minimum free energy of the in vivo secondary 30 structure and the minimum free energy of the ''hard constrained'' transcript, in which the target nucleotides were required to be unpaired (6,54). By exploring a range of flanking region lengths upstream and downstream of the target site, we chose the upstream and downstream flank lengths to be 50 nucleotides for the majority of analyses. We used RNAfold from the Vienna RNA package (55) together with our SHAPE reactivity data to calculate the energy terms in ΔG ⧧ open, the RNA 35 structures and the base pairing probabilities (BPP).

Calculation of ΔG
ΔG ⧧ cutting measures the energy required to raise the initial substrate target RNA to the transition catalysis-compatible state. ΔG ⧧ cutting is given by: author/funder. All rights reserved. No reuse allowed without permission. where ΔGduplex is the binding free energy of the miRNA-target duplex and ΔGcatalysis refers to the miRISC transition catalytic state energy. The ΔGduplex energy was calculated for the miRNA sequence and the target region sequence using RNAduplex from the Vienna RNA package (55).

Plasmid Construction
For cleavage efficiency validation, the miRNA156 target sites, followed by 0 or 2 Adenines (As) and ending with a G-quadruplex (GQS) or a stem-loop (SL) were synthesized and inserted into AflII and PacI of Firefly 3'UTR in vector inter2. We labelled the GQS constructs as 0A_GQS 15 and 2A_GQS, and the stem-loop constructs as 0A_SL and 2A_SL, with the prefix indicating the number of Adenines. Antisense of miRNA156 target site constructs with the same flanking sequence were also synthesized as the control for each construct.
To carry out the AP2 mutagenesis assay, the native AP2 sequence was inserted into the pEAQ vector by Invitrogen LR reaction. Mutations were introduced using Q5 site-directed 20 mutagenesis kits (NEB, E0554S).
For the miRNA156 overexpression vector construction in AGO1 in vivo binding assay, the MIR156B genomic sequence was inserted into AscI and SacI of vector pMDC32. Primers are listed in Table S2.
Arabidopsis protoplast transformation 25 Protoplasts from the stable MIR156 over-expression line and the stable MIR172 overexpression line were prepared and transformed according to the Tape-Arabidopsis Sandwich method (61). 16 h after transformation, protoplasts were centrifuged at 100 g for 2 min. RNA was extracted with Qiagen RNeasy kit and qRT-PCR quantification was performed with Bio-Rad CFX. Primers are listed in Table S2. 30 In vivo structure validation experiments Four-week-old tobacco leaves were co-infiltrated with agrobacterium strains harboring plasmids of 0A_SL, 2A_SL, 0A_GQS or 2A_GQS. Two days after infiltration, the leaves were treated with 144 mM of the SHAPE reagent (NAI). The control group was treated with DMSO. Total RNA was extracted using RNeasy Plant Mini Kit (Qiagen) and then by on-column DNaseI 35 author/funder. All rights reserved. No reuse allowed without permission.
The ligated cDNA samples were dissolved in 10 μl of water and used for the PCR reaction. The PCR reaction contained final concentrations of 0.5 mM VIC-labelled DNA gene-specific primers (the same as that used in the reverse transcription primers except the 5' end was labelled with Vic), 0.5 mM of linker reverse primer (AGATCGACTCAGCTACACGACGC), 200 mM 15 dNTPs, 1X ThermoPol reaction buffer and 1.25U of NEB Taq DNA polymerase in 25 μl. The solution was then extracted with Phenol:Chloroform:Isoamyl (25:24:1, v/v, sigma) and incubated with Mly1 restriction enzyme, according to the manufacturer's protocol. Finally, the reaction pellets were dried and resuspended in Hi-Di formamide (Applied Biosystems/Life Technologies).
The Ned-labelled gene-specific primer (the same as that used in reverse transcription 20 primer except the 5' end was labelled with NED) was used to make sequencing ladders using linear DNA and 1 µl 5 mM ddTTP by Klenow DNA Polymerase I (New England Biolab) (62). Then, the reaction pellets were dried, resuspended in Hi-Di formamide (Applied Biosystems/Life Technologies) and run on an Applied Biosystems 3730xl Genetic Analyzer. The resulting data were analyzed using QuSHAPE (63). 25 AGO1 in vitro cleavage assay HA-tagged AGO1 WT was immuno-purified from Arabidopsis seedlings (64). The 0A_GQS and 2A_GQS designed RNAs were transcribed in vitro with T7 polymerase (NEB, 2040S) as substrates. To perform the slice assay, cleavage buffer (100mM ATP, 10mM GTP, 60mM MgCl2, 0.5M CPO4, 1mg/ml CPK) was added to 20μl beads in extraction buffer (1:1) bearing freshly 30 purified HA-AGO1 from 3g seedling on the beads' surface. 50 cps of labelled substrate was added to the reaction and incubated at 25℃. 10μl of the resultant liquid was added to 10μl 2x RNA loading buffer (95% Formamide, 0.02% SDS, 1mM EDTA, 0.02%, Bromothymol Blue, 0.01% Xylene Cyanol), denatured for 5min at 95℃ and loaded into a 1mm PAGE gel (10% acrylamides:bis 19:1, 7M Urea, 1xTBS). Then the gel was dried and exposed to a phosphor screen 35 for image analysis.
author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.21.885699 doi: bioRxiv preprint AGO1 in vivo binding assay Four-week-old tobacco leaves were co-infiltrated with agrobacterium strains harbouring plasmids of 35S:MIR156B and 35S:HA-AGO1 DAH and 0A_GQS or 2A_GQS. Two days after infiltration, the leaves were collected and ground in liquid nitrogen. The protein/RNA complexes were extracted using two volumes of IP buffer (50 mM Tris-HCl pH 7.5, 150 mM NaCl, 5% β-5 mercaptoethanol, 1 mM EDTA, 10% glycerol, 0.1% NP-40, 1 mM PMSF, and 1X complete protease inhibitor cocktail). After removing insoluble debris by centrifugation, cell extracts were incubated with anti-HA antibody (Abcam ab9110) for 1h at 4℃ with gentle mixing. The anti-HAdecorated extracts were then incubated with pre-washed protein G magnetic beads for 1h. After incubation, the beads were washed 6 times with the IP buffer. The RNA produced after co-10 immunoprecipitation was precipitated with ethanol and glycogen, and analysed by RT-PCR. The miRNA156 expression levels were analysed by miRNA RT-PCR (65). author/funder. All rights reserved. No reuse allowed without permission.