Introduction

Ribozymes (catalytic RNAs) are thought to dominate 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 human genome. They include a HDV-like cytoplasmic polyadenylation element-binding protein 3 (CPEB3) ribozyme within the transcript of a single-copy gene (12), 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. All these ribozymes were discovered by using selection-based biochemical experiments. 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. However, little is known about the above human self-cleaving ribozymes except the CPEB3 ribozyme whose single nucleotide polymorphism (SNP) was associated with an enhanced self-cleavage activity along with a poorer episodic memory (14).

Previous studies on self-cleaving ribozymes found in retrotransposons suggested their importance in the overall size, structure, and function of different genomes. Most self-cleaving ribozymes found in retrotransposons belong to two widespread ribozyme families: HDV-like ribozymes in R2 elements (10) and L1Tc retrotransposon (11) and hammerhead ribozymes (HHR) in short interspersed nuclear elements (SINEs) of Schistosomes (15), Penelope-like elements (PLEs) (16, 17) and retrozymes (18). 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 (19, 20), Bombyx mori (21), arthropods, nematodes, birds and tunicates (22, 23). 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 human genome may also play an active role in co-transcriptional processing during evolution.

Typically, self-cleaving ribozymes are in noncoding regions of the hosting genome. For example, the CPEB3 ribozyme is located inside the intron of the CPEB3 gene (14), hammerhead ribozyme found in the 3′ UTR of Clec2 genes (16), hammerhead and HDV motifs associated with retrotransposable elements (17, 18). 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 (DMS) 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) (24). Here, using the same technique, we located the functional and base-pairing regions of LINE-1 and OR4K15 ribozymes in their original sequences, which were generated from a randomly fragmented human genomic DNA Selection based biochemical experiments (12). We showed that these two ribozymes shared essentially the same base-pairing structural elements but with a circular permutation. The structural elements, called lantern ribozymes due to their shape, are homologous to the catalytic cores of the more complex twister sister self-cleaving ribozymes in term of sequence and secondary structure, suggesting a more primitive origin for lantern ribozymes. However, the homology model 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 response 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 (24). Here, the same technique was applied to the original LINE-1 ribozyme (146 nucleotides, LINE-1-ori) and OR4K15 ribozyme (140 nucleotides, OR4K15-ori). These original sequences were generated from a randomly fragmented human genomic DNA Selection based biochemical experiments (12). Thus, their functional and structural regions are unknown. In this method, high-throughput sequencing can directly measure the relative cleavage activity (RA) 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 1A shows the average RA for a given sequence position for all LINE-1-ori mutants with single mutations at the position. Mutations in most positions in LINE-1-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-1-ori ribozyme is only sensitive to the mutations in two short segments (54-71 and 83-99) where the RA can be reduced to the 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 1B, 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 self-cleaving ribozyme. Thus, we have identified the contiguous functional regions of LINE-1 ribozymes (54-99, LINE-1-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. Figures 1C and 1D highlight the base pairs within the region 54-99 of LINE-1-ori (LINE-1-rbz) and the region 70-116 of OR4K15-ori (OR4K15-rbz), respectively. The CODA result for LINE-1-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-1-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 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 (24). Several contiguous base-paired regions were detected (Figures 1C and 1D). LINE-1-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. 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 DMS of LINE-1-rbz. To further improve the signals, we employed the contiguous functional segment (54-99, LINE-1-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-1-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 (25). 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 DMS of LINE-1-ori and LINE-1-rbz. We performed the CODA analysis (24) based on the relative activities of 45,925 and 72,875 mutation variants 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) of all single mutations for the LINE-1-rbz ribozyme, which is consistent with the distribution of RA in the functional region of LINE-1-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-1-rbz, whereas RA′ values at 1G were affected by the transcription efficiencies of different nucleotides downstream the T7 promoter.

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 (26). There are some additional weak signals in the LINE-1-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 and OR4K15.

Figure 2

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-1-rbz ribozyme, which contains two stem regions (P1, P2), two internal loops, and a stem-loop (Figure 2D left).

The consistent result between LINE-1-rbz and LINE-1-ori 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 similar to each other. Both two ribozymes have two stems (P1, P2), two internal loops and a stem-loop region. Both stem-loop regions (SL2) are insensitive of ribozyme activity to mutations (Figures 1A, 1B). The internal loop regions (L1) of LINE-1-rbz and OR4K15-rbz are nearly identical except that C8 in OR4K15-rbz is replaced by U38 in LINE-1-rbz (Figure 3A). Given the conserved core regions of the LINE-1-rbz and OR4K15-rbz ribozyme, we named them lantern ribozymes because both are shaped like a Chinese lantern.

Figure 3

Consensus sequence of two lantern ribozymes. The consensus sequences (Figure 2E) for LINE-1-rbz and OR4K15-rbz were generated by R2R (27) 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 lantern ribozymes occurs between a conserved CA dinucleotide which locates inside the longer part of the internal loops (L1) linking P1 and P2. The internal loops are the regions with the highest conservation, which suggests their important roles in cleavage activity. Both stem regions (P1, 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 1A and 1B, providing additional support for the secondary structure models obtained here.

The secondary structure of OR4K15-rbz is a circular permutation of LINE-1-rbz self-cleaving ribozyme. Removing the non-functional loop regions (Figures 1A and 1B) allows us to recognize that the secondary structure of OR4K15-rbz is a circular permutated version of LINE-1-rbz. As shown in Figure 3A, after a permutation, OR4K15-rbz has two conserved internal loops with one loop identical to LINE-1-rbz and the other loop differed from LINE-1-rbz by one base. The only mismatch U38C in L1 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 lantern ribozymes were derived from our deep mutational scanning results. We obtained three hits with identical motifs all from the twister sister ribozyme family (Supplementary Table S3). Figure 3A further shows the comparison with two previously published twister sister structures (28, 29). The two internal loops of LINE-1-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) (28). Thus, it is possible to use the structure of twister sister ribozyme to build a homology model.

Homology modelling of the lantern ribozymes. To obtain the structure model of lantern ribozymes, we used template-based homology modelling with the twister sister ribozyme structure as the template. Figure 3B shows the homology modelled structure of two-stranded 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 L1 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 L1, with A11 at the cleavage site directed inwards. As shown in Figures 3C and Supplementary Figure S8, stem P1 of LINE-1-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-1-rbz is extended through trans non-canonical A12G36 and trans sugar edge-Hoogsteen A11C37. In addition to the triple base interaction in LINE-1-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 L1 in OR4K15-core is similar to L1 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-1-rbz).

