Structural bias in T4 RNA ligase-mediated 3′-adapter ligation

Ligated and random library preparation

Random RNA oligos were synthesized by Integrated DNA Technologies (Iowa, USA). To assess the sequence content of the random oligos, reactions adding a poly(A) tail to random RNA oligos G, C and U were performed using the protocol supplied by the manufacturer (New England Biolabs, Ipswich, MA, USA). Poly(C) tailing of random RNA oligo A was performed as previously described (31). All tailing reactions were incubated at 37°C for 2 h and the reactions were stopped by phenol/chloroform/isoamyl alcohol (IAA) extraction and precipitated by ethanol. After washing with 70% ethanol, the precipitated nucleic acid was resuspended in 20 µl H₂O prior to undergoing preparation for Ion Torrent sequencing. For the ligase selected libraries, each ligation reaction contained 1.4 µM ligase, 5.5 µM of adenylated SR1 adapter, 12% PEG8000, 50 mM Tris–HCl pH 7.5, 10 mM MgCl₂, 1 mM DTT and 50 ng of each random RNA oligo in a total volume of 200 µl. Reactions were incubated at 25°C for 2 h and stopped as described above.

Ion Torrent sequencing library preparation

Ligated products and tailed random oligos were reverse transcribed into cDNA using the ProtoScript™ M-MuLV Taq Reverse transcription (RT–PCR) Kit (New England Biolabs) following the protocol supplied with the kit. The primers for reverse transcription contained the Ion Torrent (Life Technologies, Carlsbad, CA, USA) ‘trP1’ sequence (CCTCTCTATGGGCAGTCGGTGAT) and sequence complimentary to the adapter sequence for ligated libraries or complimentary to the poly A or C tailed region for the random oligo libraries (Supplementary Table S1, Ligated RT, Random A,C,G or U RT). The cDNA products were amplified by 10 cycles of PCR using LongAmp® Taq master mix (New England Biolabs) with primers ‘IT Forward’, which added the Ion Torrent ‘A’ sequence, and ‘IT Reverse’ (Supplementary Table S1). PCR products were gel purified using either E-Gel® SizeSelect 2% agarose gels (Life Technologies) or 6% acrylamide gels. The purity and concentration of purified PCR products were analyzed using an Agilent 2100 Bioanalyzer (Agilent Biotechnologies, Santa Clara, CA, USA). Each purified library was diluted to a concentration of 44 pM and was prepared for Ion Torrent sequencing by using the Ion Xpress™ Template Kit (Life Technologies) and following the protocol supplied by the manufacturer. Libraries were sequenced on an Ion PGM™ using Ion 314™ or Ion 316™ chips deposited at full density.

Adenylation of DNA oligos

The DNA adapters were synthesized by Integrated DNA Technologies (Iowa, USA) with a phosphorylated 5′-end and a blocking amino group at the 3′-end. The adapters were adenylated using a 5′-DNA Adenylation Kit (New England Biolabs) as described previously (32). The adenylation reactions were stopped by adding 1 µg of proteinase K (New England Biolabs) per µl of adenylation reaction and incubated at 37°C for 30 min. DNA was further purified by two extractions with phenol/chloroform/IAA followed by ethanol precipitation. The adenylated oligos were separated from unadenylated ones in 20% Tris-borate-EDTA (TBE)–urea acrylamide gels. Bands corresponding to adenylated oligos were isolated, crushed and soaked in 1 ml water overnight at room temperature with constant rotation. After soaking, DNA was extracted from soaking solution using phenol/chloroform/IAA and precipitated by ethanol.

Ligation reactions

