Introduction

Ribozymes (catalytic RNAs) are thought to have dominated early life forms (1) until they were gradually replaced by more effective and stable proteins. One of the few that remain active is self-cleaving ribozymes found in a wide range of living organisms (2), involved in rolling circle replication of RNA genomes (3–5), biogenesis of mRNA and circRNA, gene regulation (6–9), and co-transcriptional scission of retrotransposons (10, 11).

So far, there are a few self-cleaving ribozymes found in the human genome. They include two hammerhead ribozymes (7, 12), a HDV-like cytoplasmic polyadenylation element-binding protein 3 (CPEB3) ribozyme within the transcript of a single-copy gene (8), and self-cleaving ribozymes associated with olfactory receptor OR4K15 (olfactory receptor family 4 subfamily K member 15), insulin-like growth factor 1 receptor (IGF1R), and a LINE-1 (Long Interspersed Nuclear Element-1) retrotransposon (8). All these ribozymes were discovered by using selection-based biochemical experiments (8). More recently, the hovlinc ribozyme (13) was identified by genome-wide screening against RNAs with the specific 5’-hydroxyl termini resulted from self-cleavage, followed by a combination of RNA structure prediction as well as deletion analysis and mutation-based validation for determining its secondary structure. It has been shown that single nucleotide polymorphism (SNP) in CPEB3 ribozyme was associated with an enhanced self-cleavage activity along with a poorer episodic memory (14). Inhibition of the highly conserved CPEB3 ribozyme could strengthen hippocampal-dependent long-term memory (15, 16). However, little is known about the other human self-cleaving ribozymes.

Previous studies on self-cleaving ribozymes found in retrotransposons suggested their importance in the overall size, structure, and function of different genomes. Most selfcleaving ribozymes found in retrotransposons belong to two widespread ribozyme families: HDV-like ribozymes in R2 elements (10) and LITc retrotransposon (11) and hammerhead ribozymes (HHR) in short interspersed nuclear elements (SINEs) of Schistosomes (17), Penelope-like elements (PLEs) (18, 19) and retrozymes (20). Unlike the wide occurrence of hammerhead ribozymes in different types of retrotransposons, HDV-like ribozymes usually locate in the 5’ UTR region of the retrotransposons only, and were found in different species including Drosophila melanogaster (21, 22), Bombyx mori (23), arthropods, nematodes, birds and tunicates (24, 25). These ribozymes are important for co-transcriptional processing of the full-length transcript and translation of the downstream open reading frame (ORF). This suggests that the LINE-1 ribozyme found in the human genome may also play an active role in co-transcriptional processing during evolution.

Typically, self-cleaving ribozymes are in noncoding regions of the genome. For example, the CPEB3 ribozyme is located inside the intron of the CPEB3 gene (8), hammerhead ribozyme found in the 3’ UTR of Clec2 genes (9), hammerhead and HDV motifs associated with retrotransposable elements (20, 24). However, the OR4K15 ribozyme is located at the reverse strand of the coding region of the olfactory receptor gene OR4K15, suggesting a potentially novel functional mechanism for this ribozyme.

Previously, we performed deep mutational scanning of a self-cleaving ribozyme (CPEB3) by using error-prone PCR to generate mutants in large scale. The relative activities of these mutants can be obtained by counting the cleaved and uncleaved sequences of each mutant from high-throughput sequencing data. These activities can be utilized to analyze mutational covariation and infer accurate base-pairing structures by developing a computational tool called CODA (covariation-induced deviation of activity) (26). Here, using the same technique, we located the regions essential for the catalytic activities (the functional regions) of LINE-1 and OR4K15 ribozymes in their original sequences, which were generated from selection of randomly fragmented human genomic DNA (8). We showed that these two ribozymes shared essentially the same base-pairing structural elements but with a circular permutation. The structural elements, which shaped like the Chinese lantern, are homologous to the catalytic cores of the more complex twister sister self-cleaving ribozymes in terms of sequence and secondary structure, suggesting a more primitive origin for twister sister ribozymes. However, the homology model of twister sister-like (TS-like) ribozyme generated from the twister sister ribozyme may not be the true structure for its function due to the lack of a stem loop for its stabilization and two mismatches from twister sister in the internal loops as well as different responses to mutations.

Results

Deep mutational scanning of LINE-1 and OR4K15 ribozymes.

Our previous work indicates that deep mutational scanning can lead to co-variation signals for highly accurate inference of base-pairing structures (26). Here, the same technique was applied to the original LINE-1 ribozyme (146 nucleotides, LINE-l-ori) and OR4K15 ribozyme (140 nucleotides, OR4K15- ori). These original sequences were generated from a randomly fragmented human genomic DNA Selection based biochemical experiments (8). Thus, their functional and structural regions are unknown. It needs to be noted that there is no exact match for LINE-l-ori and OR4K15-ori in the human genome (Figure 1A). This should be related to mutations accumulated during the selection. In deep mutational scanning, high-throughput sequencing can directly measure the relative cleavage activity (RA) (See Methods) of each variant in the mutant libraries of these two ribozymes according to the ratio of cleaved to uncleaved sequence reads of the mutant, relative to the same ratio of the wild-type sequence. The relative activities of single mutations at each position indicate the sensitivity of cleavage activity to the mutations in that position. Figure 1B shows the average RA for a given sequence position for all LINE-l-ori mutants with single mutations at the position. Mutations in most positions in LINE-l-ori did not lead to large changes in RA. In other words, these sequence positions are unlikely to contribute to the specific structure required for the ribozyme to be functional. The activity of the LINE-l-ori ribozyme is only sensitive to the mutations in two short segments (54–71 and 83–99) where the RA can be reduced to nearly zero. Thus, the two terminal ends and the central region between bases 72 and 82 are not that important for LINE-1 self-cleaving activity. In Figure 1C, a similar distribution of RA was observed in OR4K15-ori. Only sequence positions inside the two short segments (70–84 and 101–116) are sensitive to mutations, which suggests that only these regions are required to form the functional RNA structure of the OR4K15 selfcleaving ribozyme. Thus, we have identified the contiguous functional regions of LINE-1 ribozymes (54–99, LINE-l-rbz) and OR4K15 ribozymes (70–116, OR4K15-rbz).

Figure 1.