Like the two twister sister ribozyme structures (28, 29), 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 (29) and A7 in three-way junctional twister sister ribozyme (28) 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 A11 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-1-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 lantern 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 L1s between lantern and the twister sister ribozymes. To address the above question, we further investigated the mutational effects on internal loop L1s. Most nucleotides in the L1 region of LINE-1-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 hydrogen-bonding interactions in the catalytic region. By comparison, mutations of C62 (C54 in three-way junctional twister sister ribozyme) at the cleavage site did not make a major change on the cleavage activity in previous studies (28, 29). This can be compared to the fact that single mutations of the corresponding nucleotide in LINE-1-rbz (C10) 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-1-core. As shown in Figure 4C, mutations on C10 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 (28, 29). More interestingly, A11U in LINE-1-core showed partial cleavage activity, whereas the corresponding mutation A63U in the four-way junctional twister sister ribozyme showed complete loss of cleavage activity (28, 29). These different mutation effects indicate that LINE-1-core may not adopt exactly the same structure and catalytic mechanism as twister sister ribozyme. In other words, lantern ribozymes may not be simply a minimal version of twister sister ribozymes.

We further examined the interactions around the cleavage site C10-A11 in LINE-1-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 C10-A11 in LINE-1-core (C40-A41 in OR4K15-core) are splayed away, with A11 (A41 in OR4K15-core) directed inwards and C10 (C40 in OR4K15-core) directed outwards. In addition, H-bonded interactions between C10 (C40 in OR4K15-core) and the two nucleotides 3′ of C10 could be found (Figure 4B, Supplementary Figure S9D). These structure differences may partially explain why mutations of C10 were detrimental to the catalysis in lantern ribozymes. Moreover, we did not observe interactions between C10 (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 lantern ribozymes. Thus, the stability of the current model structure based on the template from twister sister is uncertain.

Activity confirmation and biochemical analysis of LINE-1-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-1-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²⁺ or the enzyme strand was not included in the reaction, no cleavage could be observed even after 24h’s incubation (Supplementary Figure S10). This result further confirms that the RNA cleavage in lantern ribozymes is a catalytic reaction, not a spontaneous one. Moreover, the result confirmed that the stem loop SL2 regions in LINE-1-rbz and OR4K15-rbz (Figure 2) did not participate in the catalytic activity.

Figure 5

The core constructs of the two lantern ribozymes contain 35 and 31 nucleotides, respectively. They are the shortest and the second shortest 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, 12, 30–32). We confirmed that both LINE-1-core and OR4K15-core employed the same cleavage mechanism because the analog RNAs that lacked the 2′ oxygen atom in the -1 nucleotide (dC10 for LINE-1-core, dC40 for OR4K15-core) were unable to cleave (Figure 5B). We also investigated whether lantern 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 lantern 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 lantern ribozymes’ cleavage activity. At a concentration of 1mM, 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 lantern ribozymes have an equivalent or even higher cleavage ratio with Mn²⁺ than Mg²⁺. For LINE-1-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-1-core was ~0.04 min^-1 when measured in 10mM MgCl₂ and 100mM KCl at pH 7.5. Cleavage activities of the two lantern ribozymes were highly dependent on the concentration of divalent metal ion. The steep increase in rate constants of two lantern ribozymes plateaued at Mg²⁺ concentrations above 100mM, while for Mn²⁺ it plateaued at concentrations of 10mM. Consistent with the PAGE result, LINE-1-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-1-core when the concentrations of Mn²⁺ were less than 20mM, but higher, otherwise. Thus, there may exist important difference between OR4K15-core and LINE-1-core to explain this different bias toward ions. A detailed explanation for this difference may require high-resolution structure determination of the lantern ribozymes.

