Anticodon

The genome assemblies used in this study were scanned with tRNAscan-SE v1.3 and 2.0 using the corresponding domain search modes (-E, -B or -A options). When running v1.3 on archaeal genomes, the - -ncintron option was also included to enable noncanonical intron searches (this functionality is default in 2.0). Predicted tRNA genes in archaeal, bacterial and fungal genomes were compared on the basis of similarity of genomic coordinates and tRNA identity (isotype and anticodon). Results were grouped into four categories: (i) consistent—the predicted gene has consistent identity and start and/or end positions differ by 10 nucleotides or less, (ii) isotype mismatch—the predicted gene coordinates are the same but the isotype does not match, (iii) novel—the gene is only predicted by v2.0 but not v1.3, (iv) not detected—the gene is predicted by v1.3 but not v2.0. Within the ‘consistent’ category (Table 2), only 0.69% archaeal, 6.8% bacterial, and 1.7% fungal tRNA predictions were found to have slightly different start and/or end positions. Model organisms (Table 3) were compared by predicted gene counts and program execution time. All prediction runs were conducted on Linux servers with identical configurations (dual Intel 10-core HT processors at 2.30GHz and 128 GB memory). tRNAscan-SE 2.0 utilized eight parallel threads for Infernal (20 (link)) searches (- -thread 8 option) while v1.3 only allows a single thread.

We assembled three sets of domain-specific genomic tRNA sequences from a total of over 4000 genomes using existing tRNAscan-SE 1.3 predictions in GtRNAdb Release 15 (4 (link)) plus additional predictions (also from tRNAscan-SE 1.3 for species not yet represented in GtRNAdb), representing a broad diversity of eukaryotes, bacteria, and archaea (Table 1 and Supplementary Figure S1). Before using these tRNA sequences as training sets for building domain-specific covariance models, multiple filtering steps were used to maximize quality. To avoid the inclusion of common tRNA-derived repetitive elements that exist in many eukaryotic genomes, especially those in mammals (33–35 (link)), we first selected only the eukaryotic tRNAs with a COVE score >50 bits, a threshold reflecting more conserved, canonical tRNA features. We then selected only the top 50 scoring tRNAs for each isotype per organism to avoid overrepresentation of high-scoring tRNA-derived repetitive elements which are abundant in some species (e.g. elephant shark has over 9500 tRNA^Ala scoring over 50 bits). For the bacterial tRNA training set, all genes having potential self-splicing introns were excluded to eliminate large alignment gaps (i.e. mostly over 200 nt) which can hinder efficient model creation. Similarly for archaea, pre-processing of sequence training sets was necessary. Some species within the phyla Crenarchaeota and Thaumarchaeota contain tRNA genes that are known to have multiple noncanonical introns (5 (link),36 (link),37 (link)). Atypical tRNAs such as trans-spliced tRNAs and circularized permuted tRNAs have also previously been discovered in Crenarchaeota and Nanoarchaeota (12 (link),38 (link),39 (link)). To accommodate these special archaeal features without sacrificing performance, both mature tRNA sequences (without introns) and selected atypical genes with multiple introns at different locations were included in the archaeal tRNA training sets. As a last step, anticodons of tRNA sequences in all training sets were replaced with NNN and aligned to the corresponding original domain-specific tRNA covariance models using Infernal (20 (link)). The resulting alignments were then used to generate the new set of domain-specific tRNA covariance models with Infernal.

EukHighConfidenceFilter in tRNAscan-SE 2.0 is a post-scan filtering program for better distinguishing tRNA-derived repetitive elements in large eukaryotic genomes (metazoans and plants) from ‘real’ tRNAs that function in protein translation. Three filtering stages are involved in the classification (Supplementary Figure S5). First, tRNA predictions that are labeled as possible pseudogenes are excluded from the high confidence set; these criteria were established by the prior version of tRNAscan-SE (overall score below 55 bits and one of two conditions: primary sequence score below 10 bits or secondary structure score below 5 bits). Second, predictions with any of the following attributes are removed from the high confidence set: isotype-specific model score below 70 bits, overall score below 50 bits, or secondary structure score below 10 bits. Finally, if there are >40 predicted hits remaining for any given anticodon, a dynamic score threshold is used, starting at 71 bits, rising one bit and filtering lower-scoring hits, iteratively, until the number of predictions for that anticodon is no longer over 40 or when the score threshold reaches 95 bits. The score thresholds in stages two and three were empirically determined by comparing score distributions of predictions among eukaryotic genomes with and without large numbers of tRNA-derived repetitive elements (Supplementary Figure S6). The remaining tRNA predictions are included in the high confidence set if (i) they have a consistent isotype prediction (inferred from anticodon versus the highest scoring isotype-specific model) and (ii) they have an ‘expected’ anticodon based on known decoding strategies in eukaryotes where 15 anticodons are not used (67 (link)). There is no corresponding filter for bacterial or archaeal tRNA, as they have not been found to contain large families of tRNA-derived repetitive elements.