The deep mutation data was further employed to search for the base pairs inside these two ribozymes by using CODA analysis. CODA employed Support Vector Regression (SVR) to establish an independent-mutation model and a naive Bayes classifier to separate bases paired from unpaired (26). Moreover, incorporating Monte-Carlo simulated annealing with an energy model and a CODA scoring term (CODA+MC) could further improve the coverage of the regions under-sampled by deep mutations. Figures 1D and 1E highlight the base pairs within the region 54–99 of LINE-l-ori (LINE-l-rbz) and the region 70–116 of OR4K15-ori (OR4K15- rbz), respectively. The CODA result for LINE-l-ori shows four base pairs (A14U34, A16U32, G17C31, C18G30) in one possible stem region, and another two base pairs (A4U43, A6A41) in the second possible stem region. Due to the lower mutational coverage, the CODA result for OR4K15-ori shows fewer base pairs than LINE-l-ori. Only two base pairs (U1A47, U3A45) in one possible stem region, and one lone base pair (A13U34) have been observed. The low mutation coverage for OR4K15-ori was due to the mutational bias (27,28) of error-prone PCR (Supplementary Figures S1, S2, S3 and S4).

A few but strong co-variation signals from CODA analysis can be expanded by MC simulated annealing as demonstrated previously (26). Several contiguous base-paired regions were detected (Figures 1D and 1E). LINE-l-rbz has two highly reliable stem regions with the length of 6 nt, including one non-Watson-Crick pair A6A41. For OR4K15-rbz, the functional region has two reliable stem regions with the lengths of 5 nt and 4 nt, respectively, including one non-Watson-Crick pair G14U33. Most of the base pairs discovered by MC simulated annealing are non-AU pairs, whose covariations are difficult to capture by error-prone PCR because of mutational biases (27, 28). In addition, these two stems are consistent with the single mutation profiles. For example, two single mutations G14A and U33C showed high RA, 1.55 and 1.82, respectively, suggesting the high possibility of a G14U33 pair at these positions, because these two separate mutations change the wobble GU pair to the Watson-Crick AU and GC pairs, respectively, and lead to a more stable structure and higher activity. If G14U33 is true, a natural extension of the stem region is G15C32 for another canonical Watson-Crick pair. This pairing region was missed by the CODA+MC result. There are some other base pairs that appeared in low probabilities after MC simulated annealing. They are considered as false positives because of low probabilities and lack of support from the deep mutational scanning results. The appearance of false positives is likely due to the imperfection of the experiment-based energy function employed in current MC simulated annealing.

Further validation of base pairing information by deep mutational scanning of LINE-1- rbz.

To further improve the signals, we employed the contiguous functional segment (54–99, LINE-l-rbz) for the second round of deep mutational scanning (Supplementary Figure S5). Performing the second round is necessary because the deep mutational scanning of the original LINE-1 ribozyme has a low 18.5% coverage of double mutations (Supplementary Table S1). This low coverage is not sufficient for accurate inference of the secondary structure in the full coverage. This second round employed a chemically synthesized doped library with a doping rate of 6%, rather than a mutant library generated from error-prone PCR biased toward the sequence positions with A/T nucleotide (Supplementary Tables S1, S2 and Figures S1, S2, S3 and S4). In addition, to amplify the signals of cleaved RNAs after in vitro transcription of the mutant library of LINE-l-rbz, we selectively capture cleaved RNAs by employing RtcB ligase and a 5’-desbiotin, 3’-phosphate modified linker because they only react with the 5’-hydroxyl termini that exist only in cleaved RNAs (29). The captured active mutants were further enriched by the streptavidin-based selection after ligation. The technology improvement leads to less biased mutations in terms of the sequence positions (Supplementary Figures S1 and S2) and mutation types (Supplementary Figures S3 and S4). More importantly, it achieves 99.3% and 99.9% coverage of single and double mutations, respectively, within the contiguous functional segment (Supplementary Table S1).

Consistent secondary structure from deep mutational scanning of LINE-l-ori and LINE-l-rbz.

We performed the CODA analysis (26) based on the relative activities of 45,925 and 72,875 mutation variants (no more than 3 mutations) obtained for the original sequence and functional region of the LINE-1 ribozyme, respectively. CODA detects base pairs by locating the pairs with large covariation-induced deviations of the activity of a double mutant from an independent single-mutation model. Figure 2A shows the distribution of relative activity (RA’, measured in the second round of mutational scanning) (See Methods) of all single mutations for the LINE-l-rbz ribozyme, which is consistent with the distribution of RA in the functional region of LINE-l-ori ribozyme (Figure 1A and Supplementary Figure S5). Moderate or high relative activity values (RA’ > 0.5) at the two outer ends (18C, 30G, C46) of the stem regions suggest that these nucleotides might not be necessary for the function of LINE-l-rbz, whereas RA’ values at 1G were affected by the transcription efficiencies of different nucleotides downstream the T7 promoter.

Figure 2.

Figure 2B shows that data from both the original sequence and functional region of the ribozyme reveal the existence of two stem regions. The latter clearly shows two stems with lengths of 6 nt and 5 nt, respectively, due to nearly 100% coverage of single and double mutations for the functional region. Moreover, it suggests a non-canonical pair 6A41A at the end of the first stem, although the signal is weak. Such a weak signal is expected because the non-canonical pairs are typically less stable than standard base pairings (30). There are some additional weak signals in the LINE-l-rbz result. Most of these signals are for the base pairs at a few sequence distances apart (|i-j|<6, local base pairs). We found that some of them are caused by a high relative activity (RA’) of the double mutant and are not likely to be true positives because the corresponding single mutations are not very disruptive (RA’ > 0.5). The consistency between base pairs inferred from deep mutational scanning of the original sequences and that of the identified functional regions confirmed the correct identification of functional regions for LINE-1 ribozyme.

Probabilistic CODA results were further refined by Monte Carlo simulated annealing (CODA+MC). The resulting base pairs (Figure 2C, upper triangle) removed mostly local false positives. However, the noncanonical AA pair was removed as well, likely because the energy function does not account for noncanonical pairs. A new 18C30G pair added is a natural extension of the second stem. Two separate consecutive base pairs (7G37C, 8C36G) with relatively low signals were also added in the CODA+MC result. These two base pairs are likely false positives because they do not appear in all the models generated from Monte Carlo simulated annealing with a low probability of 0.52. Taking all the information together leads to a confident secondary structure for the LINE-l-rbz ribozyme, which contains two stem regions (Pl, P2), two internal loops, and a stem-loop (Figure 2D left).