In vitro ligation reactions containing defined RNA substrates were carried out in 10 µl reactions containing 0.5 µM miRNA, 1 µM adenylated DNA adapter, 50 mM Tris–HCl pH 7.5, 10 mM MgCl₂, 1 mM DTT, 40 U of Murine RNase Inhibitor (New England Biolabs), 12.5% PEG8000 and 0.1 µM (Figures 1, 6 and 7) or 1.3 µM ligase (Figure 8). Reactions were incubated at 25°C for 2 h and stopped by adding same volume of 2×RNA loading buffer (95% formamide, 18 mM EDTA, 0.025% SDS, bromophenol blue, xylene cyanol). The ligation reactions were then loaded on a 15% TBE–urea gel to resolve ligated product, unligated RNA and unligated DNA adapter. The nucleic acid in the gel was stained with SYBR® Gold (Life Technologies) and scanned on a Typhoon™ 9400 Variable Mode Imager (GE Healthcare, NJ, USA). The intensity of each band was quantified using Quantity One software (BIO-RAD, Hercules, CA, USA) in order to determine the ligation efficiency. The amount of miRNA in the ligated product (I_{ligated miRNA}) was normalized using the following equation. I_{ligated miRNA} = I_ligated × length_miRNA/(length_miRNA + length_adapter). Ligation efficiency was calculated using the equation, ligation efficiency = I_{ligated miRNA}/I_miRNA.

Ligation reactions in the presence of small RNA mixtures contained 40 U of RNase Inhibitor (New England Biolabs), 50 mM Tris–HCl pH 7.5, 10 mM MgCl₂, 1 mM DTT, 12.5% PEG8000, 375 fmol of sRNA extracted from mouse ES cells (ES-E14TG2a, ATCC, Manassas, VA, USA), 0.75 fmol of 5′-³²P radio-labeled miRNA, 50 pmol of SR1 or SR1-R adapter, and 6.5 pmol of Rnl2tr in 10 µl reaction volume. Reactions were incubated at 25°C for 2 h and stopped as described above. The ligated and unligated radio-labeled miRNAs were separated in 15% TBE–urea gels. Gels were exposed to a storage phosphor screen (GE Healthcare, NJ, USA) and the intensities of bands were quantified using Quantity One software (BIO-RAD, Hercules, CA, USA). The ligation efficiencies of miRNAs were calculated from the equation, ligation efficiency = I_ligated/(I_ligated + I_unligated).

Bioinformatics