Discussion

The wide distribution of RNA self-cleaving activity in different species suggest that self-cleaving ribozymes play an important role in living organisms. Most self-cleaving ribozymes 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, 33). LINE-1 belongs to the family of retrotransposons, they make up about 17% of the human genome, with over 500,000 copies (34, 35). 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. 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 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 (36). 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 (24) 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 family 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 lantern ribozymes in the circular permutated form from one another (Fig. 3A). From the mutational information of LINE-1-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 family members in a circular permutated form. Currently, three classes of self-cleaving ribozymes, the hammerhead (7, 37–40), twister (12) and the hairpin (41) ribozymes are found in different species in multiple circular permutated forms. The discovery of permutated lantern ribozymes in the human genome for the first time suggests that we can find additional examples of circular permutated lantern ribozymes in other species.

More interestingly, we found that the internal loop regions of the lantern ribozymes share a high similarity with the known catalytic region of the twister sister ribozyme (2 mismatches only, Fig. 3A). The conserved cleavage sites suggest a similar structure and mechanism as shown in Fig. 3A. However, the maximum observed rate constant (k_obs) for twister sister ribozyme was ~ 5 min^{− 1} (32), compared to ~ 0.11 ± 0.01min^{− 1} for the lantern 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 (Fig. 3A) has the RA′ of 0.6. Two mismatches (A12U and U38A) in the catalytic region of the twister sister ribozyme (Fig. 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. Absence of this tertiary interaction in the lantern ribozyme may not be able to stabilize the structural form generated from homology modeling.

Indeed, we found some differences between the responses to mutations in the catalytic regions around the cleavage sites of lantern 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 (28, 29). However, unlike twister sister ribozymes, the cytosine at the cleavage site of lantern 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 lantern ribozymes show that the cytosine (C10 in LINE-1-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 (Fig. 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 lantern ribozymes and twister sister.

While the in vivo activities and functions of these two circular-permutated lantern ribozymes require further studies, their discovery may have important therapeutics and biotech implications. First, the lantern ribozymes have the simplest secondary structure with two stems and two internal loops, compared to all naturally occurring self-cleaving ribozymes. The previous minimum number of stems is three stems for hammerhead ribozyme (42). Second, this lantern 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 (43). A recent intracellular selection-based 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 (43). The lantern 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. 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-1_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-1-ori) and OR4K15 ribozyme (labelled as OR4K15-ori) followed the same protocol as our previous work (24). 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 m13f adp1 and template switching oligo TSO (Supplementary Table S4). The total cDNA was used to generate the RNA-seq library with primer P5R1 adp1 and P7R2 adp2 (Supplementary Table S4). The DNA-seq libraries were generated by amplifying the ribozyme mutant library with primer P5R1_m13f 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 numbers 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 (24) to analyze the RA of LINE-1-ori and OR4K15-ori. For Monte Carlo (MC) simulated annealing, we used the same weighting factor of 2 for both datasets as employed previously (24).

DMS and CODA analysis on the functional region of LINE-1 ribozyme

We identified the functional region of the LINE-1 ribozyme (labelled as LINE-1-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-1-rbz. Here, we used a doped mutation library instead of a library generated by error-prone PCR. The doped mutant library Rz_LINE-1_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 templates 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_m13f_ 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 5× ProtoScript II Buffer (New England Biolabs), 2 μl of 0.1 M DTT, 0.2 μl of Murine RNase Inhibitor (40 U/μl, 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_m13r and P5R1_m13f (Supplementary Table S4). While the DNA template used for transcription was used as the template for constructing the DNA-seq library with primer P5R1_m13f 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 (44) to generate the high-quality merged reads. We counted the read numbers 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_{RN Aseq} (var) N_{DN Aseq} (wt) / N_{DN Aseq} (var) N_{RN Aseq} (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-1-ori and LINE-1-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-1-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 (24) to analyze the RA′ of LINE-1-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 C10 and C40 were substituted by deoxycytosine 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, 1mM EDTA). The reaction products were then separated by denaturing (8 M urea) 15% PAGE and visualized by SYBR Gold (Thermo Fisher) staining. GeneRuler Ultra Low Range DNA Ladder (Thermo Scientific) was used in the PAGE analysis as the marker.

For the mutants of LINE-1-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, 1mM EDTA). The reaction products were then separated by denaturing (8 M urea) 15% PAGE.

Fluorescence-based kinetic analysis

The fluorescence-based kinetic analysis of the two lantern 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. In the ion-dependent assay, a final concentration of 20mM Tris-HCl solution (pH 7.5) and 100mM KCl was used. Different concentrations of MgCl₂ or MnCl₂ 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 H1. The first-order rate constant of ribozyme cleavage k_obs was calculated in a similar way as previously described (24). 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 lantern ribozymes

We matched the base-pairs of the lantern ribozymes to the most similar twister sister ribozyme (PDB: 5Y87) by map_align (45). 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 (46) developed by our group to optimize the structure model. The source code can be accessed at https://gitee.com/hongxu66/coda-map.

Acknowledgements

This work was funded by China Postdoctoral Science Foundation (fund no. 2021M702295 to Z.Z.) and the start-up fund from Shenzhen Bay Laboratory to Y.Z. This work was also supported in part by Australia 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.

Declarations

Data Availability

Illumina sequencing data for the LINE-1-ori, OR4K15-ori and LINE-1-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-1-core and OR4K15-core are available in ModelArchive (https://modelarchive.org/doi/10.5452/ma-chp4j for LINE-1-core; https://modelarchive.org/doi/10.5452/ma-bqqf8 for OR4K15-core).

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 Files

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.