Hopeful monsters: unintended sequencing of famously malformed mite mitochondrial tRNAs reveals widespread expression and processing of sense-antisense pairs.

Although tRNA structure is one of the most conserved and recognizable shapes in molecular biology, aberrant tRNAs are frequently found in the mitochondrial genomes of metazoans. The extremely degenerate structures of several mitochondrial tRNAs (mt-tRNAs) have led to doubts about their expression and function. Mites from the arachnid superorder Acariformes are predicted to have some of the shortest mt-tRNAs, with a complete loss of cloverleaf-like shape. While performing mitochondrial isolations and recently developed tRNA-seq methods in plant tissue, we inadvertently sequenced the mt-tRNAs from a common plant pest, the acariform mite Tetranychus urticae, to a high enough coverage to detect all previously annotated T. urticae tRNA regions. The results not only confirm expression, CCA-tailing and post-transcriptional base modification of these highly divergent tRNAs, but also revealed paired sense and antisense expression of multiple T. urticae mt-tRNAs. Mirrored expression of mt-tRNA genes has been hypothesized but not previously demonstrated to be common in any system. We discuss the functional roles that these divergent tRNAs could have as both decoding molecules in translation and processing signals in transcript maturation pathways, as well as how sense-antisense pairs add another dimension to the bizarre tRNA biology of mitochondrial genomes.

INTRODUCTION

As the adapter molecule between mRNA codons and amino acids, transfer RNAs (tRNAs) are a conserved feature of life. The vast majority of known tRNAs have a uniform structure comprised of acceptor and anticodon stem, D- and T-arms, a variable region and accompanying loops to form the easily recognizable clover-leaf secondary structure. An increasing number of exceptions to this structural uniformity has come from the expanding field of mitochondrial genetics. Most animal mitochondrial genomes encode 22 tRNA genes considered sufficient to decode all codons (1,2) (but see (3,4)). However, some mitochondrial tRNAs (mt-tRNAs) can deviate significantly from the canonical tRNA structure. First identified in mammalian mitochondria in the late 1970s, a tRNA-Ser lacking the D-arm has now been found to be very common in metazoan mitochondria (5–8). Following the mt-tRNA-Ser discovery, even more extensive cases of tRNA truncation were reported in nematode mitochondria where all of the mt-tRNAs lacked either the D- or the T-arm, which were lost in favor of short replacement loops consisting of only 5–8 nt (9,10). The most extreme examples of tRNA truncation have come from a class of nematodes (Enoplea) and a superorder of arachnids (the acariform mites) (11–13). Aside from the mammalian mt-tRNA-Ser lacking the D-arm, there is very limited experimental evidence for expression or aminoacylation activity in translation of these aberrant tRNA structures (2,4). The most notable support comes from the sequencing of several mt-tRNAs from the nematode Romanomermis culicivorax, including an armless tRNA that was shown to be post-transcriptionally modified with a CCA-tail (11). In the Acariformes, mt-tRNA expression data are lacking and (unlike in Enoplea) mt-tRNAs in Acariformes are not well conserved, with some tRNA genes differing even at an intraspecific level (13). The highly divergent and aberrant structure of mt-tRNAs in acariform mites poses questions about the limits of a minimal functional tRNA (12,14–15).

Mites are among the oldest and most diverse groups of terrestrial animals (16). As one of the two superorders of mites, the Acariformes include well-known representatives such as dust mites, scabies, chiggers and spider mites. Partially due to both agricultural and medical concerns, there has been increasing interest in the mitochondrial genomics of acariform mites for phylogenetics as well as development of pest control agents and treatments (15,17). One of the most notable features discovered from the sequencing of Acariformes mitochondrial genomes has been the widespread rearrangement and extreme truncation of mt-tRNAs (13–14,18–20). Going beyond the frequently described loss of a D- or T-arm, mt-tRNAs from multiple taxa in Acariformes have been predicted to completely lack both arms and produce just a simple stem-loop structure (13). In other species, no structure identifiable as a tRNA could be found for some anticodons, and they were reported as functionally lost (12,18,21–22). In addition to their extremely short lengths, some predicted acariform mt-tRNAs have mismatches in acceptor stems (12–13,20), which appear to be widespread in diverse arachnid lineages and often occur at the interface of overlapping tRNA genes (23,24). It is possible that these mismatches are edited and corrected post-transcriptionally as observed in some other eukaryotes (25–27). Regardless of whether they are edited, these mismatches create additional difficulties when trying to annotate tRNA genes and predict structure and function.

