Gel-free library preparation for next-generation RNA sequencing and small RNA quantification

Gel-free library preparation for next-generation RNA sequencing and small RNA quantification | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Gel-free library preparation for next-generation RNA sequencing and small RNA quantification Peter Dedon, Ruixi Chen, Lili Liu, Bo Cao This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7256873/v1 This work is licensed under a CC BY 4.0 License Status: Under Review Version 1 posted You are reading this latest preprint version Abstract Next-generation RNA sequencing (RNA-seq) is hampered by “primer dimer” (PD) artifacts and its quantitative performance reduced by polymerase fall-off (PF) at RNA modifications and secondary structures. Here we improved RNA-seq efficiency by incorporating (i) a post-reverse-transcription (RT) digestion of excess primers with Escherichia coli exonuclease I for PD mitigation, thus obviating gel purification during RNA-seq library preparation, and (ii) a high-processivity reverse transcriptase to increase full-length reads. A full factorial experimental design was applied to absolute quantification RNA sequencing (AQRNA-seq), the most accurate NGS-based method for quantifying small RNAs, using cDNA libraries constructed from E. coli small RNAs (> 85% tRNA) followed by sequencing, data processing, and data analysis. The novel PF and PD mitigation approaches increased AQRNA-seq sensitivity > 10-fold by minimizing PF and maximizing target RNA reads. By increasing sensitivity and obviating gel electrophoresis for removing PD, AQRNA-seq and other NGS-based RNA-seq methods can now be automated to increase throughput and reduce RNA sample size. Biological sciences/Biological techniques/Sequencing/RNA sequencing Biological sciences/Genetics/Epigenomics Biological sciences/Biochemistry/RNA RNA sequencing absolute quantification gel-free library preparation Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 9 Figure 10 Figure 11 Figure 12 Figure 13 Introduction Next-generation sequencing (NGS) is a fast, sensitive, and cost-effective technology for rapid sequencing of nucleic acids, with recent applications enabling mapping of the > 170 chemical modifications of the epitranscriptome 1 – 6 . Initially used for DNA sequencing tasks such as whole-genome sequencing, targeted sequencing, epigenetics, and metagenomics, NGS was subsequently adapted to sequencing all types of RNA (RNA-seq) 1 – 3 . Advancements in sequencing instrumentation, coupled with the development of specialized library preparation methods and data analytical pipelines, have broadened the scope of RNA-seq applications beyond measuring gene expression, such as identifying novel transcripts, mutations, and alternative splicing 7 , as well as quantitative mapping of post-transcriptional modifications 4 – 6 and single-cell resolution RNA-seq (scRNA-seq) 8 . Until recently, one notable limitation to most RNA-seq methods has been the inability to accurately measure the absolute abundance of RNA molecules within a cellular pool. Biased ligation of sequencing adaptors 9 – 11 and premature termination of reverse transcription (RT), due to polymerase fall-offs (PF) caused by post-transcriptional RNA modifications or secondary structures 4 , 12 , are among key challenges to quantitative accuracy. These challenges, however, were addressed in the development of A bsolute Q uantification RNA-seq (AQRNA-seq) 13 , 14 , a specialized RNA-seq method for quantitative profiling of all small RNA species (< 200 nt) in any cell, tissue, and organism. Quantitative accuracy is achieved by using custom DNA linkers to minimize ligation biases, a two-step ligation approach to fully capture truncated complementary DNAs (cDNAs), and an optional AlkB demethylation 15 step to reduce PF by partial enzymatic removal of several types RNA methylation, thereby increasing the proportion of full-length cDNAs. As a result, AQRNA-seq achieves linearity between sequence read count and biological copy number, as validated in studies involving a reference library of 963 miRNAs, pooled RNA oligonucleotide standards with variable lengths, and an orthogonal tRNA abundance dataset derived from quantitative Northern blotting 16 . Despite the superior quantitative performance of AQRNA-seq compared to other methods, it suffered like most other NGS-based RNA-seq approaches from limitations posed by library preparation workflow and sensitivity. For example, all NGS-based RNA-seq methods suffer from highly abundant “primer dimers” (PD), which are formed by ligation of the 5' DNA linker to the 3' end of RT primers during library preparation. Due to their small sizes, PD can form clusters on the flow cell more efficiently than target RNA molecules, thus consuming the sequencing capacity and lowering sensitivity by producing unusable reads 17 , 18 . The traditional approach to reducing the impact of PD is to resolve the cDNA products by electrophoresis followed by gel extraction (GE) to size-select target fragments 19 . However, GE suffers from highly variable accuracy of manual gel excision and inefficient extraction of target fragments 20 , which collectively reduce sensitivity. More importantly, the time-consuming manual labor involved in GE renders library preparation for AQRNA-seq and other RNA-seq methods unsuitable for automation and high-throughput applications. In addition to PD concerns, the quantitative accuracy of AQRNA-seq can be limited by high levels of truncated cDNAs persisting for many tRNAs 14 , especially those containing modifications insensitive to AlkB treatment. While being fully captured by the two-step ligation approach, these truncated cDNAs may fail to unambiguously align with the reference sequence library. We have now solved these problems by combining (1) a novel approach to PD mitigation involving digestion of excess RT primers using E. coli exonuclease I 21 immediately after RT and (2) a highly processive reverse transcriptase (RTase), SuperScript IV, that reduces PF to increase full length reads. A full factorial experimental design (Fig. 1 ) was conducted to contrast these novel approaches with their original counterparts. Our results suggest that exonuclease I is significantly more efficient in reducing PD compared to GE, increasing RNA-seq sensitivity while in the meantime obviating human involvement in gel extraction. Similarly, SuperScript IV proved superior to AlkB in mitigating PF, especially in the presence of diverse RNA modifications, enhancing the quantitative accuracy of AQRNA-seq. The use of both exonuclease I and SuperScript IV can be readily integrated into any RNA-seq workflow, allowing full automation of RNA-seq library preparation. RESULTS EXOI enhances sensitivity and enables automation of RNA-seq workflows. The presence of PD in sequencing libraries, along with the need for gel extraction for removing PD, substantially deteriorates the performance of sequencing methods, particularly those targeting small RNAs. While gel extraction is commonly used for PD mitigation during library preparation, it imposes constraints on method sensitivity and automation. Here, we introduce EXOI as a simple, cost-effective, and universally applicable approach that substantially outperforms GE in efficiency and holds potential for enhancing sensitivity and facilitating automation, as well as high-throughput implementation, across all RNA-seq methods. In AQRNA-seq, ligation of 5’ DNA linkers to the 3’ end of RT primers leads to the formation of PD 13 . Hence, the ExoI nuclease is added immediately following RT to digest excess RT primers, thereby inhibiting the formation of PD. Electrophoretic analysis (3% agarose gel) of the cDNA libraries subsequent to PCR amplification revealed a lower density of fragments below the 200-bp size marker for samples subjected to protocols V3 and V4 (incorporating EXOI), as compared to protocols V1 and V2 (incorporating GE), respectively, confirming the efficacy of EXOI in inhibiting the formation of PD (Fig. 2 a). Furthermore, the proportion of quality-filtered reads with an insert length of > 10 bp was significantly higher ( p < 0.05) for samples subjected to EXOI as compared to GE (Fig. 2 b, Supplementary Table 1 ), suggesting a superior performance of EXOI in reducing PD and short inserts. More importantly, our results showed that the incorporation of EXOI substantially increased the concentration of cDNA libraries following qPCR amplification, exceeding that achieved with GE by over 11-fold (Fig. 2 c). This highlights the potential of EXOI to significantly enhance the sensitivity of AQRNA-seq. Finally, substituting GE with EXOI eliminates the need for human intervention, thus enabling the future automation and high-throughput implementation of AQRNA-seq. SS demonstrates potential for enhancing quantitative accuracy through its superior ability in attaining full-length reads. Having addressed challenges related to sensitivity and automation, we proceeded to optimize the quantitative accuracy of AQRNA-seq, which is contingent upon precise alignments of sequence reads against the reference library. In AQRNA-seq, the synthesis of truncated cDNAs due to premature RT terminations inevitably compromises the analysis and could be detrimental for differentiating between highly similar sequences. While ALKB can be used reduce PF, one cause of premature RT terminations, thereby increasing proportions of full-length reads, its efficacy is limited to a few methyl modifications, leaving most other modifications unaffected. Therefore, we assessed SS as an alternative approach for PF mitigation, which employs SuperScript IV, a high-processivity RTase. We identified an ascending trend in the proportion of full-length reads for individual tRNAs across protocols, with the mean proportion progressively increasing from 0.11 for V1, to 0.14 for V3, then to 0.35 for V2, and finally to 0.43 for V4 (Fig. 2 d, Supplementary Table 2 ). This suggests the potential of both SS and EXOI in increasing proportions of full-length reads, with SS showing the predominant influence. A linear regression model was constructed to assess the effect of PF and PD mitigation approaches, tRNA identity, and their potential interactions on full-length read proportions. While all main factors and interactions included in the model demonstrated significance ( p < 0.05) based on F tests ( Supplementary Table 2 ), interpretation was only deemed appropriate for the interaction among the three main factors. Specifically, post hoc comparisons were performed to compare the model-reported estimated marginal means (EMM) of the proportion between ALKB and SS samples conditional to each PD mitigation approach and each tRNA ( Supplementary Table 2 ). Our results revealed a significantly ( p < 0.05) higher proportion of full-length reads resulted from SS compared to ALKB for 39 tRNAs, regardless of the PD mitigation approach used (Fig. 2 d). Importantly, for 7 of these tRNAs, the difference between EXOI and GE was significant conditional to the incorporation of SS, but not ALKB, showcasing the interactive effect between PF and PD mitigation approaches. Therefore, these results confirm the superior performance of SS over ALKB in attaining full-length reads, particularly in conjunction with EXOI, suggesting its potential for enhancing the quantitative accuracy of AQRNA-seq. Both SS and EXOI resulted in variations in quantitative results across different protocols. Our findings suggest that SS and EXOI enhance the quantitative accuracy and sensitivity of AQRNA-seq. We next sought to determine whether these enhancements would manifest as variations in quantitative results. Pairwise Pearson’s correlation coefficients (r) based on normalized tRNA abundances revealed that the incorporation of SS resulted in a lower correlation (r = 0.84) with the original AQRNA-seq protocol compared to the incorporation of EXOI (r = 0.96) (Fig. 2 e, Supplemental Table 3 ). The lowest pairwise correlation was observed between the original protocol and the most revised version, which incorporates both SS and EXOI. To objectively rank the quantitative performance of the protocols, we assessed the correlation between the percentages of tRNAs within the tRNA pool obtained from each protocol and those derived from a quantitative Northern blotting protocol 16 , considering higher correlations with this orthogonal dataset were indicative of more accurate quantification. Our results showed that both SS and EXOI contributed to enhanced correlations with the orthogonal dataset, with the highest correlation observed for protocol V4 (Fig. 2 f, Supplemental Table 4 ). This suggests the potential of both SS and EXOI to enhance the quantitative performance of the original protocol, establishing V4 as the optimal protocol for the library preparation of AQRNA-seq. The estimated abundance for individual tRNAs is more dependent on PF mitigation approach than PD mitigation approach. Having demonstrated that both SS and EXOI collectively led to an enhancement in quantitative performance, we next proceeded to assess and compare their respective impacts on the quantification of individual tRNAs. While our results revealed consistent quantification across protocols for most tRNAs, noticeable variations in estimated abundance were observed for selected tRNAs, such as Ile-GAT-1 and Phe-GAA-1 (Fig. 3 a, Supplementary Table 5 ). Hierarchical clustering based on rlog-transformed and normalized tRNA abundance formed two distinct clusters of samples by PF mitigation approach (ALKB vs SS), while not by PD mitigation approach (GE vs EXOI). Consistent with this, projecting the samples onto the principal component analysis (PCA) score plot unveiled a more pronounced segregation of samples by PF as compared to PD mitigation approach (Fig. 3 b), further indicating a greater impact of SS relative to EXOI on estimated abundances. Importantly, the PCA score plot suggested an interaction effect between PF and PD mitigation approaches, as the separation distance between GE and EXOI samples subjected to SS (comparing V2 and V4 samples) was considerably larger than those subjected to ALKB (comparing V1 and V3 samples). To gain mechanistic understanding into the variation in estimated abundance attributable to SS and EXOI, we identified individual tRNAs exhibiting differential estimated abundance (|log 2 fold change| > 1 along with p < 0.05) across protocols as affected by the incorporation of SS and/or EXOI (Fig. 3 c, Supplementary Table 5 ). While the estimated abundance remained consistent across all protocols for most tRNAs, corresponding well with the high correlations observed among the protocols, we identified 5 tRNAs (Arg-CCG-1, Ile-GAT-1, Ile2-CAT-1, Leu-TAA-1, and Val-GAC-1) that showed significantly higher estimated abundances in SS compared to ALKB samples (Fig. 3 d), and 3 tRNAs (fMet-CAT-1, Leu-CAG-2, and Pro-CGG-1) showing lower estimated abundances (Fig. 3 e). In addition, although none of the tRNAs were differentially estimated between GE and EXOI samples, the interaction effect was deemed significant for 6 tRNAs (Met-CAT-1, Ile2-CAT-2, Phe-GAA-1, Phe-GAA-2, Thr-CGT-1, and Thr-CGT-2), suggesting that differences in estimated abundance between ALKB and SS samples were further contingent upon PD mitigation approaches (Fig. 3 f). SS showed the potential to correct underestimated tRNA abundances by ALKB. The most notable difference in estimated abundance between ALKB and SS samples was observed for Ile-GAT-1 (SS/ALKB = 13.6), followed by Ile2-CAT-1 (SS/ALKB = 5) and Val-GAC-1 (SS/ALKB = 5) (Fig. 3 d, Supplementary Table 5 ). To elucidate mechanisms driving these differences, we modeled the RT process of the corresponding tRNAs by cross-referencing the position-specific percent coverage and identity data generated in this study ( Supplementary Table 6 ) with the positional information of tRNA modifications available on Modomics 22 . This is illustrated in Fig. 4 for Ile-GAT-1, for which the normalized read counts were 1,803 and 35,884 following protocols V1 and V4, respectively, suggesting a 20-fold difference ( Supplementary Table 5 ). Regardless of the protocols used, our results suggested the presence of a highly abundant 5’ tRNA fragment, as most read alignments (3’ to 5’ direction) initiated at position 59, rather than the 3’ end (Fig. 4 a,b). The RT of this fragment, however, was thwarted by the presence of 3-(3-amino-3-carboxypropyl)uridine (acp 3 U) at position 45B, manifested as the sharp downward spikes (i.e., reduced percent coverage and identity) shown in different panels of Fig. 4 . When protocol V1 was used, there was a pronounced reduction in percent coverage at position 45B, followed by considerably reduced percent coverage from position 45A onwards, a hallmark of PF at acp 3 U causing premature RT terminations (Fig. 4 a). Consequently, the cDNAs synthesized for the tRNA fragment were likely too short (i.e., ~ 12 nt from position 45C to 59) to produce sequence reads that could be aligned with the reference sequence, leading to an underestimation in abundance. Conversely, following protocol V4, while dropping drastically at position 45B, the percent coverage returned to nearly 100% for positions past 45B, indicating that the high-processivity RTase copied past acp3U, albeit causing a -1 mutation, to create longer cDNAs that could be aligned with the reference sequence (Fig. 4 b). This critical difference in how tRNA modifications may adversely impact the RT process remained consistent for Ile2-CAT-1 and Val-GAC-1, wherein acp 3 U triggered considerable premature RT termination with V1 but to a far lesser extent with V4 ( Extended Data Figs. 1 a,b and 2 a,b). Therefore, our findings illustrate the capability for SS to correct underestimation of tRNA abundances observed in samples processed with the ALKB method. A notable feature of our analysis of Ile-GAT-1 was a marked reduction in percent identity around position 69 for Ile-GAT-1 (Figs. 4 ), suggesting the presence of a modification that has not been annotated in the Modomics database or a highly stable RNA secondary structure. This Putative modification/2° structure, potentially in combination with the pseudouridine (Y) at position 65, was correlated with near-stoichiometric PF in SS samples, which was partially mitigated by the ALKB method, as demonstrated by the higher percent coverages observed from positions 76 to 60 in ALKB samples (Fig. 4 a) compared to SS (Fig. 4 b). Thus, our findings not only suggest the potential for AQRNA-seq to unearth previously unidentified RNA modifications, along with their associated mutation signatures, but also highlight situations where ALKB may outperform SS in handling specific modifications. The impact of tRNA modifications on the accuracy of tRNA quantification varied by protocol, location, and chemical structure. Since uninterrupted RT was demonstrated to be important for accurately quantifying individual tRNAs, we next characterized the potential of all modifications present on E. coli tRNAs to interfere with RT for each protocol, aiming to provide a systematic overview of the effect of each modification on quantification of the tRNA. We identified 42 tRNA isoacceptors in the RF00005 alignments ( https://rfam.org/family/RF00005 ) that were successfully matched to reference sequences in our library and contain 28 unique modifications located at 24 positions in the alignments (Fig. 5 a, Supplementary Table 7 ). While most modifications tend to cluster in proximity within the same structural domain, the modification Y, as one of the most frequent modifications, spans across multiple domains. The impact of tRNA modifications on quantification was inferred from the position-specific percent coverage data, with lower percent coverages inversely correlated with more adverse impacts. We found that protocols V1-V3 resulted in insufficient coverage (i.e., ≤ 100 reads) at the 5’-end for 8, 7, and 7 of 42 tRNAs, respectively, consistent with premature RT terminations and/or low sensitivity resulting in reduced read counts for these tRNAs ( Extended Data Fig. 3 , Supplementary Table 6 ). Conversely, protocol V4 maintained sufficient coverage (i.e., > 100 reads) toward the 5’-end of all tRNAs, suggesting limited premature RT termination due to SS for PF mitigation and increased sensitivity conferred by EXOI for PD removal. To systematically assess and compare the impact of modifications on quantitation across the four protocols, we analyzed the percent coverage data at each position with modified residues using linear regression, focusing only on tRNAs carrying modifications at the specific position. Due to insufficient coverage at selected positions across different protocols, the analysis was conducted for 19 positions with a minimum of three data points for each protocol. The percent coverage differed significantly across protocols at 18 of these positions (Fig. 5 b, Supplementary Table 8 ). Post hoc pairwise comparisons revealed a significantly higher percent coverage for V2 and V4 compared to V1 and V3 for all 18 positions, confirming the superior capability of SS over ALKB to prevent premature PF at modifications. Consequently, modifications imposed less of an adverse impact on tRNA quantification with protocols incorporating SS compared to protocols incorporating ALKB. Notably, a significantly higher percent coverage was observed for V4 compared to V2 at 7 positions and for V3 compared to V1 at 6 positions. Again, this suggests an important role for EXOI in further enhancing quantitative accuracy of the protocols, likely due to its ability to enhance sensitivity, especially for tRNAs with low abundance. For each protocol, substantial variations in percent coverage were evident at certain positions across tRNAs with modifications, suggesting that the presence of different modifications at a specific position may have differential impacts on quantification. Therefore, we further scrutinized the percent coverage in the presence of individual modifications at all possible positions. In general, the percent coverage, normalized to the maximum value at the corresponding position in the presence of modifications, may exhibit substantial variability both across modifications at a specific position and across positions where a specific modification may localize (Fig. 5 c). For instance, among the 20 tRNAs carrying 11 distinct modifications at the wobble position 34, the normalized percent coverage averaged for each modification ranged from 0.8% to 73% for V1, 5% to 93% for V2, 1% to 80% for V3, and 2% to 96% for V4. The normalized percent coverage was subsequently classified into four categories: ultra-high (75–100%), high (50–75%), low (25–50%), and ultra-low (0–25%), wherein instances falling with the ultra-low category were further investigated, as they likely presented the greatest challenges in the quantification of AQRNA-seq ( Supplementary Table 9 ). Across different protocols, V4 resulted in the highest proportion of modification-position combinations (26 out of 44) exhibiting ultra-high normalized percent coverages while the lowest proportion of combinations (7 out of 44) showing ultra-low or insufficient data. In particular, our results revealed modification-position combinations manifesting high/ultra-high normalized percent coverages with protocol V4 but ultra-low normalized percent coverages with at least one of the other protocols, including N4-acetylcytidine (ac4C) at position 34, inosine (I) at position 34, 1-methylguanosine (m1G) at position 37, and 2-thiocytidine (s2C) at position 32. Conversely, all modification-position combinations manifesting ultra-low normalized percent coverages with V4 also fell into the ultra-low category with all the other protocols. This reaffirms the superior performance of protocol V4 as compared to the other protocols in minimizing the adverse impacts of modifications on quantification, corresponding well with our findings at the positional level. Consequently, further inferences in this regard were made exclusively based on results from protocol V4. While we identified two modification-position combinations (5-methyluridine [m5U] at position 53 and Y at position 54) showing insufficient data, both instances were associated with a minor tRNA SeC-TCA-1, thereby possibly attributable to the low abundance of this tRNA. Furthermore, we identified five modification-position combinations exhibiting ultra-low normalized percent coverages, including 2'-O-methylcytidine (Cm) at position 34, 5-carboxymethylaminomethyl-2'-O-methyluridine (cmnm5Um) at position 34, 2-lysidine (k2C) at position 34, 2-methylthio-N6-isopentenyladenosine (ms2i6A) at position 37, and Y at position 32. However, upon further examination of the read alignments of the corresponding tRNAs, it became evident that the observed ultra-low normalized percent coverages of Cm, cmnm5Um, and Y could be attributed to the presence of ms2i6A at position 37 of the corresponding tRNAs, which precipitated a sharp decline in percent coverage, followed by diminished percent coverages thereafter ( Extended Data Fig. 4 ). Hence, our findings suggested a high likelihood for ms2i6A and k2C, though not necessarily the other modifications, to induce premature PF and thereby potentially impose detrimental effect on tRNA quantification. Quantitative mapping of tRNA modifications varies among the protocols. One consequence of minimizing PF to increase the accuracy and sensitivity of small RNA quantification was a reduced ability to quantitatively map tRNA modifications through position-specific read pileups. However, quantitative modification mapping of tRNA modifications is still possible for those that cause the RT to insert the wrong nucleotide to create a mutation. We tested this approach to modification mapping by quantifying the position-specific percent sequence identity (i.e., 1-mutation frequency) at all possible positions with known modifications in E. coli (Fig. 5 d). The percent identities fell below 90% for 2, 5, 2, and 6 modifications subjected to protocols V1, V2, V3, and V4, respectively, suggesting the potential for quantitative mapping of these tRNA modifications. For example, the percent identity ranged from 0.2% to 1.7% across all protocols for wobble I, consistent with its well-documented ability to induce a T-to-C mutation during RT 23 . In addition, the percent identity was notably lower with protocols V2 and V4 as compared to protocols V1 and V3 for acp3U, ms2i6A, m1G, and cmnm5Um. These results point to the utility of EXOI/SS-modified AQRNA-seq quantitative mapping of several modifications, with the potential for identifying predictive mutation signatures by base profile surrounding target positions. DISCUSSION Here we optimized the highly accurate AQRNA-seq 13 , 14 method to improve quantitative accuracy and sensitivity, while obviating the gel purification step to boost efficiency and enable automation for greater throughput and scalability. These advancements were achieved with two enzymatic changes: (i) adding a step for post-RT digestion of excess primers using E. coli exonuclease I for effective PD mitigation, and (ii) use of SuperScript IV—a high-processivity RTase—to circumvent RNA modifications during cDNA synthesis for improved PF mitigation. These changes substantially outperformed their original, labor-intensive counterparts of gel purification for PD mitigation and use of AlkB demethylase to remove a subset of RNA methylations for PF mitigation. Both EXOI and SS can be applied to almost any NGS RNA-seq method to facilitate automation. GE is a routinely used approach in NGS RNA-seq library preparation to reduce the impact of PD in sequencing libraries and is implemented in various workflows such as CUT&RUN, ChIP-Seq, MeDIP-Seq, and the original AQRNA-seq 14 , 24 – 27 . However, its efficiency varies depending on sample characteristics, experimental conditions, gel cutting techniques, and operator experience, frequently leading to reduced sensitivity due to the loss of target fragments 28 . Furthermore, the labor-intensive nature of GE hinders the automation of NGS RNA-seq methods, thus constraining their scalability and throughput. To tackle this issue, we tested EXOI—post-RT digestion of excess RT primers using ExoI nuclease (Fig. 1 )—as a more efficient and cost-effective alternative for mitigating PD. Notably, compared to GE, EXOI substantially enhances the sensitivity of AQRNA-seq through a reduced loss of target fragments and a more effective mitigation of PD (Fig. 2 b,c). Furthermore, EXOI eliminates the need for human involvement, thus facilitating automation and enabling higher-throughput small RNA quantification. While alternative approaches for PD mitigation are available, such as bead-based size selection ( https://www.cytivalifesciences.com/en/us/news-center/better-data-via-size-selection-10001?ssp=1&setlang=en&cc=US&safesearch=moderate ) and chemically modified adaptors 29 , they have significant drawbacks that limit their use. For instance, bead-based size selection lacks sufficient resolution for small size differences between PD and target fragments (e.g., 20 nt for miRNAs) and may not be applicable to low-input samples 30 . The incorporation of chemically modified adaptors, on the other hand, requires the case-by-base redesign of DNA linkers. In contrast, the simplicity of EXOI makes it compatible and easily adaptable to most library preparation workflows. Moving forward, the incorporation of EXOI in AQRNA-seq and related methods will likely bolster the speed and accuracy of quantitative investigations into all forms of biological targets. Post-transcriptional modifications are a significant cause of PF during RT, which reduces quantitative accuracy due to imprecise alignments of reads with the reference sequence library. While the ALKB method can be used to mitigate PF by enzymatically removing several simple methyl modifications such as m 6 A, m 1 A, and m 3 C 30 – 32 , most other RNA modifications not affected by this approach and demethylation efficiency is ribonucleoside-specific and highly variable due to the oxygen sensitivity of the Fe-S cluster of AlkB 14 . PF was an important contributor to the numerous truncated cDNAs synthesized for many tRNAs in the original AQRNA-seq protocol with ALKB 14 . To further enhance the quantitative accuracy of AQRNA-seq, we tested SS as an alternative approach for PF mitigation, employing a high-processivity RTase to circumvent diverse RNA modifications during cDNA synthesis 33 . Our results show that SS significantly increases the proportion of full-length reads for nearly 80% of tRNAs, contributing to a substantial enhancement in quantitative accuracy (Fig. 3 ). Cross-referencing our sequencing results with the positional information of tRNA modifications available on Modomics 22 provided insights into the superior performance of SS in averting PF due to a broad spectrum of modifications at various locations. For instance, while acp 3 U at position 45B tended to induce PF in ALKB samples, it was bypassed during RT in SS samples, resulting in either a -1 mutation or a point mutation due to mismatched nucleotide incorporation, both of which are useful for modification mapping purposes. On a broader scale, SS yielded a significantly higher percent coverage compared to ALKB in the presence of modifications at 18 out of the 19 positions examined (Fig. 5 ), corresponding well with the significant difference in PF mitigation efficacy between the two approaches. Since the ability of modifications to trigger premature RT termination varies depending on their characteristics and location, as well as the RTase used, we aimed to provide a high-resolution overview of the impact of known tRNA modifications in E. coli on the quantitative analysis with AQRNA-seq. We found that the impact varies considerably across modifications and even across different positions carrying the same modification. This suggests intricate interactions between modifications and RT, which warrant further investigation. Importantly, despite the use of the EXOI and SS methods, we identified modifications exhibiting a high propensity for PF and thus reducing accurate quantification of the corresponding tRNAs, including ms 2 i 6 A at position 37 and k 2 C at the wobble position. In E. coli , while k 2 C34 is found specifically in tRNA-Ile2-CAT-1, where it converts both the codon recognition (AUG to AUA) and amino acid specificity (from methionine to isoleucine) of the tRNA 34 , ms 2 i 6 A37 is commonly found in tRNAs decoding UNN codons where it prevents codon misreading by stabilizing codon-anticodon interactions 35 . As such, compared to k 2 C34, ms 2 i 6 A37 may represent a greater challenge to quantification due to a larger number of tRNAs affected. Consistent with our findings, substantial PF due to ms 2 i 6 A has been highlighted in other studies 36 and leveraged to achieve quantitative mapping of the modification at single-nucleotide resolution 37 . The original AQRNA-seq method and especially the present revision incorporating EXOI and SS are aimed at sensitive and accurate quantification of small RNAs and not intended for mapping the small RNA epitranscriptome. Indeed, NGS-based modification mapping by RT-induced falloff or mutation is a well-established technology 15 , 36 – 40 . In summary, we developed a highly sensitive and readily-automated NGS RNA-seq method for small RNA quantification based on the incorporation of EXOI to reduce PD and SS to increase full-length reads. By increasing sensitivity and eliminating manual intervention throughout the library preparation process, AQRNA-seq and other NGS-based RNA-seq methods can now be fully automated to increase throughput and reduce the required amount of RNA input, thus facilitating the study of the human RNome 41 . METHODS The AQRNA-seq data analytical pipeline v 2.0. The original data analytical pipeline 14 was substantially reengineered for improved accuracy, efficiency, and user experience, resulting in the development of AQRNA-seq data analytical pipeline v2.0 ( Extended Data Fig. 5 ). Simplified procedures involved in the pipeline are described in this section, and the complete pipeline is available at GitHub (LINK to be added) with a user manual including extensive annotations for individual steps. In the revised protocol, fastp v 0.23.4 42 is used to assess the quality of raw sequence reads, and Cutadapt v 4.2 43 is used to clip both the upstream and downstream adaptors from the quality-filtered reads. For each sample, the forward and reverse reads are assembled using PEAR v 0.9.10 44 based on the identification of an overlap with statistical significance, producing a unified FASTQ file containing consensus reads in the forward read direction, which is then converted to the FASTA format using a custom Bash script. Subsequently, unique reads along with their occurrences within each sample are extracted using fastxtoolkit v0.0.13 ( http://hannonlab.cshl.edu/fastx_toolkit/ ), followed by combining the results across all samples to generate a FASTA file containing all unique reads present in at least one of the samples. The FASTA file is then processed and used as the input for the BLAST (Basic Local Alignment Search Tool) analysis to align the unique reads against a custom reference sequence library encompassing members of the target RNA species (e.g., tRNA) using blastn v 2.6.0 45 . Custom R scripts are used to filter raw blastn results based on use-defined thresholds and to resolve the mapping of a single query sequence to multiple reference sequences based on e-value, bit-score, percent identity, and query coverage. Finally, the processed blastn results are cross-referenced with the occurrences of unique reads in samples to generate a raw tRNA abundance matrix, using a custom R script, wherein the read counts are shown for each subject in the reference sequence library within each sample. Detailed instructions for constructing reference sequence libraries can be found in the section “Construction of the E. coli BW25113 tRNA reference sequence library”. Bacterial strains, culturing conditions, and small RNA isolation. E. coli strain BW25113 was selected for use in this study as it has been extensively studied as a model strain and represents the parent strain of the Keio Collection 46 , thereby facilitating the comparison of results derived from this study and potential subsequent studies. The stock culture was preserved in lysogeny broth (LB; Fisher BioReagent) with 10% DMSO at -80°C. Frozen stock culture was streaked on a LB agar plate, followed by incubation at 37°C for 16–18 h. Single colonies from the freshly streaked plate were inoculated into 5 mL LB broth, followed by overnight incubation on a spinning wheel (New Brunswick Scientific) at 37°C. Bacterial cultures were then sub-cultured (1:100) in 25 mL LB broth, followed by incubation in a shaking incubator (Shon’s Scientific Refrigeration) at 37°C with 250 rpm. Upon entry into the late-exponential phase, as indicated by optical density (OD) 600 values between 0.6–0.7, the bacterial culture was aliquoted (1 mL) into sterile Eppendorf tubes. Bacterial cell pellets were harvested by centrifugation at 500 × g and 4°C for 5 min and subsequently washed by resuspension in 1 mL of 1X phosphate buffered saline (PBS; pH = 7.4; Thermo Fisher Scientific), followed by centrifugation using the same settings. Isolation of small RNAs from the bacterial cells was conducted following the user guide of the PureLink miRNA Isolation Kit (Thermo Fisher Scientific) with modifications. Briefly, bacterial cells were lysed by adding 1 mL of the TRIzol Reagent (Sigma-Aldrich) per 1 × 10 7 cells (equivalent to approximately 1 mL of bacterial culture in late exponential phase), followed by incubation at room temperature for 5 min. To isolate total RNAs, the cell lysate was incubated with 0.2 mL of chloroform (Macron Fine Chemicals) at room temperature for 3 min, followed by centrifugation at 12,000 × g and 4°C for 15 min. The upper aqueous phase containing total RNAs was then transferred into a sterile RNase-free tube. Subsequently, small RNAs were separated from large RNAs by passing the total RNAs through two Spin Cartridges to selectively bind large and small RNAs to the membrane, with 35% and 70% ethanol, respectively. Following two wash steps, small RNAs were eluted in RNase-free water by centrifugation at 16,000 × g for 1 min. Library preparation. cDNA libraries were prepared using 75 ng of E. coli small RNAs as the input, following four variations of the AQRNA-seq library preparation protocol, each with three biological replicates, as illustrated in Fig. 1 . Protocol variant 1 (V1) represents the original protocol initially published in Hu et al. 14 and later elaborated on in Chen et al. 13 , while protocol variants 2 (V2), 3 (V3), and 4 (V4) involve alternative adaptations of the original protocol with respect to the approaches for mitigating polymerase fall-off (PF) and/or primer dimers (PD). In the following sections, protocol V1 will be briefly recapitulated, highlighting the molecular processes involved in individual steps, and the alternative protocols (V2, V3, and V4) will be introduced, emphasizing the specific adaptations made in comparison to the original protocol. Full protocols with step-by-step instructions are provided in Supplementary Notes 1–4 to ensure reproducibility and comprehensive evaluation of the study. Detailed information regarding the DNA and RNA oligonucleotides used in the protocols can be found in Supplementary Table 10 . Detailed information regarding the reagents, chemicals, and kits/columns can be found in Supplementary Table 11 . Protocol V1. Small RNAs mixed with an internal standard were dephosphorylated using Shrimp Alkaline Phosphatase (New England Biolabs) through incubation at 37°C for 30 min and then at 65°C for 5 min to both deactivate the enzyme and denature the RNAs. T4 RNA Ligase 1 (New England Biolabs) was used to ligate Linker 1 to the 3’ end of RNAs through incubation at 25ºC for 2 h and then at 16°C for 16 h. To mitigate PF, the AlkB demethylase (ArrayStar) was used to remove post-transcriptional methylations potentially present on RNAs through incubation at 25°C for 2 h (i.e., the AlkB approach). Excess Linker 1 not ligated to RNAs was then deadenylated using 5’-deadenylase (New England Biolabs) through incubation at 30°C for 1 h and subsequently digested with RecJf exonuclease (New England Biolabs) through two 30-min incubations at 37°C. Reverse transcription (RT) was initiated with a 2-min incubation at 80°C to anneal the RT primers to the RNA template, followed by a 2-h incubation at 50°C to synthesize the cDNA strand using the PrimeScript™ reverse transcriptase (RTase; TaKaRa) and a 15-min incubation at 70°C to deactivate the enzyme. After the synthesis of the cDNA strand, the RNA template was hydrolyzed using sodium hydroxide (Sigma-Aldrich) through incubation at 95°C for 3 min, followed by immediate neutralization using hydrochloric acid (VWR). T4 DNA ligase (New England Biolabs) was used to ligate Linker 2 to the 3’ end of cDNAs, and the excess Linker 2 were subsequently removed following the same procedures as described for Linker 1. Sequencing adaptors were then attached to both ends of target cDNA fragments through polymerase chain reaction (PCR), and the PCR products were subsequently subject to 3% agarose gel electrophoresis to separate target fragments from unwanted PD. To mitigate the carry-over of PD into sequencing libraries, gel blocks containing target fragments were excised, followed by gel extraction to retrieve and purify the target fragments (i.e., the GE approach,. Figure 1 ). Protocol V2. This variation of the protocol incorporates a novel approach (SS) to mitigate PF during the RT reaction, a predominant cause of premature RT termination (Fig. 1 ). Specifically, the SuperScript IV RTase (Thermo Fisher Scientific) was used in substitution of the PrimeScript RTase (TaKaRa) for synthesizing the cDNA strand. SuperScript IV has demonstrated superior processivity (i.e., the ability to consecutively add nucleotides without releasing the RNA strand) as compared to other RTases 33 , indicating its potential to minimize PF induced by a broader range of post-transcriptional modifications, not just methylations, thereby facilitating the synthesis of full-length cDNAs. As a result, the AlkB demethylation step was eliminated from the protocol. Protocol V3. This variation of the protocol incorporates a novel approach (EXOI) to mitigate the carry-over of PD into the constructed cDNA libraries (Fig. 1 ). Specifically, immediately following the RT reaction, exonuclease I from E. coli (New England Biolabs) was used to digest excess RT primers in the 3’ to 5’ direction, thereby inhibiting the formation of PD (i.e., ligation of DNA linkers to the 3’ end of RT primers) from the first place. As a result, the gel extraction step was eliminated from the protocol. Protocol V4. This variation of the protocol incorporates both the novel approach to mitigate PF (SS), as previously described for protocol V2, and the novel approach to mitigate PD (EXOI), as previously described for protocol V3 (Fig. 1 ). Consequently, both the AlkB demethylation and the gel extraction steps were obviated by using this protocol. Library quality assessment and sequencing. The constructed cDNA libraries were submitted to the MIT BioMicro Center, along with custom sequencing primers ( Supplementary Table 10 ). The quality of the libraries was assessed using AATI Fragment Analysis (Agilent) and LightCycler 480 Real-Time PCR System (Roche), and 75-bp Paired-End sequencing was performed on an Illumina MiSeq platform with the v3 reagent kit (Illumina). Data analysis. C onstruction of the E. coli BW25113 tRNA reference sequence library . To maximize the accuracy for quantifying the landscape of tRNAs, we followed a custom workflow developed in our laboratory to construct the reference sequence library for use in the data analytical pipeline. Specifically, the mature tRNA sequences of E. coli BW25113 were obtained from GtRNAdb 47 , 48 , followed by the removal of predicted pseudogenes and tRNAs with unknown or undetermined isotypes. “CCA” was manually appended to sequences lacking a 3’ overhang ending with “CCA”, based on the structural alignments generated using domain-specific covariance models in COVE 48 . Subsequently, all “U” in the sequences were replaced with “T”, and identical sequences representing multiple gene copies encoding the same isodecoder were deduplicated, preserving only one sequence in the library as a representative. The final library comprised 50 unique sequences representing distinct tRNAs at the isodecoder level, alongside three control sequences for benchmarking purposes. Transfer RNA (tRNA) abundance estimation and normalization. The raw sequence files in FASTQ format were processed through the AQRNA-seq data analytical pipeline v 2.0 ( Extended Data Fig. 1 ; detailed in the section “Optimization of AQRNA-seq data analytical pipeline”) to derive the raw tRNA abundance matrix. To enable downstream comparisons across library preparation protocols while accounting for both sequencing depth and sample composition, the raw tRNA abundance was normalized using the median of ratios method, as implemented in the DESeq2 v 1.36.0 49 package in R. Efficacy assessment and comparison of PD mitigation approaches. The presence of PD in the constructed sequencing libraries leads to sequence reads lacking inserts or containing short inserts. Therefore, for each sample, we calculated the proportion of quality-filtered reads with different insert lengths by running Cutadapt v 4.2 43 iteratively, setting various thresholds for the minimum insert length. The efficacy of PD mitigation approaches was subsequently inferred from the proportion of reads with inserts longer than 10 bp. Specifically, the read proportion data was fitted with a linear regression model, using the stats v 4.3.2 50 package in R, wherein PD mitigation approach (GE vs ExoI) was included as the primary variable of interest and biological replicate as a blocking effect to account for potential experimental variations. To account for potential effects of PF mitigation approaches on read proportion or the efficacy of PD mitigation approaches, additive and interactive models including PF mitigation approach (AlkB vs SS) as an additional factor were also constructed, and the final model was determined based on F tests. Efficacy assessment and comparison of PF mitigation approaches. As PF eventually results in partial alignments of sequence reads with the reference sequences, the proportion of full-length reads could be used as a performance measure for assessing the efficacy of different PF mitigation approaches. Therefore, we calculated the proportion of full-length reads for individual tRNAs, with full-length reads defined as those aligned with the reference sequence from the 3’ to the 5’ ends, allowing for up to five missing bases at both ends. Subsequently, the full-length read proportion data was fitted with a linear regression model, using the stats v 4.3.2 50 package in R, wherein PF mitigation approach was included as the primary variable of interest and biological replicate as a blocking effect. Additive and interactive models including tRNA and PD mitigation approach as additional factors were also constructed, and the final model was determined based on F tests. In cases where interactions in model were deemed significant, post hoc comparisons were performed using the emmeans v 1.10.0 51 package in R, with the Benjamini-Hochberg (BH) method for multiple corrections. Quantitative comparison of tRNA landscapes across library preparation protocols. The concordance in the tRNA landscape obtained from different protocols was assessed by calculating the pairwise Pearson’s correlation coefficients (r) among the protocols based on the normalized abundance of individual tRNAs, averaged across biological replicates, using the PerformanceAnalytics v 2.0.4 52 package in R. To objectively rank the performance of the protocols to accurately depict the cellular tRNA landscape, we compared the normalized tRNA abundances from each protocol with an orthogonal tRNA abundance dataset 13 , which was derived using two-dimensional gel electrophoresis to separate tRNAs, followed by quantification using quantitative Northern blot. The quantitative performance of each protocol was inferred from the Pearson’s r between its quantitative results and the orthogonal dataset, with higher correlations denoting superior performance. Assessment and comparison of the effect of PF and PD mitigation approaches on the estimated abundance of tRNAs. To investigate whether variations in quantitative results across protocols could be attributed to the incorporation of different PF and PD mitigation approaches, the raw tRNA abundance data was transformed using the regularized logarithm (rlog) method implemented in the DESeq2 v 1.36.0 49 package in R, and hierarchical clustering was performed based on the Euclidean distance matrix derived from the rlog-transformed, normalized abundance data using the pheatmap v 1.0.12 53 package in R. Furthermore, principal component analysis (PCA) was performed on the same data using DESeq2 to discern whether any effects related to the use of different PF and PD mitigation approaches would result in distinct clusters among samples subject to different protocols. To further scrutinize the effect of PF and PD mitigation approaches on estimated abundance at the level of individual tRNAs, we followed the differential expression analysis workflow implemented in DESeq2 to identify tRNAs exhibiting differential estimated abundance as influenced by PF and/or PD mitigation approaches. Briefly, the workflow includes the following steps: (i) normalization of the raw abundance data using the median of ratios method; (ii) dispersion parameter estimation for individual tRNAs using empirical Bayes shrinkage; (iii) construction of a negative binomial regression model for each tRNA; (iv) estimation of logarithmic (log) fold changes using Approximate Posterior Estimation for generalized linear model 54 ; and (v) Wald tests for significant testing (with the BH method for multiple correction). In each negative binomial regression model, the normalized tRNA abundance was designated as the response variable and biological replicate was included as a blocking effect; both PF and PD mitigation approaches were included as primary factors of interest, along with the interaction term between these two factors. Position-specific percent coverage and identity calculation. To obtain a granular view on the alignment of reads with the reference sequences and to facilitate the mechanistic understanding of variations in quantitative results across protocols, we calculated position-specific percent coverage and identity for each reference sequence in each sample. Specifically, two FASTA files were generated for each reference sequence in each sample, containing the reference sequence and all its mapped reads, respectively. blastn v 2.6.0 45 was used to align the reads with the reference sequence, specifying SAM as the output format, followed by conversion of the output file into the BAM format using samtools view v 1.19.2 55 . The coverage at each position within a reference sequence across all aligned reads was extracted using samtools depth v 1.19.2 55 and subsequently divided by the maximum coverage within the reference sequence to derive the position-specific percent coverage. samtools mpileup v 1.19.2 55 was then used to extract the pileup results of the read alignments for each reference sequence, which were parsed to obtain the position-specific percent identity using a custom R script. For each position, this was calculated as the percentage of reads that contain the same base as the reference sequence. Mapping of tRNA modifications onto the reference sequences. The RF00005 alignments for all E. coli tRNAs available on Modomics 22 were downloaded. blastn v 2.6.0 40 was used to identify matches (100% query coverage; >98% percent identity) between our reference sequences used for mapping reads and the tRNAs present in the alignments. Reference sequences without a match with any tRNAs in the alignments were excluded from further analyses. For mapping tRNA modifications onto reference sequences, positions of tRNA modifications present in the alignments were exported from Modomics, specifying short name for the symbol filter to represent modification identities. These modifications were subsequently assigned to specific positions on the reference sequences by cross-referencing the original positions within the reference sequences and their corresponding positions in the alignments. Assessment of the impact of tRNA modifications on RT. As the interference of tRNA modifications with the RT process represents a predominant cause contributing to premature RT terminations, leading to suboptimal quantitative performance, we aimed to infer the adverse impact of these modifications on RT from the position-specific percent coverage data, assuming more adverse impacts would correlate with lower percent coverages. Due to the interdependence of percent coverages at a specific position with those closer to the 3’ end, it was difficult to assess the absolute impact of modifications at individual positions. Nevertheless, our data enabled a direct comparison of the impact of modifications at all possible positions across different protocols. To this end, the position-specific percent coverage data in the presence of tRNA modifications were obtained for each position within the alignments, and positions with insufficient data points (n < 3) for any of the protocols were precluded from further analyses. For each retained position, the data were fitted with a linear regression model, using the stats v 4.2.1 50 package in R, designating percent coverage as the response variable, protocol as the primary factor of interest, and biological replicate as a blocking effect. In cases where percent coverage at a specific position exhibited significant variations across different protocols, post hoc pairwise comparisons (i.e., Tukey’s HSD tests) were performed among protocols using the emmeans v 1.8.1.1 51 package in R. We next calculated the percent coverage resulting from all unique combinations of modification and position, subjected to different protocols, aiming to identify combinations imposing greatest challenges in RT and thus the accurate quantification of AQRNA-seq. For each protocol, we normalized the percent coverages to the maximum value at the corresponding position in the presence of modifications and subsequently classified the normalized percent coverages into ultra-high (75–100%), high (50–75%), low (25–50%), and ultra-low (0–25%) categories. For modification-position combinations resulting in ultra-low normalized percent coverages, the read alignments against the reference sequence of the corresponding tRNAs were manually examined to confirm their propensities to induce premature RT terminations. Assessment of the ability of AQRNA-seq in quantitatively mapping tRNA modifications. While quantitative mapping of tRNA modifications through NGS-based methods may depend on read pileups resulting from RT stops at modified nucleotides, our endeavors towards enhancing quantitative accuracy involve minimizing premature RT terminations, inherently deprecating reliance on this approach. However, certain tRNA modifications may prompt the insertion of incorrect residues as the RTase traverses them, resulting in unique mutation signatures (e.g., reduced percent identity and specific base profiles) that could be identified through subsequent analysis of sequencing data. Hence, to elucidate the potential of different protocols in the quantitative mapping of tRNA modifications, we calculated the average percent identity in the presence of each modification at every possible position, subjected to each of the protocols. Declarations Acknowledgements The authors thank the MIT BioMicro Center and its Director, Dr. Stuart Levine, for support and advice during the performance of the studies presented here. The authors gratefully acknowledge funding from the Singapore National Science Foundation under the Singapore-MIT Alliance for Research and Technology Antimicrobial Resistance Interdisciplinary Research Group (PCD), National Institutes of Health Transformative Award ES031576 (PCD), and Center grant P30-ES002109 from the National Institute of Environmental Health Sciences of the National Institutes of Health. Data availability Sequencing data are being deposited in the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) and will be available upon request prior to publication. Code availability The complete data analysis pipeline is available at GitHub (https://github.com/dedonlab/AQRNA-seq-method-optimization.git) Author contributions R.C., L.L., B.C., and P.C.D. designed the experiments. R.C., M.S.D., and P.C.D. wrote the first draft of the manuscript. R.C., D.Y., and L.L. performed experiments to develop and apply the sequencing method. R.C. developed the data processing pipeline and performed all statistical analyses. B.C. and P.C.D. oversaw study design and analysis. All authors contributed to the manuscript writing. References Wang, Z., Gerstein, M. & Snyder, M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet 10 , 57-63, (2009). Byron, S. A., Van Keuren-Jensen, K. R., Engelthaler, D. M., Carpten, J. D. & Craig, D. W. Translating RNA sequencing into clinical diagnostics: opportunities and challenges. Nat Rev Genet 17 , 257-271, (2016). Benesova, S., Kubista, M. & Valihrach, L. Small RNA-Sequencing: Approaches and Considerations for miRNA Analysis. Diagnostics (Basel) 11 , 964, (2021). Motorin, Y. & Helm, M. Methods for RNA Modification Mapping Using Deep Sequencing: Established and New Emerging Technologies. Genes (Basel) 10 , 35 (2019). Zhang, Y., Lu, L. & Li, X. Detection technologies for RNA modifications. Exp Mol Med 54 , 1601-1616, (2022). Zhang, L. S., Dai, Q. & He, C. Base-Resolution Sequencing Methods for Whole-Transcriptome Quantification of mRNA Modifications. Acc Chem Res 57 , 47-58, (2024). Cummings, B. B. et al. Improving genetic diagnosis in Mendelian disease with transcriptome sequencing. Sci Transl Med 9 , eaal5209, (2017). Haque, A., Engel, J., Teichmann, S. A. & Lonnberg, T. A practical guide to single-cell RNA-sequencing for biomedical research and clinical applications. Genome Med 9 , 75, (2017). Alon, S. et al. Barcoding bias in high-throughput multiplex sequencing of miRNA. Genome Res 21 , 1506-1511, (2011). Fuchs, R. T., Sun, Z., Zhuang, F. & Robb, G. B. Bias in ligation-based small RNA sequencing library construction is determined by adaptor and RNA structure. PLoS One 10 , e0126049, (2015). Pang, Y. L., Abo, R., Levine, S. S. & Dedon, P. C. Diverse cell stresses induce unique patterns of tRNA up- and down-regulation: tRNA-seq for quantifying changes in tRNA copy number. Nucleic Acids Res 42 , e170, (2014). Li, F. et al. Regulatory impact of RNA secondary structure across the Arabidopsis transcriptome. Plant Cell 24 , 4346-4359, (2012). Chen, R., Yim, D. & Dedon, P. C. AQRNA-seq for Quantifying Small RNAs. J Vis Exp , 2024 Feb 2;(204), (2024). Hu, J. F. et al. Quantitative mapping of the cellular small RNA landscape with AQRNA-seq. Nat Biotechnol 39 , 978-988, (2021). Dai, Q., Zheng, G., Schwartz, M. H., Clark, W. C. & Pan, T. Selective Enzymatic Demethylation of N(2) ,N(2) -Dimethylguanosine in RNA and Its Application in High-Throughput tRNA Sequencing. Angew Chem Int Ed Engl 56 , 5017-5020, (2017). Dong, H., Nilsson, L. & Kurland, C. G. Co-variation of tRNA abundance and codon usage in Escherichia coli at different growth rates. J Mol Biol 260 , 649-663, (1996). Rychlik, W. Selection of primers for polymerase chain reaction. Mol Biotechnol 3 , 129-134, (1995). Brownie, J. et al. The elimination of primer-dimer accumulation in PCR. Nucleic Acids Res 25 , 3235-3241, (1997). Illumina, Inc. Adapter dimers causes, effects, and how to remove them. https://knowledge.illumina.com/library-preparation/general/library-preparation-general-troubleshooting-list/000001911, (2023). New England Biolabs, Inc. Six Tips for a Perfect Gel Extraction. https://www.neb.com/en-us/tools-and-resources/usage-guidelines/six-tips-for-a-perfect-gel-extraction#, (2024). Keijzers, G., Liu, D. & Rasmussen, L. J. Exonuclease 1 and its versatile roles in DNA repair. Crit Rev Biochem Mol Biol 51 , 440-451, (2016). Cappannini, A. et al. MODOMICS: a database of RNA modifications and related information. 2023 update. Nucleic Acids Res 52 , D239-D244, (2024). Nordmann, P. L., Makris, J. C. & Reznikoff, W. S. Inosine induced mutations. Mol Gen Genet 214 , 62-67, (1988). Weber, M. et al. Chromosome-wide and promoter-specific analyses identify sites of differential DNA methylation in normal and transformed human cells. Nat Genet 37 , 853-862, (2005). Schmidt, D. et al. ChIP-seq: using high-throughput sequencing to discover protein-DNA interactions. Methods 48 , 240-248, (2009). Meers, M. P., Bryson, T. D., Henikoff, J. G. & Henikoff, S. Improved CUT&RUN chromatin profiling tools. Elife 8 , e46314, (2019). Bronner, I. F., Quail, M. A., Turner, D. J. & Swerdlow, H. Improved Protocols for Illumina Sequencing. Curr Protoc Hum Genet 80 , 18 12 11-42, (2014). Gao, X. et al. A reassessment of several erstwhile methods for isolating DNA fragments from agarose gels. 3 Biotech 11 , 138, (2021). Shore, S., Henderson, J. M. & McCaffrey, A. P. CleanTag Adapters Improve Small RNA Next-Generation Sequencing Library Preparation by Reducing Adapter Dimers. Methods Mol Biol 1712 , 145-161, (2018). Chen, Z. et al. Transfer RNA demethylase ALKBH3 promotes cancer progression via induction of tRNA-derived small RNAs. Nucleic Acids Res 47 , 2533-2545, (2019). Falnes, P. O. Repair of 3-methylthymine and 1-methylguanine lesions by bacterial and human AlkB proteins. Nucleic Acids Res 32 , 6260-6267, (2004). Cozen, A. E. et al. ARM-seq: AlkB-facilitated RNA methylation sequencing reveals a complex landscape of modified tRNA fragments. Nat Methods 12 , 879-884, (2015). Zucha, D., Androvic, P., Kubista, M. & Valihrach, L. Performance Comparison of Reverse Transcriptases for Single-Cell Studies. Clin Chem 66 , 217-228, (2020). Muramatsu, T. et al. Codon and amino-acid specificities of a transfer RNA are both converted by a single post-transcriptional modification. Nature 336 , 179-181, (1988). Wilson, R. K. & Roe, B. A. Presence of the hypermodified nucleotide N6-(delta 2-isopentenyl)-2-methylthioadenosine prevents codon misreading by Escherichia coli phenylalanyl-transfer RNA. Proc Natl Acad Sci U S A 86 , 409-413, (1989). Behrens, A., Rodschinka, G. & Nedialkova, D. D. High-resolution quantitative profiling of tRNA abundance and modification status in eukaryotes by mim-tRNAseq. Mol Cell 81 , 1802-1815, (2021). Fang, Z. et al. The Transcriptome-Wide Mapping of 2-Methylthio-N(6)-isopentenyladenosine at Single-Base Resolution. J Am Chem Soc 145 , 5467-5473, (2023). Bao, Z., Li, T. & Liu, J. Determining RNA Natural Modifications and Nucleoside Analog-Labeled Sites by a Chemical/Enzyme-Induced Base Mutation Principle. Molecules 28 , 1517, (2023). Debnath, T. K. & Xhemalce, B. Deciphering RNA modifications at base resolution: from chemistry to biology. Brief Funct Genomics 20 , 77-85, (2021). Kietrys, A. M., Velema, W. A. & Kool, E. T. Fingerprints of Modified RNA Bases from Deep Sequencing Profiles. J Am Chem Soc 139 , 17074-17081, (2017). National Academies of Sciences, Engineering, and Medicine. Charting a Future for Sequencing RNA and Its Modifications: A New Era for Biology and Medicine . https://www.ncbi.nlm.nih.gov/pubmed/39159274, (2024). Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34 , i884-i890, (2018). Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J 17 , 3, (2011). Zhang, J., Kobert, K., Flouri, T. & Stamatakis, A. PEAR: a fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics 30 , 614-620, (2014). Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J Mol Biol 215 , 403-410, (1990). Baba, T. et al. Construction of Escherichia coli K-12 in-frame, single-gene knockout mutants: the Keio collection. Mol Syst Biol 2 , 2006 0008, (2006). Chan, P. P. & Lowe, T. M. GtRNAdb: a database of transfer RNA genes detected in genomic sequence. Nucleic Acids Res 37 , D93-97, (2009). Chan, P. P. & Lowe, T. M. GtRNAdb 2.0: an expanded database of transfer RNA genes identified in complete and draft genomes. Nucleic Acids Res 44 , D184-189, (2016). Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15 , 550, (2014). R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Available online at https://www.R-project.org/, (2018). Lenth, R. emmeans: estimated marginal means, aka least-squares means. R package version 1.10.0 , https://CRAN.R-project.org/package=emmeans, (2024). Peterson, B. & Carl, P. PerformanceAnalytics: econometric tools for performance and risk analysis. R package version 2.0.4. , https://CRAN.R-project.org/package=PerformanceAnalytics, (2020). Kolde, R. pheatmap: pretty heatmaps. R package version 1.0.12. , https://CRAN.R-project.org/package=pheatmap, (2019). Zhu, A., Ibrahim, J. G. & Love, M. I. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics 35 , 2084-2092, (2019). Danecek, P. et al. Twelve years of SAMtools and BCFtools. Gigascience 10 , giab008, (2021). Additional Declarations There is NO Competing Interest. Supplementary Files ExtendedDataFigureslegends.docx ChenetalSupplementaryTable102025617.xlsx Supplementary Table 10 ChenetalSupplementaryTable22025617.xlsx Supplementary Table 2 ChenetalSupplementaryTable32025617.xlsx Supplementary Table 3 ChenetalSupplementaryTable42025617.xlsx Supplementary Table 4 ChenetalExtendedDataFigure12025615.pdf Extended Data Figure 1 ChenetalSupplementaryTable62025617.xlsx Supplementary Table 6 ChenetalExtendedDataFigure52025615.pdf Extended Data Figure 5 ChenetalSupplementaryTable82025617.xlsx Supplementary Table 8 ChenetalSupplementaryTable92025617.xlsx Supplementary Table 9 ChenetalExtendedDataFigure42025615.pdf Extended Data Figure 4 ChenetalSupplementaryTable112025617.xlsx Supplementary Table 11 ChenetalSupplementaryTable52025617.xlsx Supplementary Table 5 ChenetalExtendedDataFigure22025615.pdf Extended Data Figure 2 ChenetalExtendedDataFigure32025615.pdf Extended Data Figure 3 ChenetalSupplementaryTable72025617.xlsx Supplementary Table 7 ChenetalSupplementaryNotes2025714.pdf Supplementary Notes ChenetalSupplementaryTable12025617.xlsx Supplementary Table 1 Cite Share Download PDF Status: Under Review Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7256873","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":520844194,"identity":"7f98d042-439d-4576-a6eb-bf70d7ca95c7","order_by":0,"name":"Peter Dedon","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAqklEQVRIiWNgGAWjYHACAwYeBgY5BmYQm40ELcYMzMwkaklsYCBWizn74Y0P3lTcSV/bzn+A4UPZYcJaLHvSig3nnHmWu+0wMwPjjHNEaDG4wWMmzdt2GKyFGcggSov5b6DKdDOQlr9EajEDGZ4A1sJIjBaQXyTnnDlsCHSYwcGec+mEtYBC7MObisPyZucPPnzwo8yaCIchcw4QVo+uZRSMglEwCkYBVgAAXwo5pP9PfXoAAAAASUVORK5CYII=","orcid":"https://orcid.org/0000-0003-0011-3067","institution":"Massachusetts Institute of Technology","correspondingAuthor":true,"prefix":"","firstName":"Peter","middleName":"","lastName":"Dedon","suffix":""},{"id":520844195,"identity":"3678687b-a7ce-4753-ab25-d94025c6b2f6","order_by":1,"name":"Ruixi Chen","email":"","orcid":"","institution":"Massachusetts Institute of Technology","correspondingAuthor":false,"prefix":"","firstName":"Ruixi","middleName":"","lastName":"Chen","suffix":""},{"id":520844196,"identity":"15d48d3e-ba39-4c85-afe0-bc16f78f7cc5","order_by":2,"name":"Lili Liu","email":"","orcid":"","institution":"Shanghai Jiao Tong University","correspondingAuthor":false,"prefix":"","firstName":"Lili","middleName":"","lastName":"Liu","suffix":""},{"id":520844197,"identity":"d8b1f8dc-59ab-4e80-947f-56f46fa7be8a","order_by":3,"name":"Bo Cao","email":"","orcid":"","institution":"Shanghai Jiao Tong University","correspondingAuthor":false,"prefix":"","firstName":"Bo","middleName":"","lastName":"Cao","suffix":""}],"badges":[],"createdAt":"2025-07-31 01:00:27","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7256873/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7256873/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":93034407,"identity":"e8bb774b-3a74-449b-b0c3-9b7f440b6b65","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":12824817,"visible":true,"origin":"","legend":"","description":"","filename":"ChenetalText202589.docx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/89157055918d32cb45154b54.docx"},{"id":93035384,"identity":"6f271a27-9f61-4e2a-93f9-e6301eed6d2e","added_by":"auto","created_at":"2025-10-08 10:58:41","extension":"json","order_by":1,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":5482,"visible":true,"origin":"","legend":"","description":"","filename":"COMMSBIO258781T.json","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/e3231c4504a8cba1b07bf4fa.json"},{"id":93035391,"identity":"22a481dc-573c-4843-bea6-9eb87f786afa","added_by":"auto","created_at":"2025-10-08 10:58:42","extension":"pdf","order_by":7,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":226213,"visible":true,"origin":"","legend":"","description":"","filename":"ChenetalSupplementaryNotes2025714.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/0ebe5d96ebc0d70ea830e19e.pdf"},{"id":93035385,"identity":"f0af4842-e6d1-4149-904f-c1fa843b3a7d","added_by":"auto","created_at":"2025-10-08 10:58:41","extension":"xml","order_by":8,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":154966,"visible":true,"origin":"","legend":"","description":"","filename":"COMMSBIO258781T0enriched.xml","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/7a4d1af768142ecc9b76784c.xml"},{"id":93034419,"identity":"e42c1f6b-c07f-4915-9646-56be60d9a9c9","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"png","order_by":14,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":656298,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/f8149a1c65720205b4bb70db.png"},{"id":93034443,"identity":"e4389113-2d32-4854-aa45-458e664bccc6","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"emf","order_by":15,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":10820276,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage2.emf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/4c62fe98c7920a84e6413e07.emf"},{"id":93035390,"identity":"2c223425-ea58-42a5-b42a-231ee01a5c08","added_by":"auto","created_at":"2025-10-08 10:58:41","extension":"emf","order_by":16,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1664112,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage3.emf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/4afb89da6d4836e3a70081bf.emf"},{"id":93034445,"identity":"5a59cea9-f93b-483a-bf91-218625c017c4","added_by":"auto","created_at":"2025-10-08 10:50:43","extension":"emf","order_by":17,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1591432,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage4.emf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/94c953b0af726468b525bbfa.emf"},{"id":93034442,"identity":"effeb987-0d08-4a03-af57-3f3656af62c1","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"emf","order_by":18,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":5023964,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage5.emf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/9fa4f0f172de4123dbc5dd21.emf"},{"id":93034436,"identity":"e0dae84d-bf92-4b84-9ac7-31bb44b0b097","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"png","order_by":19,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":140397,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/f69cf7988858c4679e83f00b.png"},{"id":93035773,"identity":"406a0f2f-32f2-4734-ae25-31f053628863","added_by":"auto","created_at":"2025-10-08 11:06:42","extension":"png","order_by":20,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":33992,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/1d3c7a66552656fe6d5b5a61.png"},{"id":93035394,"identity":"9a2e839f-7e8e-481f-884a-d8028eb48696","added_by":"auto","created_at":"2025-10-08 10:58:42","extension":"png","order_by":21,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":18200,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/4d148fe320b1bd7cb05f3a12.png"},{"id":93035395,"identity":"3ee7f40a-abc8-49b8-a540-d65ba8ec540b","added_by":"auto","created_at":"2025-10-08 10:58:42","extension":"png","order_by":22,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":15350,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/e4cf464bbddc8c4535eb54c5.png"},{"id":93034444,"identity":"d5981dac-7c90-4988-8e93-356b95279b35","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"png","order_by":23,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":22606,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/e68375e39bd8193c7228d0d5.png"},{"id":93034438,"identity":"9934fac1-d7a6-4f12-993a-fb1c9eff5df0","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"xml","order_by":24,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":151325,"visible":true,"origin":"","legend":"","description":"","filename":"COMMSBIO258781T0structuring.xml","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/de57f20eb34c34ce82c3d53f.xml"},{"id":93034434,"identity":"fb6eda6c-cd3e-4cd5-8b7c-f8e652ad87c9","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"html","order_by":25,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":167404,"visible":true,"origin":"","legend":"","description":"","filename":"earlyproof.html","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/17d369e810e38fd66a556ed5.html"},{"id":93034396,"identity":"4c89485f-0107-4d2b-b016-ee40d3d3bf2a","added_by":"auto","created_at":"2025-10-08 10:50:40","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1322918,"visible":true,"origin":"","legend":"\u003cp\u003eSee image above for figure legend\u0026nbsp;\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/d0f534756837a3fed602b966.png"},{"id":93034412,"identity":"3d580c0f-320a-4c7a-ab86-c263ea6ecff0","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":2172245,"visible":true,"origin":"","legend":"\u003cp\u003eSee image above for figure legend\u0026nbsp;\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/d77c9c70d1bc13bdf3e1a8cf.png"},{"id":93034395,"identity":"7068097f-3711-49ff-8d21-9476f177022a","added_by":"auto","created_at":"2025-10-08 10:50:40","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":1063748,"visible":true,"origin":"","legend":"\u003cp\u003eSee image above for figure legend\u0026nbsp;\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/51c19683d7821adf01930c60.png"},{"id":93034402,"identity":"d58244f5-5ea9-4182-b918-5f4afdd86c29","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":992432,"visible":true,"origin":"","legend":"\u003cp\u003eSee image above for figure legend\u0026nbsp;\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/09bc8a7beffccddfe23ea421.png"},{"id":93035386,"identity":"b4c1ce58-2880-49ab-a542-09c28be2af79","added_by":"auto","created_at":"2025-10-08 10:58:41","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":1762851,"visible":true,"origin":"","legend":"\u003cp\u003eSee image above for figure legend\u0026nbsp;\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/c059ab8205f8aef81e2c6894.png"},{"id":93035388,"identity":"89c06af3-c43f-41c0-913d-a7ecebb8ecc5","added_by":"auto","created_at":"2025-10-08 10:58:41","extension":"eps","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":636886,"visible":true,"origin":"","legend":" Four variations of the AQRNQ-seq library preparation protocol (designated V1-V4) were developed to achieve a full-factorial design based on the original AQNRA-seq method (V1). Protocols V1 and V3 assess reduction in RT fall-off (PF) using AlkB treatment (ALKB) to protocols V2 and V4 using SuperScript IV (SS). Protocols V1 and V2 assess reduction in primer dimers (PD) using gel extraction (GE) compared to protocols V3 and V4 using exonuclease I treatment (EXOI). Gray boxes in the center column depict the protocol workflow. For each protocol variant, an arrow traverses a series of gray boxes indicating sequential steps. The left and right columns detail the strategies for mitigating PF and reducing PD implemented in different approaches.","description":"","filename":"drawingimage1.eps","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/5bb32dd0cd215cf9f87e0f27.eps"},{"id":93034413,"identity":"d43e5666-9b9f-41bc-8551-16d582068e40","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"eps","order_by":10,"title":"Figure 10","display":"","copyAsset":false,"role":"figure","size":332624,"visible":true,"origin":"","legend":"EXOI and SS enhance the quantitative performance of AQRNA-seq.\u0026nbsp; Size distribution of cDNA libraries following PCR amplification. Representative results from one biological replicate are shown. The size of primer dimer (PD) is 175 bp, while target fragments typically range from 200 bp (two primers\u0026thinsp;+\u0026thinsp;20 bp miRNA) to 300 bp (two primers\u0026thinsp;+\u0026thinsp;120 bp 5S rRNA) in size. Proportion of quality-filtered reads by insert length for samples shown in panel . Concentration of sequencing libraries following qPCR amplification. Data are summarized as mean (bar heights)\u0026thinsp;\u0026plusmn;\u0026thinsp;standard error (error bars) across 3 biological replicates. Compact letters are used to denote statistical significance, wherein protocols sharing at least one common letter are not significantly different from each other in concentration. Proportion of full-length reads for individual tRNAs for the four library preparation protocols (V1-V4). tRNAs highlighted in red showed significantly higher proportions in samples subjected to SS compared to ALKB. tRNAs labeled with a red star showed significantly higher proportions in EXOI compared to GE samples conditional to the use of SS, but not ALKB. Pairwise correlation for quantitative results across the four protocols. Correlation between the quantitative results from each protocol and those derived from a quantitative Northern blotting protocol.","description":"","filename":"drawingimage2.eps","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/aba99aefc2fc2d1920fd2763.eps"},{"id":93034424,"identity":"9178e6f4-4be5-452c-bac5-6bb051eb244f","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"eps","order_by":11,"title":"Figure 11","display":"","copyAsset":false,"role":"figure","size":166625,"visible":true,"origin":"","legend":" () Hierarchical clustering of samples based on regularized-logarithm (rlog)-transformed tRNA abundances. () Principal component analysis (PCA) scores plot based on the same data used in panel . () tRNAs showing differential estimated abundances due to the PF mitigation approach (left panel), PD mitigation approach (middle panel), and the interaction between PF and PD mitigation approaches (right panel). () Normalized abundance for tRNAs showing higher estimated abundances in SS compared to AlkB samples. () Normalized abundance for tRNAs showing lower estimated abundances in SS compared to AlkB samples. () Normalized abundance for tRNAs for which the interaction effect was deemed significant.","description":"","filename":"drawingimage3.eps","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/8a3a2a1818df01eed88068b7.eps"},{"id":93034415,"identity":"20eaddc5-2f65-4f32-a237-164f11fbfc09","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"eps","order_by":12,"title":"Figure 12","display":"","copyAsset":false,"role":"figure","size":170779,"visible":true,"origin":"","legend":" The read alignments against the reference sequence of Ile-GAT-1 subjected to protocols V1 and V4 are illustrated in panels and , respectively. The standardized tRNA sequence positions are displayed along the x-axis of the plots, with positions 0 and 76 representing the 5\u0026rsquo; and 3\u0026rsquo; termini, respectively. Anticodons are denoted with red stars. The exact sequence of Ile-GAT-1 is presented on the top of the plots. Gaps in the sequence, as well as the yellow shadings below, correspond to positions absent in Ile-GAT-1 and thus lacking associated data. When applicable, tRNA modifications are attached to the corresponding nucleotides within the sequence. The putative RT process, inferred from the position-specific percent coverage and identity data, is illustrated in the black frame atop the plots. Refer to the main text for detailed explanations.","description":"","filename":"drawingimage4.eps","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/249939a0b6c6844ee40eae98.eps"},{"id":93034426,"identity":"84854c36-7a0c-4ebd-a0c9-a1be9cb28e42","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"eps","order_by":13,"title":"Figure 13","display":"","copyAsset":false,"role":"figure","size":291992,"visible":true,"origin":"","legend":" () Summary of site-specific tRNA modifications for the 42 tRNAs. The graph depicts the cumulative occurrence of modifications (color-coded) at each position in a generic tRNA sequence. The inset defines the color coding for each modification along with the total number of occurrences of a modification (parentheses) at all positions across 42 tRNAs in the RF00005 alignments that were matched to our reference sequence library. () The percent coverage (see Methods for definition) in the presence of tRNA modifications at positions with sufficient data points (\u0026gt;\u0026thinsp;3) for all protocols. For each position, data are summarized as mean (bar heights)\u0026thinsp;\u0026plusmn;\u0026thinsp;standard error (error bars) across tRNAs. Compact letters are used to denote statistical significance, wherein protocols sharing at least one common letter are not significantly different from each other in percent coverage. () The normalized percent coverage in the presence of each modification at all possible positions. () The percent identity in the presence of each modification at all possible positions.","description":"","filename":"drawingimage5.eps","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/6305d69a0bf2487311c0c7ad.eps"},{"id":93036618,"identity":"6cb85a20-1adf-4bd4-8579-58cfaa954411","added_by":"auto","created_at":"2025-10-08 11:36:44","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":10795429,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/1e3aa36a-2981-469b-90ef-37c40993cd77.pdf"},{"id":93034393,"identity":"0c97830e-4fc5-4b4b-95eb-169965042ac5","added_by":"auto","created_at":"2025-10-08 10:50:40","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":14566,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cbr\u003e\u003c/p\u003e","description":"","filename":"ExtendedDataFigureslegends.docx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/a1916548d479d89b900fe4d7.docx"},{"id":93035383,"identity":"663878ed-529d-4779-8136-dad4c9dc197e","added_by":"auto","created_at":"2025-10-08 10:58:41","extension":"xlsx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":12028,"visible":true,"origin":"","legend":"Supplementary Table 10","description":"","filename":"ChenetalSupplementaryTable102025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/7809ba0e7496dcb0bfc5abd6.xlsx"},{"id":93035381,"identity":"e9b1e00c-7161-4e40-a0cd-fbe5652c33ad","added_by":"auto","created_at":"2025-10-08 10:58:40","extension":"xlsx","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":41004,"visible":true,"origin":"","legend":"\u003cp\u003eSupplementary Table 2\u003c/p\u003e","description":"","filename":"ChenetalSupplementaryTable22025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/c66b84f9d71ccc4a18e06538.xlsx"},{"id":93034404,"identity":"2d6ecca2-c932-436c-a3b5-8ff3ac0286ef","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"xlsx","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":16375,"visible":true,"origin":"","legend":"\u003cp\u003eSupplementary Table 3\u003c/p\u003e","description":"","filename":"ChenetalSupplementaryTable32025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/38e44550fa1a40983e8e10ee.xlsx"},{"id":93034409,"identity":"994b4889-e39e-4377-baaf-f3ba978d115a","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"xlsx","order_by":5,"title":"","display":"","copyAsset":false,"role":"supplement","size":13657,"visible":true,"origin":"","legend":"\u003cp\u003eSupplementary Table 4\u003c/p\u003e","description":"","filename":"ChenetalSupplementaryTable42025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/555dd0443225ad54467f8b39.xlsx"},{"id":93034417,"identity":"93812128-7b7a-4598-a2c9-3707f046a815","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"pdf","order_by":6,"title":"","display":"","copyAsset":false,"role":"supplement","size":341050,"visible":true,"origin":"","legend":"Extended Data Figure 1","description":"","filename":"ChenetalExtendedDataFigure12025615.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/ab746078e2dc6c4e87ac7e60.pdf"},{"id":93034399,"identity":"f6b0e62b-e385-435d-be5a-65ac8ec6acbc","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"xlsx","order_by":7,"title":"","display":"","copyAsset":false,"role":"supplement","size":1274946,"visible":true,"origin":"","legend":"Supplementary Table 6","description":"","filename":"ChenetalSupplementaryTable62025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/5efc88d499655ea052e9b40a.xlsx"},{"id":93034431,"identity":"044c3a67-1523-4c4c-8827-3a92cb86829c","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"pdf","order_by":8,"title":"","display":"","copyAsset":false,"role":"supplement","size":890923,"visible":true,"origin":"","legend":"Extended Data Figure 5","description":"","filename":"ChenetalExtendedDataFigure52025615.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/d335e5f689c27168a01ddbb7.pdf"},{"id":93034423,"identity":"d2924a99-16d5-43ca-9b40-3d5312dcffba","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"xlsx","order_by":9,"title":"","display":"","copyAsset":false,"role":"supplement","size":15978,"visible":true,"origin":"","legend":"Supplementary Table 8","description":"","filename":"ChenetalSupplementaryTable82025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/18ce2b098fbe1fc028e84b1c.xlsx"},{"id":93034421,"identity":"12922c09-5609-47b4-a05f-ceadcd472fef","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"xlsx","order_by":10,"title":"","display":"","copyAsset":false,"role":"supplement","size":12213,"visible":true,"origin":"","legend":"\u003cp\u003eSupplementary Table 9\u003c/p\u003e","description":"","filename":"ChenetalSupplementaryTable92025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/0d3f76152f75b81523bb0ea4.xlsx"},{"id":93034427,"identity":"1f527822-b745-475c-9446-e00f4cda12b2","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"pdf","order_by":11,"title":"","display":"","copyAsset":false,"role":"supplement","size":1127915,"visible":true,"origin":"","legend":"Extended Data Figure 4","description":"","filename":"ChenetalExtendedDataFigure42025615.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/6b01468c233de26543a370ea.pdf"},{"id":93034410,"identity":"718ac130-5ee3-43dd-8051-d03fc0d93a41","added_by":"auto","created_at":"2025-10-08 10:50:41","extension":"xlsx","order_by":12,"title":"","display":"","copyAsset":false,"role":"supplement","size":24655,"visible":true,"origin":"","legend":"\u003cp\u003eSupplementary Table 11\u003c/p\u003e","description":"","filename":"ChenetalSupplementaryTable112025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/98fcd8d557e168f0d1ee7e4a.xlsx"},{"id":93035392,"identity":"034626c7-c28d-49e2-ba40-f69570ba7d62","added_by":"auto","created_at":"2025-10-08 10:58:42","extension":"xlsx","order_by":13,"title":"","display":"","copyAsset":false,"role":"supplement","size":46045,"visible":true,"origin":"","legend":"\u003cp\u003eSupplementary Table 5\u003c/p\u003e","description":"","filename":"ChenetalSupplementaryTable52025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/9919316eb9955ee064329557.xlsx"},{"id":93035389,"identity":"ce6f93b0-7637-4d4a-be9f-e4b881e73bfc","added_by":"auto","created_at":"2025-10-08 10:58:41","extension":"pdf","order_by":14,"title":"","display":"","copyAsset":false,"role":"supplement","size":341418,"visible":true,"origin":"","legend":"\u003cp\u003eExtended Data Figure 2\u003c/p\u003e","description":"","filename":"ChenetalExtendedDataFigure22025615.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/9ec324af3e4f65733d225bef.pdf"},{"id":93034439,"identity":"a9410960-f8fa-4450-ac07-2924acfee343","added_by":"auto","created_at":"2025-10-08 10:50:42","extension":"pdf","order_by":15,"title":"","display":"","copyAsset":false,"role":"supplement","size":474792,"visible":true,"origin":"","legend":"\u003cp\u003eExtended Data Figure 3\u003c/p\u003e","description":"","filename":"ChenetalExtendedDataFigure32025615.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/fb5013be3fa91c9a4546bc6e.pdf"},{"id":93035387,"identity":"19a8434a-3353-4d80-abe9-502ed96e54c6","added_by":"auto","created_at":"2025-10-08 10:58:41","extension":"xlsx","order_by":16,"title":"","display":"","copyAsset":false,"role":"supplement","size":13050,"visible":true,"origin":"","legend":"\u003cp\u003eSupplementary Table 7\u003c/p\u003e","description":"","filename":"ChenetalSupplementaryTable72025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/68037d1feb5d324c6dd4b862.xlsx"},{"id":93035396,"identity":"f4133dfb-0022-4da3-ad74-b81ea5965867","added_by":"auto","created_at":"2025-10-08 10:58:42","extension":"pdf","order_by":17,"title":"","display":"","copyAsset":false,"role":"supplement","size":226213,"visible":true,"origin":"","legend":"\u003cp\u003eSupplementary Notes\u003c/p\u003e","description":"","filename":"ChenetalSupplementaryNotes2025714.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/5ce9a140af0e70267ec7b301.pdf"},{"id":93035397,"identity":"fecc0d9e-53dc-4fa2-ac91-1972532dc4c9","added_by":"auto","created_at":"2025-10-08 10:58:42","extension":"xlsx","order_by":22,"title":"","display":"","copyAsset":false,"role":"supplement","size":14901,"visible":true,"origin":"","legend":"Supplementary Table 1","description":"","filename":"ChenetalSupplementaryTable12025617.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7256873/v1/fd984f7d6db3e3a4932c894a.xlsx"}],"financialInterests":"There is \u003cb\u003eNO\u003c/b\u003e Competing Interest.","formattedTitle":"Gel-free library preparation for next-generation RNA sequencing and small RNA quantification","fulltext":[{"header":"Introduction","content":"\u003cp\u003eNext-generation sequencing (NGS) is a fast, sensitive, and cost-effective technology for rapid sequencing of nucleic acids, with recent applications enabling mapping of the \u0026gt;\u0026thinsp;170 chemical modifications of the epitranscriptome\u003csup\u003e\u003cspan additionalcitationids=\"CR2 CR3 CR4 CR5\" citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e. Initially used for DNA sequencing tasks such as whole-genome sequencing, targeted sequencing, epigenetics, and metagenomics, NGS was subsequently adapted to sequencing all types of RNA (RNA-seq)\u003csup\u003e\u003cspan additionalcitationids=\"CR2\" citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e. Advancements in sequencing instrumentation, coupled with the development of specialized library preparation methods and data analytical pipelines, have broadened the scope of RNA-seq applications beyond measuring gene expression, such as identifying novel transcripts, mutations, and alternative splicing\u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u003c/sup\u003e, as well as quantitative mapping of post-transcriptional modifications\u003csup\u003e\u003cspan additionalcitationids=\"CR5\" citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e and single-cell resolution RNA-seq (scRNA-seq)\u003csup\u003e\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eUntil recently, one notable limitation to most RNA-seq methods has been the inability to accurately measure the absolute abundance of RNA molecules within a cellular pool. Biased ligation of sequencing adaptors\u003csup\u003e\u003cspan additionalcitationids=\"CR10\" citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e and premature termination of reverse transcription (RT), due to polymerase fall-offs (PF) caused by post-transcriptional RNA modifications or secondary structures\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e,\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e\u003c/sup\u003e, are among key challenges to quantitative accuracy. These challenges, however, were addressed in the development of \u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eA\u003c/span\u003ebsolute \u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003eQ\u003c/span\u003euantification RNA-seq (AQRNA-seq)\u003csup\u003e\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e,\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e, a specialized RNA-seq method for quantitative profiling of all small RNA species (\u0026lt;\u0026thinsp;200 nt) in any cell, tissue, and organism. Quantitative accuracy is achieved by using custom DNA linkers to minimize ligation biases, a two-step ligation approach to fully capture truncated complementary DNAs (cDNAs), and an optional AlkB demethylation\u003csup\u003e\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e step to reduce PF by partial enzymatic removal of several types RNA methylation, thereby increasing the proportion of full-length cDNAs. As a result, AQRNA-seq achieves linearity between sequence read count and biological copy number, as validated in studies involving a reference library of 963 miRNAs, pooled RNA oligonucleotide standards with variable lengths, and an orthogonal tRNA abundance dataset derived from quantitative Northern blotting\u003csup\u003e\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eDespite the superior quantitative performance of AQRNA-seq compared to other methods, it suffered like most other NGS-based RNA-seq approaches from limitations posed by library preparation workflow and sensitivity. For example, all NGS-based RNA-seq methods suffer from highly abundant \u0026ldquo;primer dimers\u0026rdquo; (PD), which are formed by ligation of the 5' DNA linker to the 3' end of RT primers during library preparation. Due to their small sizes, PD can form clusters on the flow cell more efficiently than target RNA molecules, thus consuming the sequencing capacity and lowering sensitivity by producing unusable reads\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e,\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e. The traditional approach to reducing the impact of PD is to resolve the cDNA products by electrophoresis followed by gel extraction (GE) to size-select target fragments\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. However, GE suffers from highly variable accuracy of manual gel excision and inefficient extraction of target fragments\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e, which collectively reduce sensitivity. More importantly, the time-consuming manual labor involved in GE renders library preparation for AQRNA-seq and other RNA-seq methods unsuitable for automation and high-throughput applications. In addition to PD concerns, the quantitative accuracy of AQRNA-seq can be limited by high levels of truncated cDNAs persisting for many tRNAs\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e, especially those containing modifications insensitive to AlkB treatment. While being fully captured by the two-step ligation approach, these truncated cDNAs may fail to unambiguously align with the reference sequence library.\u003c/p\u003e\u003cp\u003eWe have now solved these problems by combining (1) a novel approach to PD mitigation involving digestion of excess RT primers using \u003cem\u003eE. coli\u003c/em\u003e exonuclease I\u003csup\u003e21\u003c/sup\u003e immediately after RT and (2) a highly processive reverse transcriptase (RTase), SuperScript IV, that reduces PF to increase full length reads. A full factorial experimental design (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e) was conducted to contrast these novel approaches with their original counterparts. Our results suggest that exonuclease I is significantly more efficient in reducing PD compared to GE, increasing RNA-seq sensitivity while in the meantime obviating human involvement in gel extraction. Similarly, SuperScript IV proved superior to AlkB in mitigating PF, especially in the presence of diverse RNA modifications, enhancing the quantitative accuracy of AQRNA-seq.\u0026nbsp;The use of both exonuclease I and SuperScript IV can be readily integrated into any RNA-seq workflow, allowing full automation of RNA-seq library preparation.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e"},{"header":"RESULTS","content":"\u003cp\u003e\u003cb\u003eEXOI enhances sensitivity and enables automation of RNA-seq workflows.\u003c/b\u003e The presence of PD in sequencing libraries, along with the need for gel extraction for removing PD, substantially deteriorates the performance of sequencing methods, particularly those targeting small RNAs. While gel extraction is commonly used for PD mitigation during library preparation, it imposes constraints on method sensitivity and automation. Here, we introduce EXOI as a simple, cost-effective, and universally applicable approach that substantially outperforms GE in efficiency and holds potential for enhancing sensitivity and facilitating automation, as well as high-throughput implementation, across all RNA-seq methods. In AQRNA-seq, ligation of 5\u0026rsquo; DNA linkers to the 3\u0026rsquo; end of RT primers leads to the formation of PD\u003csup\u003e13\u003c/sup\u003e. Hence, the ExoI nuclease is added immediately following RT to digest excess RT primers, thereby inhibiting the formation of PD. Electrophoretic analysis (3% agarose gel) of the cDNA libraries subsequent to PCR amplification revealed a lower density of fragments below the 200-bp size marker for samples subjected to protocols V3 and V4 (incorporating EXOI), as compared to protocols V1 and V2 (incorporating GE), respectively, confirming the efficacy of EXOI in inhibiting the formation of PD (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea). Furthermore, the proportion of quality-filtered reads with an insert length of \u0026gt;\u0026thinsp;10 bp was significantly higher (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05) for samples subjected to EXOI as compared to GE (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb, \u003cb\u003eSupplementary Table\u0026nbsp;1\u003c/b\u003e), suggesting a superior performance of EXOI in reducing PD and short inserts. More importantly, our results showed that the incorporation of EXOI substantially increased the concentration of cDNA libraries following qPCR amplification, exceeding that achieved with GE by over 11-fold (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec). This highlights the potential of EXOI to significantly enhance the sensitivity of AQRNA-seq.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eFinally, substituting GE with EXOI eliminates the need for human intervention, thus enabling the future automation and high-throughput implementation of AQRNA-seq.\u003c/p\u003e\u003cp\u003e\u003cb\u003eSS demonstrates potential for enhancing quantitative accuracy through its superior ability in attaining full-length reads.\u003c/b\u003e Having addressed challenges related to sensitivity and automation, we proceeded to optimize the quantitative accuracy of AQRNA-seq, which is contingent upon precise alignments of sequence reads against the reference library. In AQRNA-seq, the synthesis of truncated cDNAs due to premature RT terminations inevitably compromises the analysis and could be detrimental for differentiating between highly similar sequences. While ALKB can be used reduce PF, one cause of premature RT terminations, thereby increasing proportions of full-length reads, its efficacy is limited to a few methyl modifications, leaving most other modifications unaffected. Therefore, we assessed SS as an alternative approach for PF mitigation, which employs SuperScript IV, a high-processivity RTase. We identified an ascending trend in the proportion of full-length reads for individual tRNAs across protocols, with the mean proportion progressively increasing from 0.11 for V1, to 0.14 for V3, then to 0.35 for V2, and finally to 0.43 for V4 (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ed, \u003cb\u003eSupplementary Table\u0026nbsp;2\u003c/b\u003e). This suggests the potential of both SS and EXOI in increasing proportions of full-length reads, with SS showing the predominant influence. A linear regression model was constructed to assess the effect of PF and PD mitigation approaches, tRNA identity, and their potential interactions on full-length read proportions. While all main factors and interactions included in the model demonstrated significance (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05) based on F tests (\u003cb\u003eSupplementary Table\u0026nbsp;2\u003c/b\u003e), interpretation was only deemed appropriate for the interaction among the three main factors. Specifically, post hoc comparisons were performed to compare the model-reported estimated marginal means (EMM) of the proportion between ALKB and SS samples conditional to each PD mitigation approach and each tRNA (\u003cb\u003eSupplementary Table\u0026nbsp;2\u003c/b\u003e). Our results revealed a significantly (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05) higher proportion of full-length reads resulted from SS compared to ALKB for 39 tRNAs, regardless of the PD mitigation approach used (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ed). Importantly, for 7 of these tRNAs, the difference between EXOI and GE was significant conditional to the incorporation of SS, but not ALKB, showcasing the interactive effect between PF and PD mitigation approaches. Therefore, these results confirm the superior performance of SS over ALKB in attaining full-length reads, particularly in conjunction with EXOI, suggesting its potential for enhancing the quantitative accuracy of AQRNA-seq.\u003c/p\u003e\u003cp\u003e\u003cb\u003eBoth SS and EXOI resulted in variations in quantitative results across different protocols.\u003c/b\u003e Our findings suggest that SS and EXOI enhance the quantitative accuracy and sensitivity of AQRNA-seq.\u0026nbsp;We next sought to determine whether these enhancements would manifest as variations in quantitative results. Pairwise Pearson\u0026rsquo;s correlation coefficients (r) based on normalized tRNA abundances revealed that the incorporation of SS resulted in a lower correlation (r\u0026thinsp;=\u0026thinsp;0.84) with the original AQRNA-seq protocol compared to the incorporation of EXOI (r\u0026thinsp;=\u0026thinsp;0.96) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ee, \u003cb\u003eSupplemental Table\u0026nbsp;3\u003c/b\u003e). The lowest pairwise correlation was observed between the original protocol and the most revised version, which incorporates both SS and EXOI. To objectively rank the quantitative performance of the protocols, we assessed the correlation between the percentages of tRNAs within the tRNA pool obtained from each protocol and those derived from a quantitative Northern blotting protocol\u003csup\u003e\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u003c/sup\u003e, considering higher correlations with this orthogonal dataset were indicative of more accurate quantification. Our results showed that both SS and EXOI contributed to enhanced correlations with the orthogonal dataset, with the highest correlation observed for protocol V4 (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ef, \u003cb\u003eSupplemental Table\u0026nbsp;4\u003c/b\u003e). This suggests the potential of both SS and EXOI to enhance the quantitative performance of the original protocol, establishing V4 as the optimal protocol for the library preparation of AQRNA-seq.\u003c/p\u003e\u003cp\u003e\u003cb\u003eThe estimated abundance for individual tRNAs is more dependent on PF mitigation approach than PD mitigation approach.\u003c/b\u003e Having demonstrated that both SS and EXOI collectively led to an enhancement in quantitative performance, we next proceeded to assess and compare their respective impacts on the quantification of individual tRNAs. While our results revealed consistent quantification across protocols for most tRNAs, noticeable variations in estimated abundance were observed for selected tRNAs, such as Ile-GAT-1 and Phe-GAA-1 (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea, \u003cb\u003eSupplementary Table\u0026nbsp;5\u003c/b\u003e). Hierarchical clustering based on rlog-transformed and normalized tRNA abundance formed two distinct clusters of samples by PF mitigation approach (ALKB vs SS), while not by PD mitigation approach (GE vs EXOI). Consistent with this, projecting the samples onto the principal component analysis (PCA) score plot unveiled a more pronounced segregation of samples by PF as compared to PD mitigation approach (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb), further indicating a greater impact of SS relative to EXOI on estimated abundances. Importantly, the PCA score plot suggested an interaction effect between PF and PD mitigation approaches, as the separation distance between GE and EXOI samples subjected to SS (comparing V2 and V4 samples) was considerably larger than those subjected to ALKB (comparing V1 and V3 samples). To gain mechanistic understanding into the variation in estimated abundance attributable to SS and EXOI, we identified individual tRNAs exhibiting differential estimated abundance (|log\u003csub\u003e2\u003c/sub\u003e fold change| \u0026gt;\u0026thinsp;1 along with \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05) across protocols as affected by the incorporation of SS and/or EXOI (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ec, \u003cb\u003eSupplementary Table\u0026nbsp;5\u003c/b\u003e). While the estimated abundance remained consistent across all protocols for most tRNAs, corresponding well with the high correlations observed among the protocols, we identified 5 tRNAs (Arg-CCG-1, Ile-GAT-1, Ile2-CAT-1, Leu-TAA-1, and Val-GAC-1) that showed significantly higher estimated abundances in SS compared to ALKB samples (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ed), and 3 tRNAs (fMet-CAT-1, Leu-CAG-2, and Pro-CGG-1) showing lower estimated abundances (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ee). In addition, although none of the tRNAs were differentially estimated between GE and EXOI samples, the interaction effect was deemed significant for 6 tRNAs (Met-CAT-1, Ile2-CAT-2, Phe-GAA-1, Phe-GAA-2, Thr-CGT-1, and Thr-CGT-2), suggesting that differences in estimated abundance between ALKB and SS samples were further contingent upon PD mitigation approaches (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ef).\u003c/p\u003e\u003cp\u003e\u003cb\u003eSS showed the potential to correct underestimated tRNA abundances by ALKB.\u003c/b\u003e The most notable difference in estimated abundance between ALKB and SS samples was observed for Ile-GAT-1 (SS/ALKB\u0026thinsp;=\u0026thinsp;13.6), followed by Ile2-CAT-1 (SS/ALKB\u0026thinsp;=\u0026thinsp;5) and Val-GAC-1 (SS/ALKB\u0026thinsp;=\u0026thinsp;5) (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ed, \u003cb\u003eSupplementary Table\u0026nbsp;5\u003c/b\u003e). To elucidate mechanisms driving these differences, we modeled the RT process of the corresponding tRNAs by cross-referencing the position-specific percent coverage and identity data generated in this study (\u003cb\u003eSupplementary Table\u0026nbsp;6\u003c/b\u003e) with the positional information of tRNA modifications available on Modomics\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e\u003c/sup\u003e. This is illustrated in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e for Ile-GAT-1, for which the normalized read counts were 1,803 and 35,884 following protocols V1 and V4, respectively, suggesting a 20-fold difference (\u003cb\u003eSupplementary Table\u0026nbsp;5\u003c/b\u003e). Regardless of the protocols used, our results suggested the presence of a highly abundant 5\u0026rsquo; tRNA fragment, as most read alignments (3\u0026rsquo; to 5\u0026rsquo; direction) initiated at position 59, rather than the 3\u0026rsquo; end (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea,b). The RT of this fragment, however, was thwarted by the presence of 3-(3-amino-3-carboxypropyl)uridine (acp\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003eU) at position 45B, manifested as the sharp downward spikes (i.e., reduced percent coverage and identity) shown in different panels of Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e. When protocol V1 was used, there was a pronounced reduction in percent coverage at position 45B, followed by considerably reduced percent coverage from position 45A onwards, a hallmark of PF at acp\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003eU causing premature RT terminations (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea). Consequently, the cDNAs synthesized for the tRNA fragment were likely too short (i.e., ~\u0026thinsp;12 nt from position 45C to 59) to produce sequence reads that could be aligned with the reference sequence, leading to an underestimation in abundance. Conversely, following protocol V4, while dropping drastically at position 45B, the percent coverage returned to nearly 100% for positions past 45B, indicating that the high-processivity RTase copied past acp3U, albeit causing a -1 mutation, to create longer cDNAs that could be aligned with the reference sequence (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb). This critical difference in how tRNA modifications may adversely impact the RT process remained consistent for Ile2-CAT-1 and Val-GAC-1, wherein acp\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003eU triggered considerable premature RT termination with V1 but to a far lesser extent with V4 (\u003cb\u003eExtended Data\u003c/b\u003e Figs.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea,b and \u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea,b). Therefore, our findings illustrate the capability for SS to correct underestimation of tRNA abundances observed in samples processed with the ALKB method.\u003c/p\u003e\u003cp\u003eA notable feature of our analysis of Ile-GAT-1 was a marked reduction in percent identity around position 69 for Ile-GAT-1 (Figs.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e), suggesting the presence of a modification that has not been annotated in the Modomics database or a highly stable RNA secondary structure. This Putative modification/2\u0026deg; structure, potentially in combination with the pseudouridine (Y) at position 65, was correlated with near-stoichiometric PF in SS samples, which was partially mitigated by the ALKB method, as demonstrated by the higher percent coverages observed from positions 76 to 60 in ALKB samples (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea) compared to SS (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb). Thus, our findings not only suggest the potential for AQRNA-seq to unearth previously unidentified RNA modifications, along with their associated mutation signatures, but also highlight situations where ALKB may outperform SS in handling specific modifications.\u003c/p\u003e\u003cp\u003e\u003cb\u003eThe impact of tRNA modifications on the accuracy of tRNA quantification varied by protocol, location, and chemical structure.\u003c/b\u003e Since uninterrupted RT was demonstrated to be important for accurately quantifying individual tRNAs, we next characterized the potential of all modifications present on \u003cem\u003eE. coli\u003c/em\u003e tRNAs to interfere with RT for each protocol, aiming to provide a systematic overview of the effect of each modification on quantification of the tRNA. We identified 42 tRNA isoacceptors in the RF00005 alignments (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://rfam.org/family/RF00005\u003c/span\u003e\u003cspan address=\"https://rfam.org/family/RF00005\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) that were successfully matched to reference sequences in our library and contain 28 unique modifications located at 24 positions in the alignments (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e5\u003c/span\u003ea, \u003cb\u003eSupplementary Table\u0026nbsp;7\u003c/b\u003e). While most modifications tend to cluster in proximity within the same structural domain, the modification Y, as one of the most frequent modifications, spans across multiple domains. The impact of tRNA modifications on quantification was inferred from the position-specific percent coverage data, with lower percent coverages inversely correlated with more adverse impacts. We found that protocols V1-V3 resulted in insufficient coverage (i.e., \u0026le; 100 reads) at the 5\u0026rsquo;-end for 8, 7, and 7 of 42 tRNAs, respectively, consistent with premature RT terminations and/or low sensitivity resulting in reduced read counts for these tRNAs (\u003cb\u003eExtended Data\u003c/b\u003e Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e, \u003cb\u003eSupplementary Table\u0026nbsp;6\u003c/b\u003e). Conversely, protocol V4 maintained sufficient coverage (i.e., \u0026gt;\u0026thinsp;100 reads) toward the 5\u0026rsquo;-end of all tRNAs, suggesting limited premature RT termination due to SS for PF mitigation and increased sensitivity conferred by EXOI for PD removal. To systematically assess and compare the impact of modifications on quantitation across the four protocols, we analyzed the percent coverage data at each position with modified residues using linear regression, focusing only on tRNAs carrying modifications at the specific position. Due to insufficient coverage at selected positions across different protocols, the analysis was conducted for 19 positions with a minimum of three data points for each protocol. The percent coverage differed significantly across protocols at 18 of these positions (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e5\u003c/span\u003eb, \u003cb\u003eSupplementary Table\u0026nbsp;8\u003c/b\u003e). Post hoc pairwise comparisons revealed a significantly higher percent coverage for V2 and V4 compared to V1 and V3 for all 18 positions, confirming the superior capability of SS over ALKB to prevent premature PF at modifications. Consequently, modifications imposed less of an adverse impact on tRNA quantification with protocols incorporating SS compared to protocols incorporating ALKB. Notably, a significantly higher percent coverage was observed for V4 compared to V2 at 7 positions and for V3 compared to V1 at 6 positions. Again, this suggests an important role for EXOI in further enhancing quantitative accuracy of the protocols, likely due to its ability to enhance sensitivity, especially for tRNAs with low abundance.\u003c/p\u003e\u003cp\u003eFor each protocol, substantial variations in percent coverage were evident at certain positions across tRNAs with modifications, suggesting that the presence of different modifications at a specific position may have differential impacts on quantification. Therefore, we further scrutinized the percent coverage in the presence of individual modifications at all possible positions. In general, the percent coverage, normalized to the maximum value at the corresponding position in the presence of modifications, may exhibit substantial variability both across modifications at a specific position and across positions where a specific modification may localize (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e5\u003c/span\u003ec). For instance, among the 20 tRNAs carrying 11 distinct modifications at the wobble position 34, the normalized percent coverage averaged for each modification ranged from 0.8% to 73% for V1, 5% to 93% for V2, 1% to 80% for V3, and 2% to 96% for V4. The normalized percent coverage was subsequently classified into four categories: ultra-high (75\u0026ndash;100%), high (50\u0026ndash;75%), low (25\u0026ndash;50%), and ultra-low (0\u0026ndash;25%), wherein instances falling with the ultra-low category were further investigated, as they likely presented the greatest challenges in the quantification of AQRNA-seq (\u003cb\u003eSupplementary Table\u0026nbsp;9\u003c/b\u003e). Across different protocols, V4 resulted in the highest proportion of modification-position combinations (26 out of 44) exhibiting ultra-high normalized percent coverages while the lowest proportion of combinations (7 out of 44) showing ultra-low or insufficient data. In particular, our results revealed modification-position combinations manifesting high/ultra-high normalized percent coverages with protocol V4 but ultra-low normalized percent coverages with at least one of the other protocols, including N4-acetylcytidine (ac4C) at position 34, inosine (I) at position 34, 1-methylguanosine (m1G) at position 37, and 2-thiocytidine (s2C) at position 32. Conversely, all modification-position combinations manifesting ultra-low normalized percent coverages with V4 also fell into the ultra-low category with all the other protocols. This reaffirms the superior performance of protocol V4 as compared to the other protocols in minimizing the adverse impacts of modifications on quantification, corresponding well with our findings at the positional level. Consequently, further inferences in this regard were made exclusively based on results from protocol V4. While we identified two modification-position combinations (5-methyluridine [m5U] at position 53 and Y at position 54) showing insufficient data, both instances were associated with a minor tRNA SeC-TCA-1, thereby possibly attributable to the low abundance of this tRNA. Furthermore, we identified five modification-position combinations exhibiting ultra-low normalized percent coverages, including 2'-O-methylcytidine (Cm) at position 34, 5-carboxymethylaminomethyl-2'-O-methyluridine (cmnm5Um) at position 34, 2-lysidine (k2C) at position 34, 2-methylthio-N6-isopentenyladenosine (ms2i6A) at position 37, and Y at position 32. However, upon further examination of the read alignments of the corresponding tRNAs, it became evident that the observed ultra-low normalized percent coverages of Cm, cmnm5Um, and Y could be attributed to the presence of ms2i6A at position 37 of the corresponding tRNAs, which precipitated a sharp decline in percent coverage, followed by diminished percent coverages thereafter (\u003cb\u003eExtended Data\u003c/b\u003e Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). Hence, our findings suggested a high likelihood for ms2i6A and k2C, though not necessarily the other modifications, to induce premature PF and thereby potentially impose detrimental effect on tRNA quantification.\u003c/p\u003e\u003cp\u003e\u003cb\u003eQuantitative mapping of tRNA modifications varies among the protocols.\u003c/b\u003e One consequence of minimizing PF to increase the accuracy and sensitivity of small RNA quantification was a reduced ability to quantitatively map tRNA modifications through position-specific read pileups. However, quantitative modification mapping of tRNA modifications is still possible for those that cause the RT to insert the wrong nucleotide to create a mutation. We tested this approach to modification mapping by quantifying the position-specific percent sequence identity (i.e., 1-mutation frequency) at all possible positions with known modifications in \u003cem\u003eE. coli\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e5\u003c/span\u003ed). The percent identities fell below 90% for 2, 5, 2, and 6 modifications subjected to protocols V1, V2, V3, and V4, respectively, suggesting the potential for quantitative mapping of these tRNA modifications. For example, the percent identity ranged from 0.2% to 1.7% across all protocols for wobble I, consistent with its well-documented ability to induce a T-to-C mutation during RT\u003csup\u003e23\u003c/sup\u003e. In addition, the percent identity was notably lower with protocols V2 and V4 as compared to protocols V1 and V3 for acp3U, ms2i6A, m1G, and cmnm5Um. These results point to the utility of EXOI/SS-modified AQRNA-seq quantitative mapping of several modifications, with the potential for identifying predictive mutation signatures by base profile surrounding target positions.\u003c/p\u003e"},{"header":"DISCUSSION","content":"\u003cp\u003eHere we optimized the highly accurate AQRNA-seq\u003csup\u003e\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e,\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e method to improve quantitative accuracy and sensitivity, while obviating the gel purification step to boost efficiency and enable automation for greater throughput and scalability. These advancements were achieved with two enzymatic changes: (i) adding a step for post-RT digestion of excess primers using \u003cem\u003eE. coli\u003c/em\u003e exonuclease I for effective PD mitigation, and (ii) use of SuperScript IV\u0026mdash;a high-processivity RTase\u0026mdash;to circumvent RNA modifications during cDNA synthesis for improved PF mitigation. These changes substantially outperformed their original, labor-intensive counterparts of gel purification for PD mitigation and use of AlkB demethylase to remove a subset of RNA methylations for PF mitigation. Both EXOI and SS can be applied to almost any NGS RNA-seq method to facilitate automation.\u003c/p\u003e\u003cp\u003eGE is a routinely used approach in NGS RNA-seq library preparation to reduce the impact of PD in sequencing libraries and is implemented in various workflows such as CUT\u0026amp;RUN, ChIP-Seq, MeDIP-Seq, and the original AQRNA-seq\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e,\u003cspan additionalcitationids=\"CR25 CR26\" citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e\u003c/sup\u003e. However, its efficiency varies depending on sample characteristics, experimental conditions, gel cutting techniques, and operator experience, frequently leading to reduced sensitivity due to the loss of target fragments\u003csup\u003e\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e\u003c/sup\u003e. Furthermore, the labor-intensive nature of GE hinders the automation of NGS RNA-seq methods, thus constraining their scalability and throughput. To tackle this issue, we tested EXOI\u0026mdash;post-RT digestion of excess RT primers using ExoI nuclease (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e)\u0026mdash;as a more efficient and cost-effective alternative for mitigating PD. Notably, compared to GE, EXOI substantially enhances the sensitivity of AQRNA-seq through a reduced loss of target fragments and a more effective mitigation of PD (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb,c). Furthermore, EXOI eliminates the need for human involvement, thus facilitating automation and enabling higher-throughput small RNA quantification. While alternative approaches for PD mitigation are available, such as bead-based size selection (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.cytivalifesciences.com/en/us/news-center/better-data-via-size-selection-10001?ssp=1\u0026amp;setlang=en\u0026amp;cc=US\u0026amp;safesearch=moderate\u003c/span\u003e\u003cspan address=\"https://www.cytivalifesciences.com/en/us/news-center/better-data-via-size-selection-10001?ssp=1\u0026amp;setlang=en\u0026amp;cc=US\u0026amp;safesearch=moderate\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) and chemically modified adaptors\u003csup\u003e\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e\u003c/sup\u003e, they have significant drawbacks that limit their use. For instance, bead-based size selection lacks sufficient resolution for small size differences between PD and target fragments (e.g., 20 nt for miRNAs) and may not be applicable to low-input samples\u003csup\u003e\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e. The incorporation of chemically modified adaptors, on the other hand, requires the case-by-base redesign of DNA linkers. In contrast, the simplicity of EXOI makes it compatible and easily adaptable to most library preparation workflows. Moving forward, the incorporation of EXOI in AQRNA-seq and related methods will likely bolster the speed and accuracy of quantitative investigations into all forms of biological targets.\u003c/p\u003e\u003cp\u003ePost-transcriptional modifications are a significant cause of PF during RT, which reduces quantitative accuracy due to imprecise alignments of reads with the reference sequence library. While the ALKB method can be used to mitigate PF by enzymatically removing several simple methyl modifications such as m\u003csup\u003e6\u003c/sup\u003eA, m\u003csup\u003e1\u003c/sup\u003eA, and m\u003csup\u003e3\u003c/sup\u003eC\u003csup\u003e\u003cspan additionalcitationids=\"CR31\" citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e, most other RNA modifications not affected by this approach and demethylation efficiency is ribonucleoside-specific and highly variable due to the oxygen sensitivity of the Fe-S cluster of AlkB\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e. PF was an important contributor to the numerous truncated cDNAs synthesized for many tRNAs in the original AQRNA-seq protocol with ALKB\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e. To further enhance the quantitative accuracy of AQRNA-seq, we tested SS as an alternative approach for PF mitigation, employing a high-processivity RTase to circumvent diverse RNA modifications during cDNA synthesis\u003csup\u003e\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e\u003c/sup\u003e. Our results show that SS significantly increases the proportion of full-length reads for nearly 80% of tRNAs, contributing to a substantial enhancement in quantitative accuracy (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). Cross-referencing our sequencing results with the positional information of tRNA modifications available on Modomics\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e\u003c/sup\u003e provided insights into the superior performance of SS in averting PF due to a broad spectrum of modifications at various locations. For instance, while acp\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003eU at position 45B tended to induce PF in ALKB samples, it was bypassed during RT in SS samples, resulting in either a -1 mutation or a point mutation due to mismatched nucleotide incorporation, both of which are useful for modification mapping purposes. On a broader scale, SS yielded a significantly higher percent coverage compared to ALKB in the presence of modifications at 18 out of the 19 positions examined (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e5\u003c/span\u003e), corresponding well with the significant difference in PF mitigation efficacy between the two approaches.\u003c/p\u003e\u003cp\u003eSince the ability of modifications to trigger premature RT termination varies depending on their characteristics and location, as well as the RTase used, we aimed to provide a high-resolution overview of the impact of known tRNA modifications in \u003cem\u003eE. coli\u003c/em\u003e on the quantitative analysis with AQRNA-seq.\u0026nbsp;We found that the impact varies considerably across modifications and even across different positions carrying the same modification. This suggests intricate interactions between modifications and RT, which warrant further investigation. Importantly, despite the use of the EXOI and SS methods, we identified modifications exhibiting a high propensity for PF and thus reducing accurate quantification of the corresponding tRNAs, including ms\u003csup\u003e2\u003c/sup\u003ei\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003eA at position 37 and k\u003csup\u003e2\u003c/sup\u003eC at the wobble position. In \u003cem\u003eE. coli\u003c/em\u003e, while k\u003csup\u003e2\u003c/sup\u003eC34 is found specifically in tRNA-Ile2-CAT-1, where it converts both the codon recognition (AUG to AUA) and amino acid specificity (from methionine to isoleucine) of the tRNA\u003csup\u003e\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e\u003c/sup\u003e, ms\u003csup\u003e2\u003c/sup\u003ei\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003eA37 is commonly found in tRNAs decoding UNN codons where it prevents codon misreading by stabilizing codon-anticodon interactions\u003csup\u003e\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e\u003c/sup\u003e. As such, compared to k\u003csup\u003e2\u003c/sup\u003eC34, ms\u003csup\u003e2\u003c/sup\u003ei\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003eA37 may represent a greater challenge to quantification due to a larger number of tRNAs affected. Consistent with our findings, substantial PF due to ms\u003csup\u003e2\u003c/sup\u003ei\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003eA has been highlighted in other studies\u003csup\u003e\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e\u003c/sup\u003e and leveraged to achieve quantitative mapping of the modification at single-nucleotide resolution\u003csup\u003e\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e\u003c/sup\u003e. The original AQRNA-seq method and especially the present revision incorporating EXOI and SS are aimed at sensitive and accurate quantification of small RNAs and not intended for mapping the small RNA epitranscriptome. Indeed, NGS-based modification mapping by RT-induced falloff or mutation is a well-established technology\u003csup\u003e\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e,\u003cspan additionalcitationids=\"CR37 CR38 CR39\" citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eIn summary, we developed a highly sensitive and readily-automated NGS RNA-seq method for small RNA quantification based on the incorporation of EXOI to reduce PD and SS to increase full-length reads. By increasing sensitivity and eliminating manual intervention throughout the library preparation process, AQRNA-seq and other NGS-based RNA-seq methods can now be fully automated to increase throughput and reduce the required amount of RNA input, thus facilitating the study of the human RNome\u003csup\u003e\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e"},{"header":"METHODS","content":"\u003cp\u003e\u003cb\u003eThe AQRNA-seq data analytical pipeline v 2.0.\u003c/b\u003e The original data analytical pipeline\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e was substantially reengineered for improved accuracy, efficiency, and user experience, resulting in the development of AQRNA-seq data analytical pipeline v2.0 (\u003cb\u003eExtended Data\u003c/b\u003e Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e5\u003c/span\u003e). Simplified procedures involved in the pipeline are described in this section, and the complete pipeline is available at GitHub (LINK to be added) with a user manual including extensive annotations for individual steps. In the revised protocol, fastp v 0.23.4\u003csup\u003e42\u003c/sup\u003e is used to assess the quality of raw sequence reads, and Cutadapt v 4.2\u003csup\u003e43\u003c/sup\u003e is used to clip both the upstream and downstream adaptors from the quality-filtered reads. For each sample, the forward and reverse reads are assembled using PEAR v 0.9.10\u003csup\u003e44\u003c/sup\u003e based on the identification of an overlap with statistical significance, producing a unified FASTQ file containing consensus reads in the forward read direction, which is then converted to the FASTA format using a custom Bash script. Subsequently, unique reads along with their occurrences within each sample are extracted using fastxtoolkit v0.0.13 (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://hannonlab.cshl.edu/fastx_toolkit/\u003c/span\u003e\u003cspan address=\"http://hannonlab.cshl.edu/fastx_toolkit/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e), followed by combining the results across all samples to generate a FASTA file containing all unique reads present in at least one of the samples. The FASTA file is then processed and used as the input for the BLAST (Basic Local Alignment Search Tool) analysis to align the unique reads against a custom reference sequence library encompassing members of the target RNA species (e.g., tRNA) using blastn v 2.6.0\u003csup\u003e45\u003c/sup\u003e. Custom R scripts are used to filter raw blastn results based on use-defined thresholds and to resolve the mapping of a single query sequence to multiple reference sequences based on e-value, bit-score, percent identity, and query coverage. Finally, the processed blastn results are cross-referenced with the occurrences of unique reads in samples to generate a raw tRNA abundance matrix, using a custom R script, wherein the read counts are shown for each subject in the reference sequence library within each sample. Detailed instructions for constructing reference sequence libraries can be found in the section \u0026ldquo;Construction of the \u003cem\u003eE. coli\u003c/em\u003e BW25113 tRNA reference sequence library\u0026rdquo;.\u003c/p\u003e\u003cp\u003e\u003cb\u003eBacterial strains, culturing conditions, and small RNA isolation.\u003c/b\u003e \u003cem\u003eE. coli\u003c/em\u003e strain BW25113 was selected for use in this study as it has been extensively studied as a model strain and represents the parent strain of the Keio Collection\u003csup\u003e\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e\u003c/sup\u003e, thereby facilitating the comparison of results derived from this study and potential subsequent studies. The stock culture was preserved in lysogeny broth (LB; Fisher BioReagent) with 10% DMSO at -80\u0026deg;C. Frozen stock culture was streaked on a LB agar plate, followed by incubation at 37\u0026deg;C for 16\u0026ndash;18 h. Single colonies from the freshly streaked plate were inoculated into 5 mL LB broth, followed by overnight incubation on a spinning wheel (New Brunswick Scientific) at 37\u0026deg;C. Bacterial cultures were then sub-cultured (1:100) in 25 mL LB broth, followed by incubation in a shaking incubator (Shon\u0026rsquo;s Scientific Refrigeration) at 37\u0026deg;C with 250 rpm. Upon entry into the late-exponential phase, as indicated by optical density (OD) 600 values between 0.6\u0026ndash;0.7, the bacterial culture was aliquoted (1 mL) into sterile Eppendorf tubes. Bacterial cell pellets were harvested by centrifugation at 500 \u0026times; g and 4\u0026deg;C for 5 min and subsequently washed by resuspension in 1 mL of 1X phosphate buffered saline (PBS; pH\u0026thinsp;=\u0026thinsp;7.4; Thermo Fisher Scientific), followed by centrifugation using the same settings.\u003c/p\u003e\u003cp\u003eIsolation of small RNAs from the bacterial cells was conducted following the user guide of the PureLink miRNA Isolation Kit (Thermo Fisher Scientific) with modifications. Briefly, bacterial cells were lysed by adding 1 mL of the TRIzol Reagent (Sigma-Aldrich) per 1 \u0026times; 10\u003csup\u003e7\u003c/sup\u003e cells (equivalent to approximately 1 mL of bacterial culture in late exponential phase), followed by incubation at room temperature for 5 min. To isolate total RNAs, the cell lysate was incubated with 0.2 mL of chloroform (Macron Fine Chemicals) at room temperature for 3 min, followed by centrifugation at 12,000 \u0026times; g and 4\u0026deg;C for 15 min. The upper aqueous phase containing total RNAs was then transferred into a sterile RNase-free tube. Subsequently, small RNAs were separated from large RNAs by passing the total RNAs through two Spin Cartridges to selectively bind large and small RNAs to the membrane, with 35% and 70% ethanol, respectively. Following two wash steps, small RNAs were eluted in RNase-free water by centrifugation at 16,000 \u0026times; g for 1 min.\u003c/p\u003e\u003cp\u003e\u003cb\u003eLibrary preparation.\u003c/b\u003e cDNA libraries were prepared using 75 ng of \u003cem\u003eE. coli\u003c/em\u003e small RNAs as the input, following four variations of the AQRNA-seq library preparation protocol, each with three biological replicates, as illustrated in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e. Protocol variant 1 (V1) represents the original protocol initially published in Hu et al.\u003csup\u003e14\u003c/sup\u003e and later elaborated on in Chen et al.\u003csup\u003e13\u003c/sup\u003e, while protocol variants 2 (V2), 3 (V3), and 4 (V4) involve alternative adaptations of the original protocol with respect to the approaches for mitigating polymerase fall-off (PF) and/or primer dimers (PD). In the following sections, protocol V1 will be briefly recapitulated, highlighting the molecular processes involved in individual steps, and the alternative protocols (V2, V3, and V4) will be introduced, emphasizing the specific adaptations made in comparison to the original protocol. Full protocols with step-by-step instructions are provided in \u003cb\u003eSupplementary Notes 1\u0026ndash;4\u003c/b\u003e to ensure reproducibility and comprehensive evaluation of the study. Detailed information regarding the DNA and RNA oligonucleotides used in the protocols can be found in \u003cb\u003eSupplementary Table\u0026nbsp;10\u003c/b\u003e. Detailed information regarding the reagents, chemicals, and kits/columns can be found in \u003cb\u003eSupplementary Table\u0026nbsp;11\u003c/b\u003e.\u003c/p\u003e\u003cp\u003e\u003cem\u003eProtocol V1.\u003c/em\u003e Small RNAs mixed with an internal standard were dephosphorylated using Shrimp Alkaline Phosphatase (New England Biolabs) through incubation at 37\u0026deg;C for 30 min and then at 65\u0026deg;C for 5 min to both deactivate the enzyme and denature the RNAs. T4 RNA Ligase 1 (New England Biolabs) was used to ligate Linker 1 to the 3\u0026rsquo; end of RNAs through incubation at 25\u0026ordm;C for 2 h and then at 16\u0026deg;C for 16 h. To mitigate PF, the AlkB demethylase (ArrayStar) was used to remove post-transcriptional methylations potentially present on RNAs through incubation at 25\u0026deg;C for 2 h (i.e., the AlkB approach). Excess Linker 1 not ligated to RNAs was then deadenylated using 5\u0026rsquo;-deadenylase (New England Biolabs) through incubation at 30\u0026deg;C for 1 h and subsequently digested with RecJf exonuclease (New England Biolabs) through two 30-min incubations at 37\u0026deg;C. Reverse transcription (RT) was initiated with a 2-min incubation at 80\u0026deg;C to anneal the RT primers to the RNA template, followed by a 2-h incubation at 50\u0026deg;C to synthesize the cDNA strand using the PrimeScript\u0026trade; reverse transcriptase (RTase; TaKaRa) and a 15-min incubation at 70\u0026deg;C to deactivate the enzyme. After the synthesis of the cDNA strand, the RNA template was hydrolyzed using sodium hydroxide (Sigma-Aldrich) through incubation at 95\u0026deg;C for 3 min, followed by immediate neutralization using hydrochloric acid (VWR). T4 DNA ligase (New England Biolabs) was used to ligate Linker 2 to the 3\u0026rsquo; end of cDNAs, and the excess Linker 2 were subsequently removed following the same procedures as described for Linker 1. Sequencing adaptors were then attached to both ends of target cDNA fragments through polymerase chain reaction (PCR), and the PCR products were subsequently subject to 3% agarose gel electrophoresis to separate target fragments from unwanted PD. To mitigate the carry-over of PD into sequencing libraries, gel blocks containing target fragments were excised, followed by gel extraction to retrieve and purify the target fragments (i.e., the GE approach,. Figure\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e).\u003c/p\u003e\u003cp\u003e\u003cem\u003eProtocol V2.\u003c/em\u003e This variation of the protocol incorporates a novel approach (SS) to mitigate PF during the RT reaction, a predominant cause of premature RT termination (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). Specifically, the SuperScript IV RTase (Thermo Fisher Scientific) was used in substitution of the PrimeScript RTase (TaKaRa) for synthesizing the cDNA strand. SuperScript IV has demonstrated superior processivity (i.e., the ability to consecutively add nucleotides without releasing the RNA strand) as compared to other RTases\u003csup\u003e\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e\u003c/sup\u003e, indicating its potential to minimize PF induced by a broader range of post-transcriptional modifications, not just methylations, thereby facilitating the synthesis of full-length cDNAs. As a result, the AlkB demethylation step was eliminated from the protocol.\u003c/p\u003e\u003cp\u003e\u003cem\u003eProtocol V3.\u003c/em\u003e This variation of the protocol incorporates a novel approach (EXOI) to mitigate the carry-over of PD into the constructed cDNA libraries (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). Specifically, immediately following the RT reaction, exonuclease I from \u003cem\u003eE. coli\u003c/em\u003e (New England Biolabs) was used to digest excess RT primers in the 3\u0026rsquo; to 5\u0026rsquo; direction, thereby inhibiting the formation of PD (i.e., ligation of DNA linkers to the 3\u0026rsquo; end of RT primers) from the first place. As a result, the gel extraction step was eliminated from the protocol.\u003c/p\u003e\u003cp\u003e\u003cem\u003eProtocol V4.\u003c/em\u003e This variation of the protocol incorporates both the novel approach to mitigate PF (SS), as previously described for protocol V2, and the novel approach to mitigate PD (EXOI), as previously described for protocol V3 (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). Consequently, both the AlkB demethylation and the gel extraction steps were obviated by using this protocol.\u003c/p\u003e\u003cp\u003e\u003cb\u003eLibrary quality assessment and sequencing.\u003c/b\u003e The constructed cDNA libraries were submitted to the MIT BioMicro Center, along with custom sequencing primers (\u003cb\u003eSupplementary Table\u0026nbsp;10\u003c/b\u003e). The quality of the libraries was assessed using AATI Fragment Analysis (Agilent) and LightCycler 480 Real-Time PCR System (Roche), and 75-bp Paired-End sequencing was performed on an Illumina MiSeq platform with the v3 reagent kit (Illumina).\u003c/p\u003e\u003cdiv id=\"Sec5\" class=\"Section2\"\u003e\u003ch2\u003eData analysis.\u003c/h2\u003e\u003cp\u003eC\u003cem\u003eonstruction of the E. coli BW25113 tRNA reference sequence library\u003c/em\u003e. To maximize the accuracy for quantifying the landscape of tRNAs, we followed a custom workflow developed in our laboratory to construct the reference sequence library for use in the data analytical pipeline. Specifically, the mature tRNA sequences of \u003cem\u003eE. coli\u003c/em\u003e BW25113 were obtained from GtRNAdb\u003csup\u003e\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e,\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e\u003c/sup\u003e, followed by the removal of predicted pseudogenes and tRNAs with unknown or undetermined isotypes. \u0026ldquo;CCA\u0026rdquo; was manually appended to sequences lacking a 3\u0026rsquo; overhang ending with \u0026ldquo;CCA\u0026rdquo;, based on the structural alignments generated using domain-specific covariance models in COVE\u003csup\u003e\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e\u003c/sup\u003e. Subsequently, all \u0026ldquo;U\u0026rdquo; in the sequences were replaced with \u0026ldquo;T\u0026rdquo;, and identical sequences representing multiple gene copies encoding the same isodecoder were deduplicated, preserving only one sequence in the library as a representative. The final library comprised 50 unique sequences representing distinct tRNAs at the isodecoder level, alongside three control sequences for benchmarking purposes.\u003c/p\u003e\u003cp\u003e\u003cem\u003eTransfer RNA (tRNA) abundance estimation and normalization.\u003c/em\u003e The raw sequence files in FASTQ format were processed through the AQRNA-seq data analytical pipeline v 2.0 (\u003cb\u003eExtended Data\u003c/b\u003e Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e; detailed in the section \u0026ldquo;Optimization of AQRNA-seq data analytical pipeline\u0026rdquo;) to derive the raw tRNA abundance matrix. To enable downstream comparisons across library preparation protocols while accounting for both sequencing depth and sample composition, the raw tRNA abundance was normalized using the median of ratios method, as implemented in the DESeq2 v 1.36.0\u003csup\u003e49\u003c/sup\u003e package in R.\u003c/p\u003e\u003cp\u003e\u003cem\u003eEfficacy assessment and comparison of PD mitigation approaches.\u003c/em\u003e The presence of PD in the constructed sequencing libraries leads to sequence reads lacking inserts or containing short inserts. Therefore, for each sample, we calculated the proportion of quality-filtered reads with different insert lengths by running Cutadapt v 4.2\u003csup\u003e43\u003c/sup\u003e iteratively, setting various thresholds for the minimum insert length. The efficacy of PD mitigation approaches was subsequently inferred from the proportion of reads with inserts longer than 10 bp. Specifically, the read proportion data was fitted with a linear regression model, using the stats v 4.3.2\u003csup\u003e50\u003c/sup\u003e package in R, wherein PD mitigation approach (GE vs ExoI) was included as the primary variable of interest and biological replicate as a blocking effect to account for potential experimental variations. To account for potential effects of PF mitigation approaches on read proportion or the efficacy of PD mitigation approaches, additive and interactive models including PF mitigation approach (AlkB vs SS) as an additional factor were also constructed, and the final model was determined based on F tests.\u003c/p\u003e\u003cp\u003e\u003cem\u003eEfficacy assessment and comparison of PF mitigation approaches.\u003c/em\u003e As PF eventually results in partial alignments of sequence reads with the reference sequences, the proportion of full-length reads could be used as a performance measure for assessing the efficacy of different PF mitigation approaches. Therefore, we calculated the proportion of full-length reads for individual tRNAs, with full-length reads defined as those aligned with the reference sequence from the 3\u0026rsquo; to the 5\u0026rsquo; ends, allowing for up to five missing bases at both ends. Subsequently, the full-length read proportion data was fitted with a linear regression model, using the stats v 4.3.2\u003csup\u003e50\u003c/sup\u003e package in R, wherein PF mitigation approach was included as the primary variable of interest and biological replicate as a blocking effect. Additive and interactive models including tRNA and PD mitigation approach as additional factors were also constructed, and the final model was determined based on F tests. In cases where interactions in model were deemed significant, post hoc comparisons were performed using the emmeans v 1.10.0\u003csup\u003e51\u003c/sup\u003e package in R, with the Benjamini-Hochberg (BH) method for multiple corrections.\u003c/p\u003e\u003cp\u003e\u003cem\u003eQuantitative comparison of tRNA landscapes across library preparation protocols.\u003c/em\u003e The concordance in the tRNA landscape obtained from different protocols was assessed by calculating the pairwise Pearson\u0026rsquo;s correlation coefficients (r) among the protocols based on the normalized abundance of individual tRNAs, averaged across biological replicates, using the PerformanceAnalytics v 2.0.4\u003csup\u003e52\u003c/sup\u003e package in R. To objectively rank the performance of the protocols to accurately depict the cellular tRNA landscape, we compared the normalized tRNA abundances from each protocol with an orthogonal tRNA abundance dataset\u003csup\u003e\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e\u003c/sup\u003e, which was derived using two-dimensional gel electrophoresis to separate tRNAs, followed by quantification using quantitative Northern blot. The quantitative performance of each protocol was inferred from the Pearson\u0026rsquo;s r between its quantitative results and the orthogonal dataset, with higher correlations denoting superior performance.\u003c/p\u003e\u003cp\u003e\u003cem\u003eAssessment and comparison of the effect of PF and PD mitigation approaches on the estimated abundance of tRNAs.\u003c/em\u003e To investigate whether variations in quantitative results across protocols could be attributed to the incorporation of different PF and PD mitigation approaches, the raw tRNA abundance data was transformed using the regularized logarithm (rlog) method implemented in the DESeq2 v 1.36.0\u003csup\u003e49\u003c/sup\u003e package in R, and hierarchical clustering was performed based on the Euclidean distance matrix derived from the rlog-transformed, normalized abundance data using the pheatmap v 1.0.12\u003csup\u003e53\u003c/sup\u003e package in R. Furthermore, principal component analysis (PCA) was performed on the same data using DESeq2 to discern whether any effects related to the use of different PF and PD mitigation approaches would result in distinct clusters among samples subject to different protocols.\u003c/p\u003e\u003cp\u003eTo further scrutinize the effect of PF and PD mitigation approaches on estimated abundance at the level of individual tRNAs, we followed the differential expression analysis workflow implemented in DESeq2 to identify tRNAs exhibiting differential estimated abundance as influenced by PF and/or PD mitigation approaches. Briefly, the workflow includes the following steps: (i) normalization of the raw abundance data using the median of ratios method; (ii) dispersion parameter estimation for individual tRNAs using empirical Bayes shrinkage; (iii) construction of a negative binomial regression model for each tRNA; (iv) estimation of logarithmic (log) fold changes using Approximate Posterior Estimation for generalized linear model\u003csup\u003e\u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e54\u003c/span\u003e\u003c/sup\u003e; and (v) Wald tests for significant testing (with the BH method for multiple correction). In each negative binomial regression model, the normalized tRNA abundance was designated as the response variable and biological replicate was included as a blocking effect; both PF and PD mitigation approaches were included as primary factors of interest, along with the interaction term between these two factors.\u003c/p\u003e\u003cp\u003e\u003cem\u003ePosition-specific percent coverage and identity calculation.\u003c/em\u003e To obtain a granular view on the alignment of reads with the reference sequences and to facilitate the mechanistic understanding of variations in quantitative results across protocols, we calculated position-specific percent coverage and identity for each reference sequence in each sample. Specifically, two FASTA files were generated for each reference sequence in each sample, containing the reference sequence and all its mapped reads, respectively. blastn v 2.6.0\u003csup\u003e45\u003c/sup\u003e was used to align the reads with the reference sequence, specifying SAM as the output format, followed by conversion of the output file into the BAM format using samtools view v 1.19.2\u003csup\u003e55\u003c/sup\u003e. The coverage at each position within a reference sequence across all aligned reads was extracted using samtools depth v 1.19.2\u003csup\u003e55\u003c/sup\u003e and subsequently divided by the maximum coverage within the reference sequence to derive the position-specific percent coverage. samtools mpileup v 1.19.2\u003csup\u003e55\u003c/sup\u003e was then used to extract the pileup results of the read alignments for each reference sequence, which were parsed to obtain the position-specific percent identity using a custom R script. For each position, this was calculated as the percentage of reads that contain the same base as the reference sequence.\u003c/p\u003e\u003cp\u003e\u003cem\u003eMapping of tRNA modifications onto the reference sequences.\u003c/em\u003e The RF00005 alignments for all \u003cem\u003eE. coli\u003c/em\u003e tRNAs available on Modomics\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e\u003c/sup\u003e were downloaded. blastn v 2.6.0\u003csup\u003e40\u003c/sup\u003e was used to identify matches (100% query coverage; \u0026gt;98% percent identity) between our reference sequences used for mapping reads and the tRNAs present in the alignments. Reference sequences without a match with any tRNAs in the alignments were excluded from further analyses. For mapping tRNA modifications onto reference sequences, positions of tRNA modifications present in the alignments were exported from Modomics, specifying short name for the symbol filter to represent modification identities. These modifications were subsequently assigned to specific positions on the reference sequences by cross-referencing the original positions within the reference sequences and their corresponding positions in the alignments.\u003c/p\u003e\u003cp\u003e\u003cem\u003eAssessment of the impact of tRNA modifications on RT.\u003c/em\u003e As the interference of tRNA modifications with the RT process represents a predominant cause contributing to premature RT terminations, leading to suboptimal quantitative performance, we aimed to infer the adverse impact of these modifications on RT from the position-specific percent coverage data, assuming more adverse impacts would correlate with lower percent coverages. Due to the interdependence of percent coverages at a specific position with those closer to the 3\u0026rsquo; end, it was difficult to assess the absolute impact of modifications at individual positions. Nevertheless, our data enabled a direct comparison of the impact of modifications at all possible positions across different protocols. To this end, the position-specific percent coverage data in the presence of tRNA modifications were obtained for each position within the alignments, and positions with insufficient data points (n\u0026thinsp;\u0026lt;\u0026thinsp;3) for any of the protocols were precluded from further analyses. For each retained position, the data were fitted with a linear regression model, using the stats v 4.2.1\u003csup\u003e50\u003c/sup\u003e package in R, designating percent coverage as the response variable, protocol as the primary factor of interest, and biological replicate as a blocking effect. In cases where percent coverage at a specific position exhibited significant variations across different protocols, post hoc pairwise comparisons (i.e., Tukey\u0026rsquo;s HSD tests) were performed among protocols using the emmeans v 1.8.1.1\u003csup\u003e51\u003c/sup\u003e package in R.\u003c/p\u003e\u003cp\u003eWe next calculated the percent coverage resulting from all unique combinations of modification and position, subjected to different protocols, aiming to identify combinations imposing greatest challenges in RT and thus the accurate quantification of AQRNA-seq.\u0026nbsp;For each protocol, we normalized the percent coverages to the maximum value at the corresponding position in the presence of modifications and subsequently classified the normalized percent coverages into ultra-high (75\u0026ndash;100%), high (50\u0026ndash;75%), low (25\u0026ndash;50%), and ultra-low (0\u0026ndash;25%) categories. For modification-position combinations resulting in ultra-low normalized percent coverages, the read alignments against the reference sequence of the corresponding tRNAs were manually examined to confirm their propensities to induce premature RT terminations.\u003c/p\u003e\u003cp\u003e\u003cem\u003eAssessment of the ability of AQRNA-seq in quantitatively mapping tRNA modifications.\u003c/em\u003e While quantitative mapping of tRNA modifications through NGS-based methods may depend on read pileups resulting from RT stops at modified nucleotides, our endeavors towards enhancing quantitative accuracy involve minimizing premature RT terminations, inherently deprecating reliance on this approach. However, certain tRNA modifications may prompt the insertion of incorrect residues as the RTase traverses them, resulting in unique mutation signatures (e.g., reduced percent identity and specific base profiles) that could be identified through subsequent analysis of sequencing data. Hence, to elucidate the potential of different protocols in the quantitative mapping of tRNA modifications, we calculated the average percent identity in the presence of each modification at every possible position, subjected to each of the protocols.\u003c/p\u003e\u003c/div\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgements\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors thank the MIT BioMicro Center and its Director, Dr. Stuart Levine, for support and advice during the performance of the studies presented here. The authors gratefully acknowledge funding from the Singapore National Science Foundation under the Singapore-MIT Alliance for Research and Technology Antimicrobial Resistance Interdisciplinary Research Group (PCD), National Institutes of Health Transformative Award ES031576 (PCD), and Center grant P30-ES002109 from the National Institute of Environmental Health Sciences of the National Institutes of Health.\u003c/p\u003e\n\n\u003cp\u003e\u003cstrong\u003eData availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eSequencing data are being deposited in the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) and will be available upon request prior to publication.\u003c/p\u003e\n\n\u003cp\u003e\u003cstrong\u003eCode availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe complete data analysis pipeline is available at GitHub (https://github.com/dedonlab/AQRNA-seq-method-optimization.git)\u003c/p\u003e\n\n\u003cp\u003e\u003cstrong\u003eAuthor contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eR.C., L.L., B.C., and P.C.D. designed the experiments. R.C., M.S.D., and P.C.D. wrote the first draft of the manuscript. R.C., D.Y., and L.L. performed experiments to develop and apply the sequencing method. R.C. developed the data processing pipeline and performed all statistical analyses. B.C. and P.C.D. oversaw study design and analysis. All authors contributed to the manuscript writing.\u003c/p\u003e\n"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eWang, Z., Gerstein, M. \u0026amp; Snyder, M. RNA-Seq: a revolutionary tool for transcriptomics. \u003cem\u003eNat Rev Genet\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 57-63, (2009).\u003c/li\u003e\n\u003cli\u003eByron, S. A., Van Keuren-Jensen, K. R., Engelthaler, D. M., Carpten, J. D. \u0026amp; Craig, D. W. Translating RNA sequencing into clinical diagnostics: opportunities and challenges. \u003cem\u003eNat Rev Genet\u003c/em\u003e \u003cstrong\u003e17\u003c/strong\u003e, 257-271, (2016).\u003c/li\u003e\n\u003cli\u003eBenesova, S., Kubista, M. \u0026amp; Valihrach, L. Small RNA-Sequencing: Approaches and Considerations for miRNA Analysis. \u003cem\u003eDiagnostics (Basel)\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 964, (2021).\u003c/li\u003e\n\u003cli\u003eMotorin, Y. \u0026amp; Helm, M. Methods for RNA Modification Mapping Using Deep Sequencing: Established and New Emerging Technologies. \u003cem\u003eGenes (Basel)\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 35 (2019).\u003c/li\u003e\n\u003cli\u003eZhang, Y., Lu, L. \u0026amp; Li, X. Detection technologies for RNA modifications. \u003cem\u003eExp Mol Med\u003c/em\u003e \u003cstrong\u003e54\u003c/strong\u003e, 1601-1616, (2022).\u003c/li\u003e\n\u003cli\u003eZhang, L. S., Dai, Q. \u0026amp; He, C. Base-Resolution Sequencing Methods for Whole-Transcriptome Quantification of mRNA Modifications. \u003cem\u003eAcc Chem Res\u003c/em\u003e \u003cstrong\u003e57\u003c/strong\u003e, 47-58, (2024).\u003c/li\u003e\n\u003cli\u003eCummings, B. B.\u003cem\u003e et al.\u003c/em\u003e Improving genetic diagnosis in Mendelian disease with transcriptome sequencing. \u003cem\u003eSci Transl Med\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, eaal5209, (2017).\u003c/li\u003e\n\u003cli\u003eHaque, A., Engel, J., Teichmann, S. A. \u0026amp; Lonnberg, T. A practical guide to single-cell RNA-sequencing for biomedical research and clinical applications. \u003cem\u003eGenome Med\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 75, (2017).\u003c/li\u003e\n\u003cli\u003eAlon, S.\u003cem\u003e et al.\u003c/em\u003e Barcoding bias in high-throughput multiplex sequencing of miRNA. \u003cem\u003eGenome Res\u003c/em\u003e \u003cstrong\u003e21\u003c/strong\u003e, 1506-1511, (2011).\u003c/li\u003e\n\u003cli\u003eFuchs, R. T., Sun, Z., Zhuang, F. \u0026amp; Robb, G. B. Bias in ligation-based small RNA sequencing library construction is determined by adaptor and RNA structure. \u003cem\u003ePLoS One\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, e0126049, (2015).\u003c/li\u003e\n\u003cli\u003ePang, Y. L., Abo, R., Levine, S. S. \u0026amp; Dedon, P. C. Diverse cell stresses induce unique patterns of tRNA up- and down-regulation: tRNA-seq for quantifying changes in tRNA copy number. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e42\u003c/strong\u003e, e170, (2014).\u003c/li\u003e\n\u003cli\u003eLi, F.\u003cem\u003e et al.\u003c/em\u003e Regulatory impact of RNA secondary structure across the Arabidopsis transcriptome. \u003cem\u003ePlant Cell\u003c/em\u003e \u003cstrong\u003e24\u003c/strong\u003e, 4346-4359, (2012).\u003c/li\u003e\n\u003cli\u003eChen, R., Yim, D. \u0026amp; Dedon, P. C. AQRNA-seq for Quantifying Small RNAs. \u003cem\u003eJ Vis Exp\u003c/em\u003e, 2024 Feb 2;(204), (2024).\u003c/li\u003e\n\u003cli\u003eHu, J. F.\u003cem\u003e et al.\u003c/em\u003e Quantitative mapping of the cellular small RNA landscape with AQRNA-seq. \u003cem\u003eNat Biotechnol\u003c/em\u003e \u003cstrong\u003e39\u003c/strong\u003e, 978-988, (2021).\u003c/li\u003e\n\u003cli\u003eDai, Q., Zheng, G., Schwartz, M. H., Clark, W. C. \u0026amp; Pan, T. Selective Enzymatic Demethylation of N(2) ,N(2) -Dimethylguanosine in RNA and Its Application in High-Throughput tRNA Sequencing. \u003cem\u003eAngew Chem Int Ed Engl\u003c/em\u003e \u003cstrong\u003e56\u003c/strong\u003e, 5017-5020, (2017).\u003c/li\u003e\n\u003cli\u003eDong, H., Nilsson, L. \u0026amp; Kurland, C. G. Co-variation of tRNA abundance and codon usage in Escherichia coli at different growth rates. \u003cem\u003eJ Mol Biol\u003c/em\u003e \u003cstrong\u003e260\u003c/strong\u003e, 649-663, (1996).\u003c/li\u003e\n\u003cli\u003eRychlik, W. Selection of primers for polymerase chain reaction. \u003cem\u003eMol Biotechnol\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 129-134, (1995).\u003c/li\u003e\n\u003cli\u003eBrownie, J.\u003cem\u003e et al.\u003c/em\u003e The elimination of primer-dimer accumulation in PCR. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e25\u003c/strong\u003e, 3235-3241, (1997).\u003c/li\u003e\n\u003cli\u003eIllumina, Inc. Adapter dimers causes, effects, and how to remove them. https://knowledge.illumina.com/library-preparation/general/library-preparation-general-troubleshooting-list/000001911, (2023).\u003c/li\u003e\n\u003cli\u003eNew England Biolabs, Inc. Six Tips for a Perfect Gel Extraction. https://www.neb.com/en-us/tools-and-resources/usage-guidelines/six-tips-for-a-perfect-gel-extraction#, (2024).\u003c/li\u003e\n\u003cli\u003eKeijzers, G., Liu, D. \u0026amp; Rasmussen, L. J. Exonuclease 1 and its versatile roles in DNA repair. \u003cem\u003eCrit Rev Biochem Mol Biol\u003c/em\u003e \u003cstrong\u003e51\u003c/strong\u003e, 440-451, (2016).\u003c/li\u003e\n\u003cli\u003eCappannini, A.\u003cem\u003e et al.\u003c/em\u003e MODOMICS: a database of RNA modifications and related information. 2023 update. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e52\u003c/strong\u003e, D239-D244, (2024).\u003c/li\u003e\n\u003cli\u003eNordmann, P. L., Makris, J. C. \u0026amp; Reznikoff, W. S. Inosine induced mutations. \u003cem\u003eMol Gen Genet\u003c/em\u003e \u003cstrong\u003e214\u003c/strong\u003e, 62-67, (1988).\u003c/li\u003e\n\u003cli\u003eWeber, M.\u003cem\u003e et al.\u003c/em\u003e Chromosome-wide and promoter-specific analyses identify sites of differential DNA methylation in normal and transformed human cells. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e37\u003c/strong\u003e, 853-862, (2005).\u003c/li\u003e\n\u003cli\u003eSchmidt, D.\u003cem\u003e et al.\u003c/em\u003e ChIP-seq: using high-throughput sequencing to discover protein-DNA interactions. \u003cem\u003eMethods\u003c/em\u003e \u003cstrong\u003e48\u003c/strong\u003e, 240-248, (2009).\u003c/li\u003e\n\u003cli\u003eMeers, M. P., Bryson, T. D., Henikoff, J. G. \u0026amp; Henikoff, S. Improved CUT\u0026amp;RUN chromatin profiling tools. \u003cem\u003eElife\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, e46314, (2019).\u003c/li\u003e\n\u003cli\u003eBronner, I. F., Quail, M. A., Turner, D. J. \u0026amp; Swerdlow, H. Improved Protocols for Illumina Sequencing. \u003cem\u003eCurr Protoc Hum Genet\u003c/em\u003e \u003cstrong\u003e80\u003c/strong\u003e, 18 12 11-42, (2014).\u003c/li\u003e\n\u003cli\u003eGao, X.\u003cem\u003e et al.\u003c/em\u003e A reassessment of several erstwhile methods for isolating DNA fragments from agarose gels. \u003cem\u003e3 Biotech\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 138, (2021).\u003c/li\u003e\n\u003cli\u003eShore, S., Henderson, J. M. \u0026amp; McCaffrey, A. P. CleanTag Adapters Improve Small RNA Next-Generation Sequencing Library Preparation by Reducing Adapter Dimers. \u003cem\u003eMethods Mol Biol\u003c/em\u003e \u003cstrong\u003e1712\u003c/strong\u003e, 145-161, (2018).\u003c/li\u003e\n\u003cli\u003eChen, Z.\u003cem\u003e et al.\u003c/em\u003e Transfer RNA demethylase ALKBH3 promotes cancer progression via induction of tRNA-derived small RNAs. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e47\u003c/strong\u003e, 2533-2545, (2019).\u003c/li\u003e\n\u003cli\u003eFalnes, P. O. Repair of 3-methylthymine and 1-methylguanine lesions by bacterial and human AlkB proteins. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e32\u003c/strong\u003e, 6260-6267, (2004).\u003c/li\u003e\n\u003cli\u003eCozen, A. E.\u003cem\u003e et al.\u003c/em\u003e ARM-seq: AlkB-facilitated RNA methylation sequencing reveals a complex landscape of modified tRNA fragments. \u003cem\u003eNat Methods\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 879-884, (2015).\u003c/li\u003e\n\u003cli\u003eZucha, D., Androvic, P., Kubista, M. \u0026amp; Valihrach, L. Performance Comparison of Reverse Transcriptases for Single-Cell Studies. \u003cem\u003eClin Chem\u003c/em\u003e \u003cstrong\u003e66\u003c/strong\u003e, 217-228, (2020).\u003c/li\u003e\n\u003cli\u003eMuramatsu, T.\u003cem\u003e et al.\u003c/em\u003e Codon and amino-acid specificities of a transfer RNA are both converted by a single post-transcriptional modification. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e336\u003c/strong\u003e, 179-181, (1988).\u003c/li\u003e\n\u003cli\u003eWilson, R. K. \u0026amp; Roe, B. A. Presence of the hypermodified nucleotide N6-(delta 2-isopentenyl)-2-methylthioadenosine prevents codon misreading by Escherichia coli phenylalanyl-transfer RNA. \u003cem\u003eProc Natl Acad Sci U S A\u003c/em\u003e \u003cstrong\u003e86\u003c/strong\u003e, 409-413, (1989).\u003c/li\u003e\n\u003cli\u003eBehrens, A., Rodschinka, G. \u0026amp; Nedialkova, D. D. High-resolution quantitative profiling of tRNA abundance and modification status in eukaryotes by mim-tRNAseq. \u003cem\u003eMol Cell\u003c/em\u003e \u003cstrong\u003e81\u003c/strong\u003e, 1802-1815, (2021).\u003c/li\u003e\n\u003cli\u003eFang, Z.\u003cem\u003e et al.\u003c/em\u003e The Transcriptome-Wide Mapping of 2-Methylthio-N(6)-isopentenyladenosine at Single-Base Resolution. \u003cem\u003eJ Am Chem Soc\u003c/em\u003e \u003cstrong\u003e145\u003c/strong\u003e, 5467-5473, (2023).\u003c/li\u003e\n\u003cli\u003eBao, Z., Li, T. \u0026amp; Liu, J. Determining RNA Natural Modifications and Nucleoside Analog-Labeled Sites by a Chemical/Enzyme-Induced Base Mutation Principle. \u003cem\u003eMolecules\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 1517, (2023).\u003c/li\u003e\n\u003cli\u003eDebnath, T. K. \u0026amp; Xhemalce, B. Deciphering RNA modifications at base resolution: from chemistry to biology. \u003cem\u003eBrief Funct Genomics\u003c/em\u003e \u003cstrong\u003e20\u003c/strong\u003e, 77-85, (2021).\u003c/li\u003e\n\u003cli\u003eKietrys, A. M., Velema, W. A. \u0026amp; Kool, E. T. Fingerprints of Modified RNA Bases from Deep Sequencing Profiles. \u003cem\u003eJ Am Chem Soc\u003c/em\u003e \u003cstrong\u003e139\u003c/strong\u003e, 17074-17081, (2017).\u003c/li\u003e\n\u003cli\u003eNational Academies of Sciences, Engineering, and Medicine. \u003cem\u003eCharting a Future for Sequencing RNA and Its Modifications: A New Era for Biology and Medicine\u003c/em\u003e. https://www.ncbi.nlm.nih.gov/pubmed/39159274, (2024).\u003c/li\u003e\n\u003cli\u003eChen, S., Zhou, Y., Chen, Y. \u0026amp; Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e34\u003c/strong\u003e, i884-i890, (2018).\u003c/li\u003e\n\u003cli\u003eMartin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. \u003cem\u003eEMBnet J\u003c/em\u003e \u003cstrong\u003e17\u003c/strong\u003e, 3, (2011).\u003c/li\u003e\n\u003cli\u003eZhang, J., Kobert, K., Flouri, T. \u0026amp; Stamatakis, A. PEAR: a fast and accurate Illumina Paired-End reAd mergeR. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e30\u003c/strong\u003e, 614-620, (2014).\u003c/li\u003e\n\u003cli\u003eAltschul, S. F., Gish, W., Miller, W., Myers, E. W. \u0026amp; Lipman, D. J. Basic local alignment search tool. \u003cem\u003eJ Mol Biol\u003c/em\u003e \u003cstrong\u003e215\u003c/strong\u003e, 403-410, (1990).\u003c/li\u003e\n\u003cli\u003eBaba, T.\u003cem\u003e et al.\u003c/em\u003e Construction of Escherichia coli K-12 in-frame, single-gene knockout mutants: the Keio collection. \u003cem\u003eMol Syst Biol\u003c/em\u003e \u003cstrong\u003e2\u003c/strong\u003e, 2006 0008, (2006).\u003c/li\u003e\n\u003cli\u003eChan, P. P. \u0026amp; Lowe, T. M. GtRNAdb: a database of transfer RNA genes detected in genomic sequence. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e37\u003c/strong\u003e, D93-97, (2009).\u003c/li\u003e\n\u003cli\u003eChan, P. P. \u0026amp; Lowe, T. M. GtRNAdb 2.0: an expanded database of transfer RNA genes identified in complete and draft genomes. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e44\u003c/strong\u003e, D184-189, (2016).\u003c/li\u003e\n\u003cli\u003eLove, M. I., Huber, W. \u0026amp; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. \u003cem\u003eGenome Biol\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 550, (2014).\u003c/li\u003e\n\u003cli\u003eR Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Available online at https://www.R-project.org/, (2018).\u003c/li\u003e\n\u003cli\u003eLenth, R. \u003cem\u003eemmeans: estimated marginal means, aka least-squares means. R package version 1.10.0\u003c/em\u003e, https://CRAN.R-project.org/package=emmeans, (2024).\u003c/li\u003e\n\u003cli\u003ePeterson, B. \u0026amp; Carl, P. \u003cem\u003ePerformanceAnalytics: econometric tools for performance and risk analysis. R package version 2.0.4.\u003c/em\u003e, https://CRAN.R-project.org/package=PerformanceAnalytics, (2020).\u003c/li\u003e\n\u003cli\u003eKolde, R. \u003cem\u003epheatmap: pretty heatmaps. R package version 1.0.12.\u003c/em\u003e, https://CRAN.R-project.org/package=pheatmap, (2019).\u003c/li\u003e\n\u003cli\u003eZhu, A., Ibrahim, J. G. \u0026amp; Love, M. I. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e35\u003c/strong\u003e, 2084-2092, (2019).\u003c/li\u003e\n\u003cli\u003eDanecek, P.\u003cem\u003e et al.\u003c/em\u003e Twelve years of SAMtools and BCFtools. \u003cem\u003eGigascience\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, giab008, (2021).\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"nature-portfolio","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"","title":"Nature Portfolio","twitterHandle":"","acdcEnabled":false,"dfaEnabled":false,"editorialSystem":"ejp","reportingPortfolio":"","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"RNA sequencing, absolute quantification, gel-free library preparation","lastPublishedDoi":"10.21203/rs.3.rs-7256873/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7256873/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eNext-generation RNA sequencing (RNA-seq) is hampered by \u0026ldquo;primer dimer\u0026rdquo; (PD) artifacts and its quantitative performance reduced by polymerase fall-off (PF) at RNA modifications and secondary structures. Here we improved RNA-seq efficiency by incorporating (i) a post-reverse-transcription (RT) digestion of excess primers with \u003cem\u003eEscherichia coli\u003c/em\u003e exonuclease I for PD mitigation, thus obviating gel purification during RNA-seq library preparation, and (ii) a high-processivity reverse transcriptase to increase full-length reads. A full factorial experimental design was applied to absolute quantification RNA sequencing (AQRNA-seq), the most accurate NGS-based method for quantifying small RNAs, using cDNA libraries constructed from \u003cem\u003eE. coli\u003c/em\u003e small RNAs (\u0026gt;\u0026thinsp;85% tRNA) followed by sequencing, data processing, and data analysis. The novel PF and PD mitigation approaches increased AQRNA-seq sensitivity\u0026thinsp;\u0026gt;\u0026thinsp;10-fold by minimizing PF and maximizing target RNA reads. By increasing sensitivity and obviating gel electrophoresis for removing PD, AQRNA-seq and other NGS-based RNA-seq methods can now be automated to increase throughput and reduce RNA sample size.\u003c/p\u003e","manuscriptTitle":"Gel-free library preparation for next-generation RNA sequencing and small RNA quantification","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-10-08 10:50:36","doi":"10.21203/rs.3.rs-7256873/v1","editorialEvents":[],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"communications-biology","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"commsbio","sideBox":"Learn more about [Communications Biology](http://www.nature.com/commsbio/)","snPcode":"","submissionUrl":"","title":"Communications Biology","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"ejp","reportingPortfolio":"Communications Series","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"ece1be25-0697-4c6b-8043-e297a9793df6","owner":[],"postedDate":"October 8th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"under-review","subjectAreas":[{"id":55365199,"name":"Biological sciences/Biological techniques/Sequencing/RNA sequencing"},{"id":55365200,"name":"Biological sciences/Genetics/Epigenomics"},{"id":55365201,"name":"Biological sciences/Biochemistry/RNA"}],"tags":[],"updatedAt":"2026-04-16T12:27:07+00:00","versionOfRecord":[],"versionCreatedAt":"2025-10-08 10:50:36","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-7256873","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7256873","identity":"rs-7256873","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}