The total DNA was isolated from a single specimen and used as the PCR template. Primers were designed to produce amplicons overlapping by approximately 100 bp. PCR products were sequenced bidirectionally by the Bio-Transduction Lab in Wuhan using the Sanger method and the same set of primers that were used for the amplification (Additional file 1: Table S1). Sequenced data were quality-inspected using sequencing chromatograms; 20–30 bases at both ends were deleted from each amplicon. The mitogenome was assembled manually using DNAstar v7.1 [55 ], ensuring that overlaps between fragments were identical. PCGs and rRNAs were approximately located using DNAstar and MITOS [56 (link)] and then manually fine-tuned according to the orthologous sequences using BLAST and BLASTx [57 (link)] and PhyloSuite [58 , 59 (link)]. tRNAs were annotated using MITOS [56 (link)] and ARWEN [60 (link)] tools, but some tRNAs had to be identified manually via comparisons with related species. The secondary structure of tRNAs was further studied using ARWEN. Tools such as tRNAscan, which searches for a complete cloverleaf structure, are not suitable for the identification of tRNAs that exhibit nonstandard secondary structures. ARWEN was designed especially for this purpose: it first identifies only the most conserved domain, the anticodon stem, and subsequently searches for the presence of D-stem and T-stem structures, and the search for an acceptor stem then provides specificity. Due to this high sensitivity, ARWEN is also prone to producing a substantial false discovery rate [47 (link)]. The assembled circular mitogenome was visualised using OGDRAW [61 (link)].

tRNA-derived fragment and tiRNA sequencing was performed as previously described (25 (link)). Briefly, total RNA was extracted using TRIzol (Invitrogen) from three pairs of HASMCs with or without PDGF-BB (50 ng/ml) treatment for 24 h and quantified by a NanoDrop (Thermo Fisher Scientific, USA). tRNA-derived fragments (tRF and tiRNA) are heavily modified by RNA modifications that interfere with small RNA-seq library construction. RNA modifications need to be removed before library construction, which is the biggest difference from common small RNA sequencing. RNA modifications were removed by the rtStar™ tRF and tiRNA Pretreatment Kit (Arraystar, Rockville, MD, USA). A sequencing library specific for tRFs was constructed from the pretreated total RNA by the following steps: (1) 3′ and 5′-adapter ligations and m1A and m3C demethylation; (2) cDNA synthesis; (3) PCR amplification; and (4) to differentiate tRFs with full-length tRNA, sequencing libraries are size-selected for the RNA biotypes to be sequenced using an automated gel cutter. Size of 134–160 bp PCR-amplified fragments (corresponding to a 14–40 nt small RNA size range) were selected, while the full length tRNA is 76–90 nt. The libraries were denatured into single-stranded DNA molecules, captured on Illumina flow cells, amplified in situ as sequencing clusters and sequenced for 50 cycles on an Illumina NextSeq 500 system (Illumina, CA, USA) using a NeXTs 500/550 V2 kit (#FC-404-2005, Illumina, CA, USA). The sequencing data has been uploaded to the GEO database (GEO accession number GSE212737).
The abundance of tRF and tiRNA was evaluated using their sequencing counts and was normalized as counts per million of total aligned reads (CPM). Differentially expressed tRFs and tiRNAs analyses was performed with R package edgeR (33 (link)). Fold change (cut-off 1.5), p-value (cut-off 0.05) were used for screening differentially expressed tRFs and tiRNAs. The number of subtypes tRFs, which the CPM of the sample or the average CPM of the group is not less than 20, can be counted against tRNA isodecoders which share the same anticodon but have differences in their body sequence. The stacked bar chart is plotted with R barplot package.