Programs like tRNAscan-SE (28) and ARWEN (29) that are used to infer tRNA gene presence from genomic data frequently fail to detect the highly degenerate mt-tRNAs in Acariformes taxa. The true loss of a mt-tRNA gene would likely require that a nuclear-encoded tRNA be imported into the mitochondrial matrix to as a functional replacement (30) – a phenomenon that is ubiquitous in some eukaryotic lineages (e.g. plants (4,31)) but thought to be rare in metazoans (2). Doubts about the extensive import of tRNAs into animal mitochondria have led to the revisiting and manual inspection of noncoding space in multiple Acariformes mitochondrial genomes, resulting in predictions for additional, but sometimes not all, mt-tRNA genes (12–13,15,22). In a more distantly related Sarcoptidae mite, the partial sequencing of two truncated mt-tRNA genes thought to be lost (32) suggests that genomic data alone is often insufficient to detect aberrant tRNA genes in the mitochondrial genomes of multiple lineages.

An additional complication of proposed mt-tRNA losses relates to the multi-functionality of tRNAs in mitochondrial gene expression. Mt-tRNAs have long been thought to be processing signals for mRNA and rRNA transcript maturation. Under what is known as the tRNA punctuation model, protein-coding and rRNA genes are often flanked by one or more tRNA genes in mitochondrial genomes, effectively co-opting the recognition and cleavage by the tRNA-interacting enzymes RNase P and RNAse Z to produce excised mRNA or rRNAs from polycistronic mitochondrial transcripts (33,34). Mt-tRNAs have also been implicated as origins of mitochondrial DNA replication, where hybridization between the tRNA and the complementary DNA gene sequence initiates replication (35,36). Thus, the complete loss of some mt-tRNA genes may necessitate additional evolutionary mechanisms to maintain functional mitochondrial gene expression.

Extensive work had been done to predict tRNA candidates as well as assign functionality solely based on genomic data in Acariformes mites (12), but there remains a lack of experimental data for their expression and processing. This scarcity of expression data may be partly attributed to longstanding technical difficulties in sequencing tRNAs (tRNA-seq), due to tRNAs being incredibly recalcitrant to reverse transcription (37), a step necessary in the vast majority of high-throughput RNA-seq methods. tRNAs are the most extensively modified RNAs known (38), and chemical modifications at the Watson–Crick face of tRNA bases can interfere with reverse transcriptase activity. The stalling, skipping, or disassociation of a reverse transcriptase can result in termination of cDNA synthesis or the misincorporation of incorrect nucleotides in the cDNA at the corresponding modified base positions (39–41). Additionally, the tightly base-paired 3′- and 5′-termini of tRNAs can inhibit adapter ligation, which is necessary for the priming of reverse transcription in RNA-seq library preparation (42). A major breakthrough in tRNA-seq came from the application of the demethylating enzyme AlkB to remove certain reverse transcription-inhibiting modifications prior to tRNA-seq library construction (43,44). Additional improvements in the sequencing of full-length tRNA molecules included the development of adapters that take advantage of the unique structure of mature tRNAs (45,46). All mature tRNAs have protruding, unpaired bases at the 3′-termini, composed of the discriminator base and the trinucleotide CCA (47). This CCA ‘tail’ is post-transcriptionally added by tRNA nucleotidyltransferases and has been considered a hallmark of tRNA functionality as the site of aminoacylation (48). YAMAT-seq is a recently developed technique based on adapters that are complementary to these protruding nucleotides, aiding adapter ligation and resulting in the capture and sequencing of full-length, mature tRNAs (46).

Here, we report the sequencing of mt-tRNAs from the acariform spider mite and common plant pest Tetranychus urticae, which was the inadvertent outcome of applying mitochondrial enrichment techniques and targeted tRNA-seq methods to infested plant tissue. This fortuitus capture of the famously divergent mt-tRNAs from T. urticae provides unprecedented insight into the expression and post-transcriptional modifications of some of the most aberrant tRNAs ever found.

MATERIALS AND METHODS

Mitochondrial isolation and RNA extraction

Mitochondrial isolations were performed on T. urticae infested leaf tissue from the angiosperm Silene vulgaris. Plants were grown in a greenhouse with supplemental lighting (16-h/8-h light/dark cycle) in the Colorado State University greenhouse, Fort Collins, CO, USA. All mitochondrial isolations were done in triplicate. Each replicate used 75 g of leaf tissue from a dozen 8-month-old plants. Leaf tissue was disrupted in a Nutri Ninja Blender for 2 × 2 s short bursts, and 1 × 4 s blending in 350 ml of a homogenization buffer containing: 0.3 M sucrose, 5 mM tetrasodium pyrophosphate, 2 mM ethylenediaminetetraacetic acid, 10 mM KH2PO4, 1% PVP-40, 1% bovine serum albumin, 20 mM ascorbic acid and 5 mM cysteine, pH 7.5-KOH.