The consistent result between LINE-l-rbz and LINE-l-ori suggested that reliable ribozyme structures could be inferred by deep mutational scanning. This allowed us to use OR4K15-ori to directly infer the final inferred secondary structure for the functional region of OR4K15. The secondary structures for these two ribozymes (Figure 2D) are surprisingly similar to each other. Both ribozymes have two stems (Pl, P2), two internal loops and a stem-loop region. Both stem-loop regions (SL2) are insensitive of ribozyme activity to mutations (Figures 1B, 1C). The internal loop regions (LI) of LINE-l-rbz and OR4K15-rbz are nearly identical except that C8 in OR4K15-rbz is replaced by U38 in LINE-l-rbz (Figure 3A).

Figure 3.

Consensus sequence of LINE-l-rbz and OR4K15-rbz.

The consensus sequences (Figure 2E) for LINE-l-rbz and OR4K15-rbz were generated by R2R (31) based on the multiple sequence alignment of 1394 and 621 mutants with RA ≥ 0.5 from deep mutational scanning. Figure 2E overlays the consensus sequences onto the secondary structure models from the CODA+MC analysis. The self-cleavage of the two ribozymes occurs between a conserved CA dinucleotide which locates inside the longer part of the internal loops (LI) Unking P1 and P2. The internal loops are the regions with the highest conservation, which suggests their important roles in cleavage activity. Both stem regions (Pl, P2) are not very conserved because of covariations. For example, in Figure 2E (right), 1U47A can be replaced by 1C47G, and 14G33U will form 14C33G after double mutations. Covariation of all CG pairs was missing in OR4K15-rbz due to the mutational biases. The stem-loop regions show the lowest sequence identity, consistent with its insensitivity of catalytic activity to mutations in Figures 1B and 1C, providing additional support for the secondary structure models obtained here.

The secondary structure of OR4K15-rbz is a circular permutation of LINE-l-rbz selfcleaving ribozyme.

Removing the peripheral loop regions (Figures 1B and 1C) allows us to recognize that the secondary structure of OR4K15-rbz is a circular permutated version of LINE-l-rbz. As shown in Figure 3A, after a permutation, OR4K15-rbz has two conserved internal loops with one loop identical to LINE-l-rbz and the other loop differed from LINE-l-rbz by one base. The only mismatch U38C in LI has the RA’ of 0.6, suggesting that the mismatch is not disruptive to the functional structure of the ribozyme. As the simplest ribozymes reported so far, it is important to know if these two ribozymes share the same motifs with the known ribozyme families.