The 3′-adapter, or homopolymer tail sequences and 5′-constant regions were trimmed off using Galaxy (33–35). Only trimmed reads that were 21 nt in length were considered in subsequent analyses. RNA CONTRAfold (http://contra.stanford.edu/contrafold/index.html) was used for RNA secondary structure prediction. Default settings were used for the prediction. To predict RNA and adapter cofold structures, the Vienna RNAcofold (http://rna.tbi.univie.ac.at/cgi-bin/RNAcofold.cgi) was used with the default setting of minimum free energy algorithms and folding temperature at 25°C (36,37). The algorithm of Vienna RNAcofold prediction is based on the minimum free energy model (38). In our analysis, the 1999 Turner Model was used as the energy parameter during prediction (38).

Mouse ES cell small RNAs preparation

Total RNA was extracted from mouse embryonic stem cells (ES-E14TG2a, ATCC) using TRI Reagent (Sigma-Aldrich, MO, USA). The total RNA was subjected to Dnase I (New England Biolabs) digestion at 37°C for 30 min. RNA was further purified by acid-phenol:choloform (Life Technologies) extraction and ethanol precipitation. Pellets were resuspended in H₂O and the integrity of total RNA was assessed by checking the ribosomal RNAs in 1% agarose gels. Small RNAs <40 nt were isolated by flashPAGE™ fractionation (Life Technologies) and precipitated by 1.5 volume of isoproponal, 1/10 volume of 3 M sodium acetate and 25 µg of linear acrylamide (Amresco, Solon, OH, USA). After precipitation and washing, the small RNAs were resuspended in water and the RNA concentration was determined by Qubit Fluorometer (Life Technologies). Each pmol of small RNAs were further treated with 1 U of alkaline phosphatase, Calf Intestinal (CIP, New England Biolabs) at 37°C for 1 h. The reaction was stopped by acid–phenol:choloform extraction and small RNAs were collected by ethanol precipitation. After resuspension in H₂O, the concentration of purified small RNAs was determined by Quibit Fluorometer (Life Technologies).

The efficiency of 3′-adapter ligation varies between different miRNAs

To illustrate the variation in 3′-adapter ligation efficiency for different miRNAs, we selected 25 miRNAs from miRBase (30) and performed ligations using Rnl2tr (Figure 1A). As shown in Figure 1B, the ligation efficiency of a given miRNA to a pre-adenylated DNA adapter is highly variable, ranging from ligation of as little as 0% of the input to as much as 100%. These results confirm that there is significant bias in miRNA 3′-adapter ligation reactions that cannot be easily explained by miRNA primary sequence or predicted secondary structure (Supplementary Table S2).

Ligation and HTS of oligonucleotide pools to study ligation bias

To study the bias of T4 RNA ligases, we designed an in vitro ligation selection assay that uses a pool of random RNA oligos to which a 3′-adapter is ligated (Figure 2). Following reverse transcription and amplification, the sequence of the ligated products was determined using the Ion Torrent sequencing platform (29). The random RNA oligo pool contained equimolar amounts of four random RNA oligos, which consisted of the same constant 21 nt region at the 5′-end, followed by a 20 nt random region, and a U, C, G or A at the 3′-end (Supplementary Table S1). A fixed nucleotide at the 3′-end of oligo is required for oligo synthesis. Therefore, it was necessary to hand mix an equimolar amount of four random oligos to generate an oligo pool containing 21 random positions.

To assess the frequency of nucleotides at each randomized position in the oligo pool, the random RNA oligos were subjected to Ion Torrent sequencing without undergoing adapter ligation (Figure 2). Each oligo was tailed with poly nucleotides using poly(A) polymerase. The random RNA oligo ending with U, C or G was tailed with ribonucleotide A, while the random RNA oligo ending with A was tailed with ribonucleotide C in order to distinguish the last A from the tailed region. We performed poly-adenylation reactions under conditions that minimized 3′-end nucleotide bias (39). These conditions resulted in essentially complete tailing of the input oligos (data not shown). After tailing, the RNAs were reverse transcribed and the cDNAs were used for Ion Torrent sequencing library preparation. After sequencing, 29 737 sequences from each of the four random oligos were pooled to generate a random input library with 118 948 reads in total (Supplementary Table S3).

Our ligation selection reactions used 9.0 × 10¹³ molecules of random oligos in order to cover all the 4.4 × 10¹² possible sequence combinations from 21 random positions. For each library, a ligation reaction was performed using an excess of RNA ligase and 5′-pre-adenylated DNA adapter (SR1) over random RNA oligos. The SR1 adapter was blocked at its 3′-end by an amino group to prevent it from participation intra- and inter-molecular ligation with other SR1 molecules (25,28,39). Ligated products were reverse transcribed into cDNA and sequenced on the Ion PGM™. Rnl1 and four variants of Rnl2tr were tested in the selection assay. After quality control and adapter trimming, we obtained 10⁵–10⁶ sequences for each library. The summary of sequence reads and quality control of libraries are shown in Supplementary Table S3.

T4 RNA ligases do not show appreciable bias in RNA substrates at the primary sequence level

We first looked for evidence that the ligases have any primary sequence preference within the 21 nt random region of the RNA substrates. To do so, we calculated the frequency of each nucleotide at each position from 10⁵ to 10⁶ sequences in each library (Supplementary Table S4). In Figure 3A and Supplementary Figure 1, the raw nucleotide frequency each position in random region was plotted in enoLOGO format (40). The frequency of a particular nucleotide at each position is proportional to the size of the letter representing that base. The distribution of nucleotide frequencies in ligated libraries is very similar to that of the random input library. For the random input library, the nucleotide frequencies from positions 1 to 19 shares a similar distribution with A and U being close to 25%, but slightly less C (19–20%) and higher G (30–34%) (Supplementary Table S4). The frequencies at position 20 showed a slightly different distribution pattern compared to positions 1–19 with 22% A, 22% C, 29% G and 26% U. The nucleotide percentage of position 21 was 25% for all 4 nt reflecting equal numbers of sequences that were combined from four random oligo libraries.

We then normalized the nucleotide frequency in the ligated libraries to that of the random input library and determined the enrichment of nucleotides at every position using a previously described method (41). Briefly, the relative ratio of each nucleotide n at position p in each ligase selected library was determined from the nucleotide frequency in the ligated RNA pool versus the frequency in the random input pool by using the equation:

The RN_np ratio at each position was normalized to 1 by equation:

The value of (RN_np −0.25) was plotted according to the nucleotide positions for each ligase (Figure 3B). If (RN_np −0.25) of a nucleotide n at position p is equal to 0, it indicates the ligase does not have any preference for nucleotide n at position p. If (RN_np −0.25) is greater or less than 0, it means that nucleotide n is preferred or not preferred at position p, respectively. As shown in Figure 3B, Rnl1 and the four variants of Rnl2tr show minimal preference for any particular nucleotide at any particular position in our RNA substrates. We interpret these observations to mean that the primary sequence of RNA substrate has minimal impact on the ligation efficiency.

Ligation efficiency is affected by the secondary structure within an RNA substrate

Given the striking differences we observed in miRNA 3′-adapter ligation efficiency and having failed to observe significant sequence bias in our sequencing experiments, we next asked whether T4 RNA ligases prefer certain secondary structures within RNA substrates as suggested previously (25). Considering the function of T4 Rnl1 is to repair a break in the anticodon loop of tRNA^Lys, it is possible that ligases such as Rnl1 prefer RNA substrates similar in structure to its biological substrate. To explore the possible secondary structure preferences of ligases, sequences from each library were analyzed by CONTRAfold to predict their secondary structures (42). After folding, structural predictions of the 42 nt sequences were sorted into groups based on the number of unpaired nucleotides at their 3′-end (Figure 4). If the 3′-end nucleotide is predicted to be paired, there are ‘0’ unpaired nucleotides and if no pairing was predicted within the RNA, there are ‘42’ unpaired nucleotides. The percentage of each group in a specific library was calculated and the enrichment of that structure group was determined by comparing its percentage to that in the random input library using the following equation, ‘(Observed–Expected)/Expected’. ‘Observed’ is the percentage of a specific group in a ligated library and ‘Expected’ is its corresponding percentage in the random input library. If the calculated value is positive or negative, it means the structure is over- or under-represented in the ligated library, respectively. As shown in Figure 4, it is clear that all tested ligases share a similar preference. Specifically, RNAs with fewer than three unpaired nucleotides at the 3′-end are under-represented in the ligated libraries while RNAs with three or more unpaired nucleotides at the 3′-end appear at their expected frequency or are over-represented. Thus, RNA ligases generally prefer RNA with a relatively accessible 3′-end and one source of bias in ligation is secondary structure at the 3′-end of an RNA substrate.

Analysis of RNA and adapter cofolding

When comparing predicted structure at the 3′-end to the ligation efficiency of an RNA (Figure 1 and Supplementary Table S2), it is clear that the preference of ligases for RNAs lacking 3′-end secondary structure cannot completely explain the ligation bias we observed. We hypothesized that another possible source of bias could result from the interaction between the RNA substrate and the adapter. Thus, we predicted the cofold structures of our sequenced library members with the SR1 adapter to examine the correlation between cofold structures and ligation efficiencies.

Using the Vienna RNAcofold algorithm (36,37), we cofolded all sequences from all sequenced libraries with the SR1 adapter and classified the cofold structures according to the regional secondary structure at the ligation junction. These structures are summarized and presented in dot-bracket notation and schematic representation (Figure 5A). Brackets indicate paired nucleotides and dots represent unpaired nucleotides. An ampersand symbol, ‘&’, denotes the ligation junction. Structure class 2, 6, 10 and 14 are distinct from the rest of the classes because the RNA and adapter do not form heterodimers.

Only 8 out of the 16 structure classes were present in all sequenced libraries as shown in the distribution plot in Figure 5B. Structure classes 5, 7 and 13 are the three most abundant and they make up 19.2–33.1% of the total number of classes in each library. Structure classes 1, 3, 9, 11 and 15 were less abundant and ranged from 1.1% to 8.6% in each library. We did not observe structure classes 2, 4, 6, 8, 10, 12, 14 and 16 in any of our sequenced libraries.

To assess the cause of the absent structure classes in our sequenced libraries, and exclude that the absence was caused by a lack of sequence coverage, we generated random libraries in silico. Ten simulated random libraries, each consisting of 300 000 random sequences, were generated by a computer containing the same 5′-constant sequence and 21-nt random region. Each sequence was then cofolded with the SR1 adapter and the average distribution of cofold structures from the 10 libraries was compared to our sequenced random input library (Figure 5B). Distribution of the simulated random libraries is similar to that of our sequenced random library confirming that our sequenced random library indeed represents a random pool. The absence of observed cofold structure classes is therefore not due to the sequence coverage nor bias introduced by poly(A) polymerase. In addition, the absent cofold structure classes do not result from the presence of 5′-constant region in the RNA (Supplementary Figure S2) since the presence or absence of the 5′-constant region did not affect their presence in our cofold predictions. The absent cofold structure classes are structures where the RNA and adapter are predicted to not form heterodimers or where adapters are internally structured. Since the SR1 adapter does not form internal secondary structure or homodimers according to secondary structure prediction (Figure 6A), it makes sense that these structure classes would not be represented when RNAs are cofolded with the SR1 adapter. These observations lead us to conclude that the distribution of cofold structures is largely due to the SR1 adapter and limitations in what structures that it can form.

Enrichment of RNA-adapter cofold structures in the ligated libraries

To determine whether the ligases prefer or discriminate against a specific cofold structure, we first compared the percentage of each cofold structure class in the ligated libraries to that in the random library. We further calculated enrichment for each cofold structure class using ‘(Observed–Expected)/Expected’, in which ‘Observed’ is the percentage of a specific cofold structure class from a ligated library and ‘Expected’ is the percentage of the corresponding class from the random library. When the value is greater or less than 0, it means the structure class is over- or under-represented, respectively. If the value is equal to 0, it means there is no preference for that structure class. As shown in Figure 5C, it is evident that two structure classes, 7 and 13, were slightly over-represented in the libraries of all tested ligases and none of the ligases showed much preference toward structure class 11. Three structure classes, 3, 5 and 9, were found to be under-represented in the ligated libraries, though structure 5 was only slightly so. The natural substrate of T4 Rnl1, cleaved tRNA^Lys, is predicted to cofold into structure class 5 (5). Interestingly, the ligases did not differ in their preference for structure class 5 as reflected by the representation of this class in our sequenced libraries. However, the ligases showed different preferences toward structure classes 1 and 15. For instance, class 15 was under-represented in the library using T4 Rnl1 but was over-represented when using Rnl2tr, but all of the mutant variants of T4 Rnl2 showed similar preferences to structure class 15 as Rnl1.

RNA-adapter cofold structure influences ligation efficiency

Our observations that certain cofold structure classes were over- or under-represented in our sequenced libraries prompted us to examine the influence of RNA-adapter cofolding on ligation efficiency. We designed a new adapter, SR1-S, which shares the same sequence as SR1 at the 5′-end from positions 1 to 12, but with a modified 3′-sequence so that, in contrast to SR1, SR1-S is predicted to have secondary structure. We presume that internal secondary structure in the adapter would reduce its ability to productively cofold with RNAs (Figure 6A and Supplementary Table S5). In Figure 6B, SR1-S migrated faster than the SR1 adapter, which is the same length, likely because of incomplete denaturation or renaturation. SR1-S stained more strongly than an equivalent amount of SR1 in gels, consistent with the 2× increased staining of dsDNA versus ssDNA with SYBR Gold (43). We then compared the ligation efficiencies of 18 miRNAs with both the SR1 and SR1-S adapter. As shown in Figure 6B and C, we observed large overall decreases in ligation efficiency using SR1-S. The average ligation efficiency of the miRNAs decreased from 27% with SR1 to 8% with SR1-S. The decrease in ligation was most dramatic for miRNAs which ligated efficiently with SR1, such as miR-31, let-7, miR-139-3p, miR-122, miR-253, miR-1614 and miR-2419. The fold decrease in ligation efficiency for these RNAs ranged from 3- to 22-fold. miRNAs that had low ligation efficiencies with SR1 also ligated poorly with SR1-S. The finding that we can modulate ligation efficiency by changing secondary structure within the adapter without changing primary sequence at the 5′-end is at odds with a previous report (26). This report concluded that RNA ligase bias is due to primary sequence-specific preferences at the ligation junction, specifically the first two nucleotides of the 5′-end of adapter sequence and the last two nucleotides of the 3′-end of the RNA (26). Our results, especially when considered with the nucleotide frequency results from our sequencing experiments (Figure 3), contradict this conclusion. Together, our results demonstrate that structures within and between the RNA acceptor and the adapter are more important than primary sequence in influencing ligation efficiency. These results support a model where favorable heterodimeric RNA and adapter cofold structures promote efficient ligation.

To further test our model, we tried to improve the ligation efficiency of seven miRNAs that ligate poorly with SR1 by designing a new adapter for each miRNA. In each case, the adapter was designed so that the predicted RNA-adapter cofold structure was changed from an under-represented class to an over-represented class (Table 1). The predicted cofold structures of six miRNAs with SR1 belonged to one of two under-represented cofold structure classes, either class 1 or 5. The adapters we designed for each of these miRNAs changed their predicted RNA-adapter cofold structures to the over-represented structure class 13. The seventh miRNA, miRNA-5183, was already predicted to cofold with SR1 in structure class 13, despite the fact that it ligates poorly with SR1. When we examined predicted loop size of miRNAs within structure class 13, and correlated this with ligation efficiency, we noted a trend that some loop sizes appear to be unfavorable for ligation (data not shown). The adapter we designed for miR-5183 adjusted the size of the loop in the cofold structure class 13.

Overall, the average ligation efficiency of these miRNAs increased from 2.7% with SR1 to 18.7% with the newly designed adapters (Figure 7 and Supplementary Table S1). The increase in ligation efficiency for miRNAs ranged from 2.2- to 17.1-fold. We interpret these results to indicate that a RNA-adapter pair that is predicted to cofold in an under-represented class is unfavorable for ligation. These data confirm the important role of RNA-adapter cofolding during ligation. We suggest that the method demonstrated here is useful for improving the ligation of any known RNA sequence and could be applied to situations where accurate quantification of a set of known RNAs is important.

Improved miRNA ligation efficiency using 5′-end randomized adapters

Our results support the hypothesis that cofold structures are a major factor contributing to T4 RNA ligase bias. We therefore attempted to minimize the bias using a pool of adapters with randomized 5′-ends in order to decrease the likelihood that a particular miRNA will be incompetent for ligation with a single sequence adapter. In other words, we attempt to minimize bias by supplying many adapters to increase the likelihood of more favorable cofold structures between adapters and all RNAs in the sample. We designed a pre-adenylated adapter, SR1-R, which contains the same sequence as SR1 except the first six nucleotides at 5′-end are randomized (Supplementary Table S1).

We performed ligation reactions using Rnl2tr to compare the ligation efficiency of each miRNA with SR1 and SR1-R (Figure 8A). miRNAs that ligated efficiently with SR1 were also efficiently ligated with SR1-R, while miRNAs that ligated poorly with SR1 generally showed improved ligation efficiency with SR1-R (Figure 8B). Reflecting this observation, the average ligation efficiency increased from 67% with SR1 to 78% with SR1-R. Furthermore, the index of dispersion, defined as the ratio of the variance to the mean, decreased from 0.13 with SR1 to 0.034 with SR1-R. The decreased index of dispersion indicates less divergence in ligation efficiencies among miRNAs and suggests that using a randomized SR1-R adapter reduces the bias of Rnl2tr in ligation.

While the ligations of pure miRNA substrates demonstrated the potential benefits of using randomized adapters, we further examined whether the same benefits were retained for a particular miRNA when a pool of small RNA substrates is present. In order to test this, we radio-labeled an miRNA of interest at the 5′-end with ³²P so that we could use small amounts (0.75 fmol) of the miRNA in a mixture containing small RNAs extracted from mouse embryonic stem (ES) cells (Figure 9A). Ligation reactions were performed using 500-fold excess of mouse ES cell small RNAs and excess amount of adapter compared to radio-labeled miRNA. The ligated and unligated radio-labeled miRNA were separated in 15% TBE–urea gels (Figure 9B). The average ligation efficiency increased from 36% with SR1 to 46% with SR1-R. The ligation efficiency was increased for nine miRNAs, unchanged for 13 miRNAs and somewhat decreased for two miRNAs (Figure 9C). All miRNAs with <40% ligation efficiency with SR1 exhibited improved or unchanged ligation efficiencies when ligated to SR1-R. In agreement with the results of the experiments with pure miRNAs, the ligation bias was reduced as measured by the index of dispersion, which decreased from 0.42 with SR1 to 0.23 with SR1-R. In summary, the use of randomized adapters when ligating adapters to the 3′-end of miRNAs with Rnl2tr results in generally improved ligation efficiency and decreased ligation bias.