Differential centrifugation was performed to remove nuclei, plastids and cellular debris in a Beckman Avanti JXN-30 centrifuge with a JS-24.38 swinging bucket rotor with the following centrifugation steps: i) 10 min at 500 g with max brake, ii) 10 min at 1500 g with max brake, iii) 10 min at 3000 g with max brake. After each centrifugation step, the supernatant was transferred into a clean centrifuge tube. Mitochondria were then pelleted by centrifugation for 10 min at 20 000 g with the brake off. Supernatant was discarded and the pellet was resuspended using a goat-hair paintbrush and 2 ml wash buffer containing 0.3 M sucrose, 10 mM MOPS, 1 mM EGTA, pH 7.2-KOH. A total 30 ml wash buffer was then added to the resuspended pellet and the homogenate was centrifuged for 5 min at 3000 g. Mitochondria were once again pelleted by centrifugation for 10 min at 20 000 g. The pellet was resuspended with 500 μl wash buffer and paint brush. The mitochondrial supernatant was then added to a glass Dounce homogenizer and homogenized with three strokes.

Homogenized mitochondria were then suspended on top of a Percoll gradient with the following Percoll density layers, 18, 25, 50%. The gradient was then centrifuged at 40 000 g for 45 min with the brake off. The mitochondrial band at the 25%:50% interface was then aspirated off of the gradient and diluted with 30 ml of wash buffer. The diluted mitochondria were then centrifuged at 20 000 g for 10 min. The supernatant was vacuum aspirated, and the mitochondrial pellet was resuspended in a fresh 30 ml of wash buffer and centrifuged at 10 000 g for 10 min. The supernatant was vacuum aspirated, and the mitochondrial pellet was resuspended in 1000 μl of fresh wash buffer. Resuspended mitochondria were centrifuged at 10 000 g for 10 min. Supernatant was removed with a pipette and the mitochondrial pellet immediately went into RNA extraction procedures, using 1000 μl TRIzol following TRIzol manufacturer's RNA extraction protocol.

AlkB treatment

Total mitochondrial RNA from each isolation was divided into two treatments: untreated and demethylated with the enzyme AlkB. AlkB reactions were performed using a modified version of existing protocols (43–44,49). Demethylation was performed by treating 6 μg of total mitochondrial RNA with 250 pmols of AlkB in a reaction volume of 50 μl containing: 70 μM ammonium Iron(II) sulfate hexahydrate, 0.93 mM α-ketoglutaric acid disodium salt dihydrate, 1.86 mM ascorbic acid and 46.5 mM HEPES (pH 8.0), incubated at 37°C for 60 min. The reaction was followed by a phenol–chloroform RNA extraction and ethanol precipitation with the addition of 0.08 μg of RNase-free glycogen, and resuspension in water. RNA integrity was checked on a TapeStation 2200.

YAMAT adapter ligation and Illumina library construction

All adapter and primer sequences used in library construction can be found in Supplementary Table S1. In order to remove amino acids from the mature tRNAs (deacylation), demethylated and untreated RNA was incubated in 20 mM Tris HCl (pH 9.0) at 37°C for 40 min. Following deacylation, adapter ligation was performed using modified protocols from (46,50). A 9 μl reaction volume containing 1 μg of deacylated RNA and 1 pmol of each Y-5′ adapter (4 pmols total) and 4 pmols of the Y-3′ adapter was incubated at 90°C for 2 min. A total of 1 μl of an annealing buffer containing 500 mM Tris–HCl (pH 8.0) and 100 mM MgCl2 was added to the reaction mixture and incubated for 15 min at 37°C. Ligation was performed by adding one unit of T4 RNA Ligase 2 enzyme (New England Biolabs) in 10 μl of 1× reaction buffer and incubating the reaction at 37°C for 60 min, followed by overnight incubation at 4°C.

RT of ligated RNA was performed using SuperScript IV (Invitrogen) according to the manufacturer's protocol. Briefly, 1 μl of 2 μM RT primer, and 1 μl of 10 mM dNTP mix was added to 11 μl of the deacylated RNA from each sample. The mixture was briefly vortexed, centrifuged and incubated at 65°C for 5 min. Then, 4 μl of 5× SSIV buffer, 1 μl 100 mM DTT, 1 μl RNaseOUT and 1 μl of SuperScriptIV were added to each reaction. The mixture was then incubated for 10 min at 55°C for RT and inactivated by incubating at 80°C for 10 min. The resulting cDNA was then amplified by polymerase chain reaction (PCR) in a 50 μl reaction containing 7 μl of the RT reaction, 25 μl of the NEBNext 2× PCR Master Mix, 1 μl of the PCR forward primer, 1 μl of the PCR reverse primer and 15.5 μl dH2O. Ten cycles of PCR were performed on a Bio-Rad C1000 Touch thermal cycler with an initial 1 min incubation at 98°C and 10 cycles of 30 s at 98°C, 30 s at 60°C and 30 s at 72°C, followed by 5 min at 72°C.

Size selection of the resulting PCR products was done on a BluePippin (Sage Science) with 3% agarose gel cassettes and marker Q3 following the manufacturer's protocol. The size selection parameters were set to a range of 180–231 bp, with a target of 206 bp. Size-selected products were then cleaned using solid phase reversable immobilization beads and resuspended in 10 mM Tris (pH 8.0). Libraries were dual-indexed and sequenced on an Illumina NovaSeq 6000 with paired-end, 150-bp reads.