Here we used pattern-based similarity search (RNAbob, http://eddylab.org/software.html) to search the structural patterns against the known ribozyme families in Rfam database (Supplementary Figure S6). The structural patterns of two ribozymes were derived from our deep mutational scanning results. We obtained three hits with identical motifs all from the twister sister ribozyme class (Supplementary Table S3). Figure 3A further shows the comparison with two previously published twister sister structures (32, 33). The two internal loops of LINE-l-rbz with 5 and 6 nucleotides, respectively, differ from the catalytic internal loops of the four-way junctional twister sister ribozyme only by one base in each loop (Figure 3A) (32). Given the identical lantern-shaped regions of the LINE-l-rbz and OR4K15-rbz ribozyme, we named them twister sister-like (TS-like) ribozymes. Thus, it is possible to use the structure of twister sister ribozyme to build a homology model.

Homology modelling of the TS-like ribozymes.

To obtain the structure model of TS-like ribozymes, we used template-based homology modelling with the twister sister ribozyme structure as the template. Figure 3B shows the homology modelled structure of bimolecular LINE-1-core (internal loops plus stems) built from the more identical four-way junctional twister sister ribozyme (PDB ID: 5Y87). The structure of the four-way junctional twister sister ribozyme revealed the internal loops LI as the catalytic region involving a guanine-scissile phosphate interaction (G5-C62-A63), continuous stacking interactions, additional pairings and hydrated divalent Mg²⁺ ions. Figure 3B, 3C and Supplementary Figure S7 show that the LINE-1-core model can be built as the twister sister ribozyme in the internal loops LI, with Al 1 at the cleavage site directed inwards. As shown in Figures 3C and Supplementary Figure S8, stem Pl of LINE-l-rbz is extended by the Watson-Crick U9A39 as a part of the G7(U9A39) base triple and Watson-Crick C8G40. Meanwhile, stem P2 of LINE-l-rbz is extended through trans non-canonical A12G36 and trans sugar edge-Hoogsteen Al 1C37. In addition to the triple base interaction in LINE-l-rbz, G7 is also H-bonded to both non-bridging O atoms of the phosphate connecting C37 and U38 (Figure 4A). This, along with the H-bond interaction between proR O atom of the phosphate and A39 N6, generates a stable square array of H-bonds.

Figure 4.

We also modelled the structure of OR4K15-core in the same way as LINE-1-core. As shown in Supplementary Figure S9A and S9B, the model structure of LI in OR4K15-core is similar to LI in the LINE-1-core model, with stems P1 and P2 extended by additional base-pairing interactions in the tertiary structure of OR4K15-core. The base triple interaction of G37(A9U39) also exists in OR4K15-core, with the difference that G37 is H-bonded to A9 N6, rather than U39 (Supplementary Figure S9C). Thus, the square array of H-bonds could not be observed as there was no direct interaction between A9 and the phosphate connecting C7 and C8. This modelling result suggests that the bond between A9 (A39 in LINE-1-core) and the phosphate may not be important, consistent with the fact that the replacement of AU pair by a non-canonical UU pair did not disrupt cleavage activity (RA’ =0.8 for A39U in LINE-l-rbz).

Like the two twister sister ribozyme structures (32, 33), the nucleotide U38 in the LINE-1- core model (or C8 in the OR4K15-core model) is extruded from the helix (Figure 3C, Supplementary Figure S9A). The corresponding nucleotides, A8 in three-way junctional twister sister ribozyme (33) and A7 in four-way junctional twister sister ribozyme (32) are involved in the tertiary contact with the stem loop SL4 (Figure 3A). Lacking the extra stem loops in twister sister ribozyme, the key contacts for U38 of LINE-1-core (C8 in OR4K15-core) are the H-bond interaction with Al 1 N6 (A41 N6 in OR4K15-core), and the interactions between the phosphate and G7 (G37 in OR4K15-core). No involvement of the nucleobase part (U38 in LINE-l-core, C8 in OR4K15-core) in the modelled structures is compatible to the mismatches U38C and U38A found in the comparison, as changes of the nucleobase part did not have a big influence on the H-bonded interactions involving U38 (C8 in OR4K15-core) (Figure 3A). Thus, it is not entirely clear whether the model structure of TS-like ribozymes built from twister sister will remain stable in the absence of any tertiary interactions with the stem loop SL2 or it will deviate significantly from the twister sister.

Differences in mutational responses in the Lis between the TS-like and the twister sister ribozymes.

Considering the high similarity of the internal loops, we further investigated the mutational effects on the internal loop Lis. Most nucleotides in the LI region of LINE-l-rbz are relatively conserved, as single mutations showed RA’ < 0.2 according to deep mutational scanning results (Figure 2A). The conservation is consistent with the stacking and hydrogenbonding interactions in the catalytic region. By comparison, mutations of C62 (C54 in three-way junctional twister sister ribozyme) at the cleavage site did not lead a major change on the cleavage activity in previous studies (32, 33). This can be compared to the fact that single mutations of the corresponding nucleotide in LINE-l-rbz (CIO) had relative activities as low as 0.07 in the deep mutational scanning result. To further confirm the above result, we performed cleavage assays on several mutations at the cleavage site of LINE-l-core. As shown in Figure 4C, mutations on CIO showed either partial (C10U, C10G) or complete loss (C10A) in cleavage activity. While the corresponding mutations in the four-way junctional twister sister ribozyme showed pronounced (C62U, C62A) or somewhat reduced (C62G) cleavage activity (32, 33). More interestingly, A11U in LINE-l-core showed partial cleavage activity, whereas the corresponding mutation A63U in the four-way junctional twister sister ribozyme showed complete loss of cleavage activity (32, 33). These different mutation effects indicate that LINE-l-core may not adopt exactly the same structure and catalytic mechanism as twister sister ribozyme. In other words, TS-like ribozymes may not be simply a minimal version of twister sister ribozymes.

We further examined the interactions around the cleavage site C 10-Al 1 in LINE-l-rbz (C40- A41 in OR4K15-core). The corresponding nucleotides in twister sister ribozymes (C62-A63 in four-way junctional twister sister ribozyme, C54-A55 in three-way junctional twister sister ribozyme) adopt either a splayed-apart or a base-stacked conformation, with the scissile phosphate H-bonded to the G5 (four-way junctional twister sister ribozyme) or C7 (three-way junctional twister sister ribozyme). In our modelled structures (Figure 4B, Supplementary Figure S9D), the bases at the cleavage site CIO-All in LINE-l-core (C40-A41 in OR4K15- core) are splayed away, with All (A41 in OR4K15-core) directed inwards and CIO (C40 in OR4K15-core) directed outwards. In addition, H-bonded interactions between CIO (C40 in OR4K15-core) and the two nucleotides 3’ of CIO could be found (Figure 4B, Supplementary Figure S9D). These structure differences may partially explain why mutations of CIO were detrimental to the catalysis in TS-like ribozymes. Moreover, we did not observe interactions between CIO (C40 in OR4K15-core) and G36/C37 (G6/C7 in OR4K15-core) although they were important to anchor the cytosine in twister sister ribozymes. This may be due to replacement of the non-canonical pair G5U64 (G6U57 in three-way junctional twister sister ribozyme) by A12G36 (G6A42 in OR4K15-core) in TS-like ribozymes. Thus, the stability of the current model structure based on the template from twister sister is uncertain.

Activity confinnation and biochemical analysis of LINE-l-core and OR4K15-core.

To confirm the cleavage activity of the core region (without the stem-loop linker), we obtained the segments (54–71 for LINE-1, 101–116 for OR4K15) as the substrate strands with the cleavage site and the segments (83–99 for LINE-1, 70–84 for OR4K15) as the enzyme strands (Figure 5A), as they made up the LINE-l-core and OR4K15-core. When the substrate strand was mixed with the enzyme strand, most of the substrate strand was cleaved in the presence of Mg²⁺ (Figure 5B). When Mg²⁺ was not included in the reaction, no cleavage could be observed even after 24h’s incubation (Supplementary Figure S10). This result confirms that the RNA cleavage in TS-like ribozymes is a catalytic reaction, accelerated by Mg²⁺. Moreover, the result confirmed that the stem loop SL2 regions in LINE-l-rbz and OR4K15-rbz (Figure 2) did not participate in the catalytic activity.

Figure 5.

The core constructs of the two TS-like ribozymes contain 35 and 31 nucleotides, respectively. They are the simplest and the second simplest self-cleaving ribozymes reported so far. In all previously reported self-cleaving ribozyme families, the self-cleaving reaction occurs through an internal phosphoester transfer mechanism, in which the 2’-hydroxyl group of the -1 (relative to the cleavage site) nucleotide attacks the adjacent phosphorus resulting in the release of the 5’ oxygen of+1 (relative to the cleavage site) nucleotide (2, 34–37). We confirmed that both LINE-l-core and OR4K15-core employed the same cleavage mechanism because the analogous RNAs that lacked the 2’ oxygen atom in the -1 nucleotide (dClO for LINE-l-core, dC40 for OR4K15-core) were unable to cleave (Figure 5B). We also investigated whether TS-like ribozymes can cleave when Mg²⁺ was replaced by Co(NH₃)₆³⁺ in the reaction. Co(NH₃)₆³⁺ is isosteric with Mg(H₂O)₆²⁺, but the divalent cations cannot directly participate in catalysis as the amino ligands cannot readily dissociate. In Figure 5C, a total loss of cleavage activity in Co(NH₃)₆³⁺ in the absence of Mg²⁺ indicates that the TS-like ribozymes require the binding of divalent cations not only for structure folding, but also for catalysis.

As shown in Figure 5D, we further examined the dependence of the metal ions on TS-like ribozymes’ cleavage activity. At a concentration of ImM, ribozyme cleavage can be observed with Mg²⁺, Mn²⁺, Co²⁺ and Zn²⁺ but little or none with Ca²⁺, Cu²⁺, Ba²⁺, Ni²⁺, Na⁺, K⁺, Li⁺, Cs⁺ or Rb⁺, indicating that direct participation of specific hydrated divalent metal ions is required for self-cleavage. More interestingly, we found that the TS-like ribozymes have an equivalent or even higher cleavage ratio with Mn²⁺ than Mg²⁺. For LINE-l-core, the cleaved fractions were ~57% for Mg²⁺, ~74% for Mn²⁺. For OR4K15-core, the cleaved fractions were ~9% for Mg²⁺, ~79% for Mn²⁺. We further characterized the cleavage rates of these two ribozymes under different concentrations of Mg²⁺ and Mn²⁺ by using a fluorescence resonance energy transfer (FRET)-based method, as shown in Supplementary Figure S11 and Figure S12. We used the first-order rate constant (k_obs) to represent the efficiency of the cleavage reaction. The k_obs of LINE-l-core under single-turnover condition was ~0.05 min¹ when measured in lOrnM MgCl₂ and 100mM KCl at pH 7.5 (Figure S13). Only a slightly lower value of k_obs (~0.03 min^-1) was observed for LINE-l-rbz (Figure S14). This confirms that the stem loop region SL2 does not contribute much to the cleavage activity of the TS-like ribozymes. Cleavage activities of the two TS-like ribozymes were highly dependent on the concentration of divalent metal ions. The steep increase in rate constants of two TS-like ribozymes plateaued at Mg²⁺ concentrations above 100mM, while for Mn²⁺ it plateaued at concentrations of 10mM. Consistent with the PAGE result, LINE-l-core showed higher cleavage rates than OR4K15- core when they were incubated with Mg²⁺. However, the cleavage rates of OR4K15-core were only lower than LINE-l-core when the concentrations of Mn²⁺ were less than 20mM, but higher, otherwise. Thus, there may exist important difference between OR4K15-core and LINE-l-core to explain this different bias toward ions. A detailed explanation for this difference may require high-resolution structure determination of the TS-like ribozymes.

Homology search of two TS-like ribozymes.

To locate close homologs of the two TS-like ribozymes, we performed cmsearch based on a covariance model (38) built on the sequence and secondary structural profiles. In the human genome, we got 1154 and 4 homolog sequences for LINE-l-rbz and OR4K15-rbz, respectively. For OR4K15-rbz, there was an exact match located at the reverse strand of the exon of OR4K15 gene (Figure 6A). The other 3 homologs of OR4K15-rbz belongs to the same olfactory receptor family 4 subfamily K (Figure 6C). However, there was no exact match for LINE-l-rbz (Figure 6A). Interestingly, a total of 1154 LINE-l-rbz homologs were mapped to the LINE-1 retrotransposon according to the RepeatMasker (http://www.repeatmasker.org) annotation. Figure 6B showed the distribution of LINE-l-rbz homologs in different LINE-1 subfamilies in the human genome. Only three subfamilies L1PA7, L1PA8 and L1P3 (L1PA7-9) can be considered as abundant with LINE-l-rbz homologs (>100 homologs per family). The consensus sequences of all homologs obtained are shown in Figure 6D. In order to investigate the self-cleavage activity of these homologs, we mainly focused on the mismatches in the more conserved internal loops. The major differences between the 5 consensus sequences are the mismatches in the first internal loop. The widespread A12C substitution can be found in majority of LINE-l-rbz homologs, this substitution leads to a one-base pair extension of the second stem (P2) but almost no activity (RA’: 0.03) based on our deep mutational scanning result. Then we selected 3 homologs without A12C substitution for LINE-l-rbz for in vitro cleavage assay (Figure 6E). But we didn’t observe significant cleavage activity, this might be caused by GU substitutions in the stem region. For 3 homologs of OR4K15-rbz, we only found one homolog of OR4K15 with pronounced self-cleavage activity (Figure 6F). In addition, we performed similar bioinformatic search of the TS-like ribozymes in other primate genomes. Similarly, the majority (15 out of 18) of primate genomes have a large number of LINE-1 homologs (>500) and the remaining three have essentially none. However, there was no exact match. Only one homolog has a single mutation (U38C) in the genome assembly of Gibbon (Figure S15). The majority of these homologs have 3 or more mismatches (Figure S15). For OR4K15-rbz, all representative primate genomes contain at least one exact match of the OR4K15-rbz sequence.

Figure 6.

Discussion

The wide distribution of RNA self-cleaving activity in different species suggest that selfcleaving ribozymes play an important role in living organisms. Most self-cleaving ribozymes (39) are in non-coding regions of the hosting genome. For example, the CPEB3 ribozyme in the intron of the CPEB3 gene (8), hammerhead ribozyme in the 3’ UTR of Clec2 genes (9), hammerhead and HDV motifs associated with retrotransposons (11, 40). LINE-1 belongs to the family of retrotransposons, they make up about 17% of the human genome, with over 500,000 copies (41, 42). The retrotransposition process of LINE-1 increases the copy number of repeats and generates insertional mutations in the genome, which contributes to the genetic novelties as well as genomic instability. However, we found that the LINE-l-rbz sequence shared at least two mismatches with the ancient human LINE-1 families, that led to the loss of self-cleavage activity. Considering the location of the LINE-1 ribozyme inside the 5’ UTR, it is suggested that this ribozyme might be related to the transcriptional regulation of the LINE-1 retrotransposon during evolution, like the HDV-like ribozymes. By comparison, the OR4K15 ribozyme is located at the reverse strand of the coding region of the olfactory receptor gene OR4K15, suggesting a potentially novel functional mechanism for this ribozyme as antisense transcription has been shown to be ubiquitous in mammals and may affect the gene regulation (43). However, the progress on understanding the functional roles of LINE-1 and OR4K15 self-cleaving ribozymes was slow because of lacking structural clues.

The centerpiece of structural clues is the RNA base-pairing structure resulted from the complex interplay of secondary and tertiary interactions. Previously, we developed a technique (26) for inferring the base-pairing information of the self-cleaving ribozymes by deep mutational scanning. In this paper, we applied this technique to the LINE-1 and OR4K15 ribozymes, to obtain the functional regions and the base-pairing information. As both LINE-1 and OR4K15 ribozymes were discovered from the randomly fragmented human genomic DNA selection based biochemical experiments, the original full-length (WT) sequences of the ribozymes do not necessarily represent the functional region of the RNAs. However, the functional regions can cleave when they were transcribed along with the upstream and downstream sequences, and, thus, they are more likely to be active in cellular environment.

Indeed, we found that the functional regions for both ribozymes are very short (35 and 31 nucleotides for LINE-1 and OR4K15, respectively) and surprisingly in the same class in a circular permutated form. We are confident about the secondary structures obtained because of the consistency between the results from the deep mutational study of the original sequence and the functional region of the LINE-1 ribozyme. The obtained secondary structure is further confirmed by the discovery of the functional regions of the two TS-like ribozymes in the circular permutated form from one another (Figure 3A). From the mutational information of LINE-l-rbz, the only mismatch (U38C) from OR4K15 in the catalytic region has the RA’ of 0.6. While the corresponding mismatch mutation C8U in OR4K15 has the RA of 1.54, or 0.65 for U8C. This consistence in activity reduction for the U-→C mutation further confirmed that these two ribozymes belong to the same class in a circular permutated form. Currently, three classes of self-cleaving ribozymes, the hammerhead (7, 44—47), twister (34) and the hairpin (48) ribozymes are found in different species in multiple circular permutated forms. The discovery of permutated TS-like ribozymes in the human genome for the first time suggests that we can find additional examples of circular permutated TS-like ribozymes in other species.

More interestingly, we found that the internal loop regions of the TS-like ribozymes share a high similarity with the known catalytic region of the twister sister ribozyme (2 mismatches only, Figure 3A). The conserved cleavage sites suggest a similar structure and mechanism as shown in Figure 3A. However, the maximum observed rate constant (k_obs) for twister sister ribozyme was ~5 min^-1 (37), compared to ~0.ll ± 0.01min¹ for the TS-like ribozymes. This rate difference is not mainly caused by mismatches in the catalytic core between the twister sister and LINE-1/OR4K15 ribozymes. In fact, deep mutational scanning results indicate that one mismatch (U38C) in the catalytic region of OR4K15-rbz (Figure 3A) has the RA’ of 0.6. Two mismatches (A12U and U38A) in the catalytic region of the twister sister ribozyme (Figure 3A) have the RA’ of 0.28 and 0.47, respectively, whereas the coexistence of these two mismatches only has the RA’ of 0.09. Thus, a more sophisticated structure along with long-range interactions involving the SL4 region in the twister sister ribozyme must have helped to stabilize the catalytic region for the improved catalytic activity. Similarly, previous studies have demonstrated that peripheral regions of hammerhead (49), hairpin (50) and HDV (51,52) ribozymes could greatly increase their self-cleavage activity. Given the importance of the peripheral regions, absence of this tertiary interaction in the TS-like ribozyme may not be able to fully stabilize the structural form generated from homology modelling.

Indeed, we found some differences between the responses to mutations in the catalytic regions around the cleavage sites of TS-like ribozymes and those of twister sister ribozymes. Although there are pronounced differences around the cleavage sites in previous twister sister ribozyme structures, both cytosines at the cleavage sites were found to be insensitive to mutations (32, 33). However, unlike twister sister ribozymes, the cytosine at the cleavage site of TS-like ribozyme was sensitive to mutations in our mutation assays. According to our modelled structures, the H-bonding interactions associated with the cleavage site C10-A11 (C40-A41 in OR4K15-core) in TS-like ribozymes show that the cytosine (CIO in LINE-l-core, C40 in OR4K15) is relatively conserved, which is consistent with the mutation assay results. In the modelled structure, we did not observe direct contacts helped to stabilize the scissile phosphate and the in-line alignments required for the catalytic step. This may suggest that conformational changes may be required to generate the active conformation prior to the catalytic step as in the twister sister ribozyme structures. Or this is due to inexact structures generated from homology modelling, because Mg²⁺ was not included during homology modelling. Our experimental results (Figure 5B-D) showed that Mg²⁺ is essential for both structure folding and catalysis. Absence of Mg²⁺ during modelling might lead to the conformational change of the tertiary structures. Thus, further studies are clearly needed to illustrate possible structural and mechanistic difference between TS-like ribozymes and twister sister ribozyme.

According to the bioinformatic analysis result, there are some TS-like ribozymes (one LINE-1-rbz homolog in the Gibbon genome, and some OR4K15-rbz homologs) with in vitro cleavage activity in primate genomes. Unlike the more conserved CPEB3 ribozyme which has a clear function, the function of the TS-like ribozymes is not clear, as they are not conserved, belong to the pseudogene or located at the reverse strand. While the in vivo activities and functions of these two circular-permutated TS-like ribozymes require further studies, their discovery may have important therapeutics and biotech implications. First, the TS-like ribozymes have the simplest secondary structure with two stems and two internal loops, compared to all naturally occurring self-cleaving ribozymes. Previous studies also have revealed some minimized forms of self-cleaving ribozymes, including hammerhead (19, 53) and HDV-like (54) ribozymes. However, when comparing the conserved segments, they (>= 36 nt) are not as short as the TS-like ribozymes (31 nt) found here. Second, this TS-like ribozyme after removing the stem-loop region only requires a 15-nucleotide enzyme strand and a 16-nucleotide substrate strand for its function. Thus, it can be easily modified into a trans-cleaving ribozyme, which was found useful for bioengineering and RNA-targeting therapeutics (55). A recent intracellular selectionbased study has shown that hammerhead ribozyme can be engineered into a trans-acting ribozyme with consistent gene knockdown ability on different targeted mRNA either in prokaryotic or eukaryotic cells (55). The TS-like ribozyme, the simplest one known, may have an unprecedented advantage, compared to existing trans-cleaving ribozymes.

What is also important is that this work further confirms the usefulness of deep mutational scanning and high-throughput sequencing for RNA structural inference by CODA analysis (26, 56). The method is not limited to self-cleaving ribozymes, if combined with other functional or phenotypic selection techniques. Moreover, the structural and full mutational information provided by this method could be utilized to discover additional functional RNAs and thus help us to investigate the hidden side of non-coding RNAs.

Materials and Methods

All the oligonucleotides listed in Supplementary Table S4, except for the doped library, were purchased from IDT (Integrated DNA technologies) and Genscript. The doped mutant library Rz_LINE-l_doped with a doping rate of 6%, was purchased from the Keck Oligo Synthesis Resource at Yale University. All the high-throughput sequencing experiments were performed on the Illumina HiSeq X platform by Novogene Technology Co., Ltd.

Deep mutational scanning experiments on original ribozyme sequences

The workflow for two deep mutational scanning experiments is illustrated in Supplementary Figure S13. The deep mutational scanning of the original sequences of LINE-1 ribozyme (labelled as LINE-l-ori) and OR4K15 ribozyme (labelled as OR4K15-ori) followed the same protocol as our previous work (26). Briefly, the mutant library was generated from 3-rounds of error-prone PCR, barcoded using primer Bar F and Bar R (Supplementary Table S4), and then diluted. After that, the transcribed RNA of the barcoded mutant library was reverse transcribed with RT ml3f adpl and template switching oligo TSO (Supplementary Table S4). The total cDNA was used to generate the RNA-seq library with primer P5R1 adpl and P7R2 adp2 (Supplementary Table S4). The DNA-seq libraries were generated by amplifying the ribozyme mutant library with primer P5Rl_ml3f and P7R2_t7p (Supplementary Table S4). Both DNA-seq and RNA-seq libraries were sequenced on an Illumina HiSeq X sequencer with 25% PhiX control by Novogene Technology Co., Ltd. We counted the number of reads of the cleaved and uncleaved portions of each variant by mapping their respective barcodes. The relative activity (RA) of each variant is calculated by using the equation RA(var) = N_cleaved (var) N_total (wt)/N_total (var) N_cleaved (wt).

Covariation-induced deviation of activity (CODA) analysis

We used CODA (26) to analyze the RA of LINE-l-ori and OR4K15-ori. For Monte Carlo (MC) simulated annealing, we used the same weighting factor of 2 for both datasets as employed previously (26).

Deep mutational scanning and CODA analysis on the functional region of LINE-1 ribozyme

We identified the functional region of the LINE-1 ribozyme (labelled as LINE-l-rbz) by removing the terminal sequence positions whose mutations do not affect cleavage activities. Then, the second round of deep mutational scanning experiment was applied to the LINE-l-rbz. Here, we used a doped mutation library instead of a library generated by error-prone PCR. The doped mutant library Rz_LINE-l_doped (Supplementary Table S4) was amplified by using primer T7prom and M13F (Supplementary Table S4). The PCR product was then gel-purified, quantified, and diluted. Approximately 5 × 10⁵ doped DNA molecules were amplified by T7prom and M13F primers (Supplementary Table S4) to produce enough DNA template for in vitro transcription. The transcribed and purified RNAs (10 pmol) of the mutant library were mixed with 100 pmol rM13R_5desBio_3P (Supplementary Table S4), 2 μl 10X RtcB reaction buffer (New England Biolabs), 2 μl 1mM GTP, 2 μl 10mM MnCl₂, 0.2 μl of Murine RNase Inhibitor (40 U/μl, New England Biolabs) and 1 μl RtcB RNA ligase (New England Biolabs) to a total volume of 20 μl. The reaction mixture was incubated at 37 °C for 1 h, then purified using RNA Clean & Concentrator-5 kit (Zymo Research). The purified products were mixed with 2 μl of 10 μl M RT_ml3f_ adp1 (Supplementary Table S4) and 1 μl of 10mM dNTPs in a volume of 8 μl, and then heated to 65 °C for 5min and placed on ice. Reverse transcription was initiated by adding 4 μl of 5x ProtoScript II Buffer (New England Biolabs), 2 μl of 0.1 M DTT, 0.2 μl of Murine RNase Inhibitor (40 U/pl, New England Biolabs) and 1 μl ProtoScript II RT (200 U/μl, New England Biolabs) to a total volume of 20 μl. The reaction mixture was incubated at 42°C for 1 h, and then purified by Sera-Mag Magnetic Streptavidin-coated particles (Thermo Scientific).

The purified products after streptavidin-based selection were amplified by PCR to construct the RNA-seq library with primer P7R2_ml3r and P5Rl_ml3f (Supplementary Table S4). The DNA template used for transcription was used as the template for constructing the DNA-seq library with primer P5Rl_ml3f and P7R2_t7p (Supplementary Table S4). The DNA-seq and RNA-seq libraries were also sequenced on an Illumina HiSeq X sequencer with 25% PhiX control by Novogene Technology Co., Ltd. Raw paired-end sequencing reads were filtered and then merged by using the bioinformatic tool PEAR (57) to generate the high-quality merged reads. We counted the read number of each variant in DNA-seq and RNA-seq by mapping the barcode. Then, we estimated the relative activity (RA’) of each variant by using the equation RA’(var) = N_RNAseq(var)N_DNAseq(wt)/N_DNAseq(var)N_RNAseq(wt).

This relative activity (RA’) measure is different from the relative activity (RA) when both cleaved and uncleaved sequences were obtained in RNA-seq in the first-round deep mutational scanning. However, it should be a good estimate of RA because a mutant with a higher activity (cleavage rate) should have a higher possibility to be captured by RtcB ligase. We confirmed RA’ as a good estimate of RA by comparing shared single mutations from LINE-l-ori and LINE-l-rbz and obtained a strong Pearson correlation coefficient of 0.66 (Spearman correlation coefficient of 0.769). Moreover, both RA and RA’ (18–30) show that the central region of LINE-l-rbz is insensitive to mutations (Supplementary Figure S5). Moreover, secondary structures derived from RA’ and RA are consistent with each other (see below), confirming the correct identification of structural and functional regions of LINE-1 ribozyme.

We used the same method CODA (26) to analyze the RA’ of LINE-l-rbz. For Monte Carlo (MC) simulated annealing, we used the same weighting factor of 2.

PAGE-based cleavage assays

Both ribozyme sequences were separated into two portions: the substrate part (LI_S_5F3T and OR_S_5F3T, Supplementary Table S4) and the enzymatic part (LI_E and OR_E, Supplementary Table S4). LI_S_5F3T and OR_S_5F3T were synthesized and HPLC purified by Genscript. Cytosine CIO and C40 were substituted by deoxy cytosine in LI_S_dC10 and OR_S_dC40 (Supplementary Table S4) during synthesis.

In a 20 μl reaction system, 40 pmol LI_S_5F3T or OR_S_5F3T was mixed with 2 μl 0.2 M Tris-HCl and 2 μl 10 mM metal ion stock solution unless noted otherwise. 60 pmol LI_E/OR_E (Supplementary Table S4) was added to initiate the reaction. Reactions were incubated at 37°C for 1 h, and then stopped by adding an equal volume of stop solution (95% formamide, 0.02% SDS, 0.02% bromophenol blue, 0.01% Xylene Cyanol, ImM EDTA). The reaction products were then separated by denaturing (8 M urea) 15% PAGE and visualized by SYBR Gold (Thermo Fisher) staining.

For bimolecular mutants of LINE-l-rbz and OR4K15-rbz, we separated the ribozymes into the substrate part (5’ 6-FAM labelled) and the enzymatic part in a similar way. In a 20 μl reaction system, 40 pmol substrate RNA was mixed with 2 μl 0.2 M Tris-HCl and 2 μl 10 mM Mg²⁺ stock solution. 60 pmol enzyme RNA was added to initiate the reaction. Reactions were incubated at 37°C for different time intervals, and then stopped by adding an equal volume of stop solution (95% formamide, 0.02% SDS, 0.02% bromophenol blue, 0.01% Xylene Cyanol, ImM EDTA). The reaction products were then separated by denaturing (8 M urea) 15% PAGE.

For single-stranded mutants of LINE-l-rbz and OR4K15-rbz, we incubated 60 pmol RNA in a 20 μl reaction system. Mg²⁺ stock solution was added to initiate the reaction. The final concentrations for the reaction solution were 20 mM Tris-HCl buffer (pH 7.5), 100 mM KC1, and 10 mM MgCh. Reactions were incubated at 37°C for different time intervals, and then stopped by adding an equal volume of stop solution (90% formamide, 0.02% SDS, 0.02% bromophenol blue, 0.01% Xylene Cyanol, 50mM EDTA). The reaction products were then separated by denaturing (8 M urea) 12% PAGE.

Fluorescence-based kinetic analysis

The fluorescence-based kinetic analysis of the two TS-like ribozymes was similar to our previous work on the CPEB3 ribozyme. The substrate RNAs LI_S_5F3T and OR_S_5F3T were synthesized with 5’ FAM as the fluorophore and 3’ TAMRA as the quencher. Cleavage of the substrate part by the enzymatic part will relieve the quenching and therefore generate a fluorescence signal.

For one 100 μl reaction system, 100 pmol enzyme RNA was mixed with 20 pmol substrate RNA unless noted otherwise. In the ion-dependent assay, a final concentration of 20mM Tris-HCl solution (pH 7.5) and lOOmM KC1 was used. Different concentrations of MgCh or MnCh solutions were employed to initiate the reactions.

The fluorescence (Ex 488 nm/Em 520 nm) was measured at 37°C for 3 hours using BioTek Synergy Hl, and then normalized by the control (no divalent ions) at t = 0. The first-order rate constant of ribozyme cleavage k_obs was calculated in a similar way as previously described (26). The kinetics data were fitted to , where F is the fluorescence intensity, t is the time, A is the fluorescence at completion, and B is the amplitude of the observable phase. Each data point in Supplementary Figure S11 is the average of k_obs of 3 independent reactions.

3D structure modelling of the TS-like ribozymes

We matched the base-pairs of the TS-like ribozymes to the most similar twister sister ribozyme (PDB: 5Y87) by map_align (58). And then we modelled the RNA structure by replacing, inserting, and deleting the nucleotides with household RNA modelling python script. Finally, we use RNA-BRiQ (59) developed by our group to optimize the structure model. The source code can be accessed at https://gitee.com/hongxu66/coda-map.

Secondary structure-based search of LINE-l-rbz and OR4K15-rbz

All representative primate genome assemblies in fasta format and the corresponding RepeatMasker (http://www.repeatmasker.org) annotation files were downloaded from the UCSC genomic website (http://genome.ucsc.edu). We used cmbuild from Infernal (60, 61) to build the covariance model (38) of LINE-l-rbz and OR4K15-rbz with the wild type sequence and predicted secondary structure from deep mutational scanning experiment as input. Afterwards, we used the covariance model to search against the latest release of different representative primate genome assemblies by using cmsearch from Infernal. The hits with E-value <1.0 were then used for downstream analysis.

Data Availability

Illumina sequencing data for the LINE-l-ori, OR4K15-ori and LINE-l-rbz were submitted to the NCBI Sequence Read Archive (SRA) under SRA accession number PRJNA662002 (https://www.ncbi.nlm.nih.gov/sra/PRJNA662002).

The homology modelling structures of LINE-l-core and OR4K15-core are available in ModelArchive (https://modelarchive.org/doi/10.5452/ma-chp4j for LINE-l-core; https://modelarchive.org/doi/10.5452/ma-bqqf8_for_OR4K15-core).

Competing financial interests

The authors declare no competing financial interests.

Additional Declarations:

The authors declare no competing interests.

Acknowledgements

This work was funded by China Postdoctoral Science Foundation (fund no. 2021M702295 to Z.Z.) and a grant from National Natural Science Foundation of China (Grant #22350710182) to Y.Z. This work was also supported in part by Australian Research Council DP180102060 and DP210101875 and by National Health and Medical Research Council (1121629) of Australia to Y.Z. We would like to acknowledge the use of the High-Performance Computing Cluster ‘Gowonda’ at Griffith University and the supercomputing facility at the Shenzhen Bay Laboratory to complete this research. This research/project has also been undertaken with the aid of the research cloud resources provided by the Queensland Cyber Infrastructure Foundation (QCIF) and Australian Research Data Commons (ARDC). We also gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan V GPU used for this research.

Author Contributions

JZ conceived of the study. JZ and JW designed experimental studies and wrote the manuscript; ZZ and JZ conducted experimental studies and wrote the manuscript; PX and XH developed computational methods. ZZ performed computational and bioinformatics analysis. YZ participated in the initial design, assisted in data analysis, and drafted the manuscript; and all authors read, contributed to the discussion, and approved the final manuscript.

Supplementary materials

Table S1.

Table S2.

Table S3.

Table S4.

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Figure S8.

Figure S9.

Figure S10.

Figure S11.

Figure S12.

Figure S13.

Figure S14.

Figure S15.

Figure S16.