Raw sequencing reads are available via the NCBI Sequence Read Archive under BioProject PRJNA662108. Processed reads, as well as final mapped counts and counts per million for each tRNA/stem-loop, can be accessed through NCBI’s Gene Expression Omnibus (51) under GEO Series accession number GSE162913 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE162913).

YAMAT-seq read processing and mapping

Reads were trimmed with Cutadapt (ver.1.16 (52)) using the following options: -q 10 –discard-untrimmed –nextseq-trim = 20. Forward and reverse trimmed reads were merged using BBMerge (53) (BBTool software package), with a minimum overlap of 20 bp and 0 mismatches. Identical reads were summed and collapsed into read families using the FASTQ/A Collapser tool from the FASTX-Toolkit version 0.0.13 (http://hannonlab.cshl.edu/fastx_toolkit/index.html). The mapping of reads to reference tRNA gene set (described below) was done with previously published pipeline and Perl scripts (50). The number of reads with a CCA tail was also calculated with a custom Perl script. Modified sites that differed from the reference were determined with a pre-existing pipeline (50), except MUSCLE (ver. 3.8.31, (54)) was used as the alignment tool. Pipeline scripts can be found at https://github.com/warrenjessica/YAMAT-scripts.

The final numbers of mapped reads per tRNA/stem-loop per library and expression for each transcript in counts per million reads can be found in Supplementary Table S2.

Extraction and taxonomic identification of Tetranychus reads

While applying a pipeline to remove commonly occurring environmental contaminants (e.g. soil bacteria) it became apparent that a considerable number of reads mapped to Tetranychus mitochondrial genome accessions deposited on GenBank. Preliminary analysis found that T. urticae (Genbank accession KJ729022) was the most common species. In order to extract all possible reads originating from T. urticae, reads were BLASTed (blastn, e-value of < 1e-3, low complexity regions not filtered [-dust no]) against the T. urticae mitochondrial genome and only those that produced a hit were retained for further filtering. To remove any tRNAs originating from other organisms, reads that hit to the T. urticae mitochondrial genome were then BLASTed against the GenBank nucleotide database (downloaded 28 January 2020), and the taxonomic hit information for the top hit for each read was extracted. Only reads in which the top hit had a taxonomic assignment containing the genus Tetranychus were retained. The entire mitochondrial genome of T. urticae has been included as a contaminating contig in the whole genome assembly of the flowering plant Dioscorea rotundata (Genbank accession LC219377). Therefore, reads with a top hit to this D. rotundata accession were also retained.

To determine the Tetranychus species and strain that was sequenced, the filtered Tetranychus reads were mapped to the eight NCBI RefSeq Tetranychus reference mitochondrial genomes (HM753535, KJ729017, KJ729018, KJ729019, KJ729020, KJ729021, KM111296 and NC_010526) as well the green strain of T. urticae (KJ729022) and the red strain of T. urticae (KJ729023) published by Chen et al. (13) Of these reads, the overwhelming majority (94%) showed a best hit (or a tied best it) to T. urticae strain green (KJ729022). The small minority of reads that had a better match to another Tetranychus species or strain of T. urticae (i.e. strain red) were also very close hits to T. urticae strain green and appeared to reflect variation at sites subject to reverse transcription misincorporations (see modification index methods below). This is consistent with T. urticae being a very common greenhouse pest in Colorado (55), strongly indicating that T. urticae was the primary or sole contributing Tetranychus species. Filtered reads were then mapped back to a final reference set of T. urticae strain green mt-tRNA genes (reference set determined below) requiring that any mapped read represented ≥60% coverage of a reference tRNA and had an e-value < 1e-3.

tRNA coordinate determination

The degenerate nature of mt-tRNAs in Acariformes has resulted in historical difficulties in finding and defining the coordinates (i.e. the start and stop position) of tRNA genes. In order to construct a reference set of tRNAs based on tRNA transcription data, the ‘closest’ function in bedtools (ver. 2.27.1, (56)) was used to assign reads to a gene based on any amount of overlap with a previously annotated tRNAs. Reads that had overlap with more than one tRNA gene were assigned to the gene with which the greatest overlap occurred. The most frequent start and stop position for each gene were used as the final reference coordinates for mapping (Supplementary Table S3). Additional transcribed stem-loops were found with the same protocol by identifying all reads that did not overlap with a previously annotated tRNA gene. Final references for all genes can be found in Supplementary Table S4.

Folding predictions

Folding models for tRNA-Asn, tRNA-His, tRNA-Leu2 and tRNA-Lys were generated using tRNAscan-SE 2.0 (search options -O -C -X 0.1 (28)). All other tRNAs were not detected by tRNAscan-SE. The remaining tRNA models were based on Chen et al., if available (13). Antisense tRNAs, stem-loops and tRNA-Phe and tRNA-Thr, which were sustainably shorter than previous predictions, were folded manually to maintain arm structure. Supplementary Table S5 lists multiple folding predictions including maximum free energy, centroid, tRNAscan-SE and Chen et al. (13).

RESULTS

The highly degenerate mt-tRNA complement of T. urticae is expressed and post-transcriptionally processed

The application of a gradient-based mitochondrial isolation method on leaf tissue from the angiosperm S. vulgaris resulted in the co-enrichment of the mitochondria from common plant pests inhabiting the tissue alongside the intended plant mitochondria. When these mitochondrial-enriched fractions were subjected to targeted tRNA-seq methods that preferentially capture mature tRNAs (three biological replicates processed both with and without AlkB enzymatic treatment), they produced a substantial number of reads mapping to the mt-tRNA genes of the spider mite T. urticae. Reads mapping to all previously annotated T. urticae mt-tRNA gene regions were sequenced in all three AlkB-treated libraries and in two of three untreated libraries (only tRNA-Gln, tRNA-Leu2 and tRNA-Trp were not detected in the other untreated library; Supplementary Table S6). Similar to other tRNA expression profiles described utilizing YAMAT-seq (46,50), there was a high degree of variability in coverage, with the abundance of the most and least frequently sequenced tRNAs spanning more than three orders of magnitude (Figure 1 and Supplementary Table S6). Approximately 97% of reads mapping to a tRNA reference contained a 3′ CCA tail, which is typical of the mature, processed tRNAs that are targeted by YAMAT-seq (Supplementary Table S3).

Figure 1.

Mitochondrial genome map of Tetranychus urticae with tRNA and stem-loop expression. Outer circle shows the gene organization including expressed and CCA-tailed sequences determined by this study. Light green blocks are protein coding genes, dark green are ribosomal RNAs, black are stem-loops and dark blue are mt-tRNAs. The two gene tracks represent gene transcribed on the majority strand (outer track) and minority strand (inner track). The histograms show the sense (outward-facing, light blue bars) and antisense (inward-facing, dark gray bars) detection of reads mapping to each tRNA or stem-loop. Expression is in counts per million mapped reads on a log10 scale after summing all mapped reads in the six libraries and can be found in Supplementary Table S2a. Sense and antisense expression are defined relative to prior T. urticae mitochondrial annotations (13). Diagram was created using Circos (ver. 0.69–9) (80).

The majority of transcripts corresponded to the length and start/stop coordinates predicted from the prior annotation of the sequenced T. urticae mitochondrial genome (13) with some minor truncations compared to the genomic annotations; however, a few tRNAs appear to deviate considerably from the predicted models. tRNA-Phe and tRNA-Thr were previously annotated as having a D-arm, but all transcripts mapping to these genes were considerably shorter, possibly lacking both arms depending on how these structures fold in vivo (Figure 3 and Supplementary Table S5).

Figure 3.

Predicted tRNA secondary structures continued from Figure 2. See Figure 2 legend for details.

More surprisingly, the template strand of multiple tRNAs was found to be opposite to what has been previously annotated. Many RNA-seq methods do not retain the strand information of the original RNA template due to the synthesis and adapter ligation of randomly fragmented and primed double-stranded cDNA molecules (57). However, annealing of Y-shaped YAMAT-seq adapters directly to tRNA molecules in a defined orientation (as well as the presence of the CCA tail) provides 5′- and 3′-prime information, allowing for the identification of the transcribed strand (46). Strand information is particularly useful in tRNA-seq because the structure of tRNAs as well as the complementary nature of nucleic acid base-pairing means that the reverse complement of a tRNA gene can often produce a similar, ‘mirrored’ folding structure. Therefore, both strands could be predicted to produce a functional tRNA, albeit with different anticodons. Without experimental evidence and conservation of tRNA arrangement in a mitochondrial genome, inferring a tRNA’s identity and sense-orientation from genomic data alone can require a process of elimination searching for anticodons, and in cases of ‘reciprocal pairs’ of anticodons that are complementary to each other, some guesswork. For example, in the dust mite Dermatophagoides pteronyssinus, the same region of the mitochondrial genome has been annotated to encode either tRNA-Val or its antisense complement tRNA-Tyr depending on the study (12,19), demonstrating the challenge in assigning tRNA identities to genomic sequences in which complementary expression can also form tRNA-like structures. In the case of T. urticae, no sense expression of tRNA-Cys, tRNA-Val, or tRNA-Ala was detected (Figure 1). Note that tRNA-Cys and tRNA-Ala have complementary anticodons (GCA and TGC, respectively). Therefore, the finding of exclusive antisense expression for both of these two genes effectively ‘switches’ the anticodon/tRNA identity of the genes relative to their annotations in the genome. Similarly, the antisense-only expression of the tRNA-Val could be functionally compensated by the predominantly antisense expression of tRNA-Tyr (Figure 1, as these tRNAs also have complementary anticodons (GTA/TAC).

Expression and processing of non-tRNA stem-loops

Although the vast majority of transcripts mapped to a predicted tRNA gene, there was a detectable level of reads originating from specific positions within 12S rRNA, the putative control region (D-loop), a region internal to ND5 and an antisense region internal to CYTB, all with a post-transcriptionally added CCA-tail (Supplementary Table S3 and Figure 1). Post-transcriptional addition of CCA to non-tRNAs has now been described in the mitochondria of numerous other species, including rat (58), human (59) and plants (60). All of the apparent non-tRNAs that we detected in T. urticae have a predicted stem-loop structure (Supplementary Table S5), which may make them effective substrates for a CCA-adding enzyme.

The novel detection of sense–antisense expression of tRNAs

In addition to cases where only antisense reads were found (see above), mirrored expression of both sense and antisense transcripts was also detected for a number of T. urticae mt-tRNA genes. Specifically, tRNA-Arg, tRNA-Glu, tRNA-Leu1, tRNA-Phe, tRNA-Pro, tRNA-Ser2 and tRNA-Tyr all had both sense and antisense reads mapped but with variable percentages (Figure 1 and Supplementary Table S6). Very little sense expression of tRNA-Tyr was detected, with over 98% of all the reads mapped to the reference being in the antisense direction with an anticodon of TAC (Val). Conversely, tRNA-Glu and tRNA-Leu1 had only a single antisense read sequenced and thus lack evidence for a meaningful level of antisense expression. The remaining tRNAs with both sense and antisense expression all had at least seven (and up to 348) antisense reads, representing anywhere from <1 to 43% of the mapped sequences the respective reference gene.

To assess whether the mirrored sense–antisense expression observed for T. urticae mt-tRNAs was a general phenomenon or perhaps an artefact of the tRNA-seq method, we examined the abundant mt-tRNA reads generated from the intended host plant species, S. vulgaris. We found that a total of five (imperfectly matched) reads were generated from a single plant mitochondrial gene (mt-tRNA-Trp) out of 105 867 reads that were mapped to that tRNA (representing 0.00005%). There were no antisense reads and 1 215 689 sense reads for the remaining four S. vulgaris mt-tRNA reference genes. We also analyzed reads derived from S. vulgaris plastids and found that there was a total of just 39 reads (or 0.000014% out of 2 532 488 total mapped for the 30 plastid tRNA genes) that could have an antisense orientation, and many of these reads failed a hit coverage threshold in mapping (i.e. fragments). Thus, the abundant sense–antisense expression appears to be a distinct feature of T. urticae mt-tRNAs and not a universal outcome of tRNA-seq in these heterogeneous samples. However, caution should be applied when interpreting the relative expression between sense and antisense transcripts from the same gene, as substantial reverse-transcription biases still exist in tRNA-seq datasets (46,50). Similarly, the lack of detection of a tRNA transcripts, in either orientation, should not be taken as definitive evidence of lack of expression.

Evidence for base modifications affecting reverse transcription at positions 8 and 9 and 3′-end processing of overlapping tRNAs

Misincorporation of incorrect nucleotides during the reversion transcription of tRNAs has been used to predict the modification state of tRNA bases (40,41). We considered a base confidently modified if it was misread in >30% of the reads from untreated libraries (Figures 2 and 3). Additional positions had misincorporation rates below the 30% cutoff and are listed in Supplementary Table S7. No position had a 100% misincorporation rate, but position percentages ranged up to 99% misread in all libraries (position 9T in antisense tRNA-Cys; Supplementary Table S7), to never being misread. Thymines at positions 8 and 9 and adenines at position 9 in T. urticae mt-tRNAs were the most likely to exhibit differences relative to the reference sequence (Figures 2 and 3; Supplementary Table S7 and Figure S1). The percentage of reads with the correct base at position 9A greatly increased in at least two of the three AlkB-treated libraries for tRNA-Asn, tRNA-His, tRNA-Leu2 and tRNA-Met (Supplementary Table S7 and Figure S1). AlkB treatment did not appear to be as effective in one of the replicates (replicate 2). An apparent error in quantifying RNA concentration during AlkB treatment of replicate 2 resulted in a higher ratio of RNA to AlkB and may account for some replicate behavior. The misincorporation signal was often almost eliminated (i.e. the correct base was read) after AlkB treatment in replicates 1 and 3 for position 9A in tRNA-Asn, tRNA-His, tRNA-Leu2 and tRNA-Met, whereas it was only reduced (but not eliminated) in replicate 2 (Supplementary Table S7 and Figure S1). The most frequently misread base was a 9T, and the number of reads with the correct base did not increase after AlkB treatment (i.e. an apparently AlkB-insensitive modification) (Supplementary Figure S1 and Table S7).

Figure 2.

Predicted secondary structures of Tetranychus urticae mt-tRNAs and stem-loops. Most structures presented were based on Chen et al. (13) and tRNAscan-SE (28). The primary reference sequence is based on the most highly detected YAMAT read that had 100% base identify to the reference genomic sequence. Bases were inferred to be modified if >30% of the reads at a position differed from the reference sequence in untreated libraries and are highlighted with a dashed red circle. Bases detected in sequenced reads that were not present in the genomic reference are indicated with an asterisk. Some tRNAs had frequently occurring alternative 3′- and 5′-ends which are presented next to the structure and indicated with a ‡. Folded structure diagrams were created using VARNA ver. 3.9 (81). Although YAMAT-seq provides information about the location of tRNA transcript ends, including the position of the acceptor stem, it does not assess the overall secondary structure of tRNAs, so these structures are speculative predictions.

Some transcripts also contained nucleotides (adenosines and cytosines in addition to the 3′ CCA tail) not encoded in the mitochondrial genome, including sense reads from tRNA-Gln and antisense reads from tRNA-Ser2 (Figures 2 and 3), which suggests that these nucleotides may be post-transcriptionally added. Prior tRNA gene annotations for tRNA-Gln had predicted a mispaired 4A:43C and overlap with tRNA-Leu2 (13). However, sequenced reads suggest that transcription of genomically encoded sequence ends at position 38 (Supplementary Table S3). Similarly, the antisense expression of tRNA-Ser2 appears to end at position 36 and is finished with the post-transcriptional addition of nucleotides, thereby preventing overlap with the antisense expression of tRNA-Ala (Figure 3 and Supplementary Table S3).

DISCUSSION

The predicted sequence and structure of aberrant mt-tRNAs in multiple metazoan systems has raised questions about their expression, maturation and functionality because of short sequence length, base mispairing in stems, and a lack of canonical tRNA structure (5,9,61–63). These extreme deviations from canonical tRNA structure have made identifying mt-tRNA with genomic information alone difficult, leading to some studies inferring their outright loss (14,18,21–22). Our accidental sequencing of T. urticae mt-tRNAs, an arachnid from a superorder (Acariformes) well known for degenerate mt-tRNAs, has confirmed that previously predicted tRNA genes, some as short as 44 nt at maturity, are transcribed and modified with the post-transcriptional addition of the CCA-tail (Supplementary Table S3). There is also evidence from sequencing misincorporations at specific sites (especially positions 8 and 9) that these tRNAs have post-transcriptionally modified bases. One position, 9A, was frequently misread in T-armless mt-tRNAs and appeared to be removed by AlkB treatment (Supplementary Figure S1 and Table S7). One possibility is that this site carries a 1-methyladenosine (m1A9) modification, which is known to be removed by AlkB (40,43–44) and has been found to be necessary for efficient aminoacylation and EF-Tu-binding of T-armless mt-tRNAs in nematodes (64). Misincorporations at positions 8T and 9T were most common, and we hypothesize that these positions are modified to prevent base pairing with A/T rich replacement loops, thereby maintaining the L or boomerang-like structure necessary for ribosome interaction (11,65). Some modifications in arms or replacement loops may even be preventing the formation of stem structure, thereby maintaining an armless state. We advise that any specific sites of interest should be investigated for possible contributions of sequence variation at the genomic level to apparent misincorporation patterns. Overall, the expression, end-processing, and base modification of T. urticae mt-tRNAs suggest that the mt-tRNA gene complement is functional and sufficient to decode the mitochondrial genome despite extremely degenerate structures.

However, caution is warranted when inferring a role in translation because stem-loop structures (e.g. in the D-loop and certain non-tRNA genes) were also expressed and processed, and it is unlikely that these function in translation as tRNAs. Although the post-transcriptional modification of tRNAs and especially the addition of the CCA-tail have been frequently used as evidence for functionality of tRNAs (10,11), the detected processing of apparent non-tRNA stem-loops in T. urticae and other systems (58–60) may mean that the addition of a CCA-tail alone is insufficient to assign function as a tRNA in translation. Aminoacylation and EF-Tu binding assays will be necessary to determine whether tRNAs and stem-loop structures are functional in translation or have other functions in the organelle.

Mt-tRNAs have long been hypothesized to have additional functional roles in both vertebrate and invertebrate mitochondrial expression based on co-opting tRNA-interacting enzymes (RNase P and Z) for mRNA and rRNA transcript maturation through the endonucleolytic cleavage of polycistronic RNAs, releasing transcripts immediately before and after a tRNA sequence (33,66). The sense and antisense expression of mt-tRNA-Tyr and tRNA-Arg in T. urticae occur at protein-coding boundaries on both strands, raising the possibility that tRNA genes can act in a novel form of dual-strand punctuation by being recognized by enzymes in both orientations. As such, the conservation and expression of tRNAs or tRNA-like structures in acariform mitochondria could be for the maintenance of a structure generally recognizable by enzymatic partners in transcript maturation pathways for other genes (19), as well as possible replication initiation sites (35,36). The loss of translation-related function of a mt-tRNA but the conservation of some tRNA structural features (i.e. the acceptor stem) for transcript maturation has been suggested for marsupial mt-tRNA-Lys, where a nonfunctional tRNA-Lys is retained in the mitochondrial genome as a processing signal despite its functional replacement by an imported nuclear-encoded tRNA-Lys (67). However, the antisense expression and modification of multiple mt-tRNAs in T. urticae not immediately up or downstream of any previously annotated genes suggests that expression may be maintained for functions other than as a processing signal. Alternatively, some processed sequences may simply be a byproduct of being in-between two sense expressed tRNAs. In addition, there appears to be extensive modification of some tRNAs, including potential post-transcriptional 3′ terminus ‘finishing’ through the addition of nucleotides other than the CCA-tail (Figures 2 and 3). In the case of tRNA-Gln, these nucleotide additions appear to restore a correctly base paired acceptor stem when a base mispairing exists in the genomic sequence. Taken together, the extensive 3′-terminus modification and possibly stabilizing chemical base modifications (i.e. positions 8 and 9) do suggest the conservation of pathways to maintain structure for enzymatic and ribosome interaction.

If all or most of the T. urticae mt-tRNAs are functional in translation, the structural variability of the mt-tRNAs pool raises additional questions about the enzymatic partners necessary for mt-tRNA function. Enzymes and complexes that interact with tRNAs (e.g. EF-Tu and the ribosome) often use conserved structural features such as the T-arms for substrate recognition and interaction (68,69). The coexistence of tRNAs lacking either or both arms in addition to tRNAs with conventional shape (i.e. both arms present) results in T. urticae having a mt-tRNA pool with extreme structural variability. A similar situation in nematodes has arisen where mt-tRNAs either lack a D- or T-arm such that tRNA interactions require two distinct EF-Tu enzymes, one binding only to T-armless tRNAs and the other binding to D-armless Ser-tRNAs (70,71). What coevolutionary features of nuclear-encoded enzymes facilitate both degenerate as well as variable mt-tRNA structure is Acariformes is yet to be investigated. Given that highly degenerate mt-tRNAs have arisen independently in multiple lineages, the evolution of the tRNA-interacting enzymes may offer clues into the propensity for the evolution divergent mt-tRNAs (61,69).

The surprising detection of sense–antisense expression of multiple tRNAs theoretically expands the function of T. urticae mt-tRNAs as decoding molecules in mitochondrial protein synthesis. The possibility of a single mt-tRNA gene serving as two different isoacceptors has been suggested for the nematode Romanomermis iyengari where tRNA-Ala and tRNA-Cys are predicted as superimposed reverse complements at the same genomic location (59). There have been some reports providing limited evidence of antisense mt-tRNA expression in mammals and yeast (72,73), and long noncoding RNAs (lncRNAs) in human mitochondria containing both antisense protein coding and tRNAs have been found to be substrates for RNAse P, where a mirrored tRNA sequence serves as a cleavage point for lncRNA release (74). However, dual sense–antisense expression and post transcriptional processing of mt-tRNA genes has remained a largely undetected phenomenon despite the obvious potential for tRNAs to produce mirrored secondary structures (75).

In the case of human mitochondrial lncRNAs, RNA–RNA duplexes were found to form between mRNAs and the corresponding antisense components of the lncRNAs, possibly regulating mRNA stability (74). Antisense expression of RNAs has been extensively reported in bacterial species, where the majority act by base pairing with complementary mRNAs, likely having a modulating effect on translation (76). While it is possible that the antisense mt-tRNAs in Tetranychus hybridize with their complementary tRNA counterparts, it is unlikely that the molecules sequenced in this study were extracted in that double-stranded form because the Y-shaped adapters used in YAMAT-seq require a single contiguous RNA strand with a conventional acceptor stem to produce a sequenceable library molecule. The rapid degradation of most antisense transcripts in human mitochondria suggest that the removal of nonfunctional antisense RNA is physiologically important (77). What, if any, regulatory role these mirrored tRNAs have on sense tRNA expression is yet to be determined.

If these mirrored tRNAs are functional as decoding molecules, the transcription of two tRNAs from a single tRNA gene could have consequences for both mitochondrial protein expression and evolution. One obvious outcome from dual tRNA expression could be the facilitation of mitochondrial genome reduction through the tolerance of mt-tRNA gene loss. Interestingly, sense–antisense expression of mt-tRNAs has been proposed to account for the highly pathogenic nature of mutations in mt-tRNA genes (75). Disease causing mutations are 6.5 times more frequent in mt-tRNAs than in other mitochondrial genes and this discrepancy in pathology may arise from mutations perturbing both the sense and antisense function of a tRNA as decoding molecules in translation (78). The complementary expression of tRNAs has even been suggested to represent the ancestral decoding system of life (79). As such, our findings in support of the long-hypothesized potential for mirrored sense–antisense expression of tRNAs has broad implications for the function and evolution of mitochondrial genomes.

Click here for additional data file.