Transfer RNAs (tRNAs) are short, noncoding molecules acting as intermediaries between the genetic information in nucleic acids and protein sequences. Although the mechanistic roles of tRNAs in ribosomal protein biosynthesis are well understood, they have noncanonical functions in several aspects of cellular metabolism. Plant tRNAs have been implicated in tetrapyrrole and cytokinin biosynthesis (Chery and Drouard, 2022), plant cell growth and immunity (Soprano et al., 2018), and regulation of auxin response in Arabidopsis (Leitner et al., 2015). Increased attention has also been given to tRNA-derived fragments (tRFs), a class of small RNAs produced from the enzymatic cleavage of tRNAs. Initially thought as mere tRNA degradation byproducts, tRFs have been linked to gene regulation, ribosome biogenesis, plant–pathogen interactions, and stress response in plants (Park and Kim, 2018; Alves and Nogueira, 2021; Wang et al., 2023; Panstruga and Spanu, 2024). tRFs have been reported to be involved in the RNA silencing pathway and are the major source of the transposable element SINEs (short interspersed nuclear elements; Bermudez-Santana et al., 2010; Phizicky and Hopper, 2010; Soprano et al., 2018). All tRNA genes are postulated to be derived from an ancestral “proto-tRNA” (Eigen et al., 1989), and during evolution, a tRNA repertoire was generated from gene duplication and numerous mutational events. These processes gave rise to the core and dispensable sets of tRNA genes.

Despite the growing knowledge and interest in plant tRNA biology, studies on how tRNAs are structured and organized on a genome-wide scale in plants still number too few. A survey of the content, distribution, and clustering of tRNA genes and pseudogenes in many eukaryotes, including nine genomes from the green lineage, has been reported (Bermudez-Santana et al., 2010). More recent studies have also reported the evolution of tRNA gene content in the three domains of life, involving 13 plant genomes (Santos and Del-Bem, 2023), as well as the tRNA anticodon frequency of 128 plant genomes (Mohanta et al., 2020). Databases of tRNA gene sets from hundreds of plant nuclear and organellar genomes, covering diverse families of plants, have also been developed (e.g., Cognat et al., 2022; Mokhtar and Allali, 2022), whose curators were able to provide a general survey of the tRNA gene populations of 51 and 256 plant species, respectively. However, these mostly only provided insights on the tRNA gene content of these plants, and separate studies fully utilizing the information from these databases are yet to be found. To date, the first comprehensive study that focused on tRNA gene content, structure, and distribution in plants covered both the nuclear and organellar genomes of only five angiosperms—consisting of three eudicots and two monocots—and one green alga (Michaud et al., 2011). However, given the species diversity within the flowering plants, a more comprehensive and systematic comparative study is needed to provide a global landscape of plant tRNA structure and organization. The increased availability of plant genomes will provide common patterns and taxon-specific particularities of plant tRNAs.

Compared to other eukaryotic genomes, plant genomes possess a smaller variation in the number of tRNA genes and a varying abundance of tRNA gene clusters (Bermudez-Santana et al., 2010). The following tRNA gene organization has been reported among flowering plant genomes: a predominantly A-/T-rich region spanning 50 nucleotides upstream of the tRNA gene, an upstream CAA motif and a downstream poly(T) termination signal found in most tRNAs, and intron-containing tRNA^Met and tRNA^Tyr genes (Michaud et al., 2011). Except for Arabidopsis, a similar chromosomal distribution of tRNA genes in terms of the numbers of tRNA genes per megabase of the chromosome was also reported within angiosperms, which hinted at the possibility of excessive tRNA gene duplications in some plant genomes (Michaud et al., 2011). Although a significant correlation between genome size and number of tRNA genes have been reported among 74 eukaryotic genomes (Bermudez-Santana et al., 2010), five plant genomes (Michaud et al., 2011), and eight monocot genomes (Planta et al., 2022), a more recent regression analysis involving a higher coverage of plant genomes (128 genomes) instead reported a weak correlation (Mohanta et al., 2020).

In the case of organellar genomes, previous studies also reported the lack of certain tRNA isoacceptors in some plant plastomes and mitogenomes (Michaud et al., 2011; Mohanta et al., 2020). Although possessing significantly fewer tRNA genes than the nuclear genome, the organellar genomes from mitochondria and chloroplasts can also encode their tRNAs. The chloroplast genome is assumed to encode all the tRNA species required for protein synthesis, and unlike the mitochondria, chloroplasts do not import cytosolic tRNAs (Maréchal-Drouard et al., 1993). A relaxed wobble rule might also explain the small number of organellar tRNAs that can read all codons of the universal genetic code (Crick, 1966; Percudani, 2001).

Several different sequencing-based approaches have been developed to quantify highly modified tRNAs. However, modifications on tRNAs can impair cDNA synthesis by premature reverse transcriptase (RT) stops (Pinkard et al., 2020; Padhiar et al., 2024). These methods [e.g., ARM-seq (Cozen et al., 2015), DM-TGIRT-seq (Zheng et al., 2015), YAMAT-seq (Shigematsu et al., 2017), Nano-tRNAseq (Lucas et al., 2024); see Padhiar et al. (2024) for a comprehensive review] incorporate pre-treatment of RNA before library construction and the use of modified adapters; pre-treatment of RNA produces less complex secondary structures and fewer modifications that can lead to premature RT stops (Padhiar et al., 2024). Plant tRNA expression and post-transcriptional modifications have been characterized in Arabidopsis thaliana by modifying RNA-seq to involve a demethylating enzyme and using a tRNA-specific adapter (Shigematsu et al., 2017; Warren et al., 2021). While these are promising advancements in direct tRNA sequencing, at its current state, computationally predicting tRNA genes from whole-genome sequencing data is still the preferred method in most tRNA gene studies (Chan et al., 2021).

This study compared and analyzed the tRNA gene content, structure, and organization of 69 nuclear plant genomes—including available chloroplast and mitochondrial genome counterparts ( Supplementary Table 1 ). Included in our analyses are 44 genomes from the eudicot lineage, 20 from the monocot lineage, four from the ANA clade (Amborellales, Nymphaeales, and Austrobaileyales), and one from Ceratophyllaceae, the sister clade to eudicots. The eudicot and monocot genomes were chosen to cover as much family in the flowering plant lineage; the chosen plant genomes span 32 families—two from the ANA clade (Amborellaceae and Nymphaeaceae), Ceratophyllaceae, nine from monocots, and 20 from dicots ( Supplementary Figure 1 ). Having these lineages within the scope of this study should provide a better and more inclusive analysis of tRNA genes in plants. Using the widely adopted tool tRNAscan-SE (Chan et al., 2021), tRNA genes from these genomes were computationally predicted and then filtered for a “high-confidence” set of tRNA genes discarding pseudogenes. To characterize these “high-confidence” tRNA genes, we also screened the tDNAs for regulatory sequences commonly associated with the RNA polymerase III-transcribed plant tRNA genes: the upstream TATA-box and CAA motifs (Choisne et al., 1998; Yukawa et al., 2000; Dieci et al., 2006; Michaud et al., 2011), the intragenic A and B boxes (Choisne et al., 1998; Dieci et al., 2006), and the downstream poly(T) stretches (Yukawa et al., 2000; Braglia et al., 2005; Arimbasseri and Maraia, 2015).

Comparative genomics analyses revealed that the number of nuclear tRNA genes varied mainly among the plant genomes studied, even among genomes of the same lineage. In contrast, the number of organellar tRNA genes had slight variation and was consistent regardless of plant lineage. Moreover, gene duplications in tRNA gene clusters appeared more prevalent in eudicots. All nuclear genomes were found to have a strictly conserved tRNA^Ala-AGC species and intron-containing tRNA^eMet and tRNA^Tyr genes that exhibited modest sequence conservation. Regulatory sequences found in the nuclear tRNA genes include the upstream TATA-box and CAA motifs (found upstream of 22%–32% and 78%–82% of tRNA genes detected, respectively), the intragenic A and B boxes (found in all tRNA genes detected) with general lineage- and isotype-specific motifs, and the downstream poly(T) termination signals (found downstream of 67%–72% of tRNA genes detected). Overall, this study revealed differences and patterns acquired by plant genomes throughout evolution and can serve as a foundation for further studies on plant tRNA gene function and regulation.

Plants with sequenced chloroplast, mitochondrial, or nuclear genomes ( Supplementary Table 1 ) were used for the comparative analysis of tRNA gene content, structure, and organization. Aquilegia coerulea and Acorus americanus were included in the analysis as these are members of the basal-most eudicot clade and the sister lineage to all other monocots (Filiault et al., 2018; Givnish et al., 2018), respectively. Amborella trichopoda, Nymphaea colorata, Nymphaea thermarum, and Euryale ferox under the ANA (Amborellales, Nymphaeales, and Austrobaileyales) clade are sisters to all other angiosperms. Ceratophyllum demersum belongs to the species-poor lineage of Ceratophyllales and is sister to eudicots (Yang et al., 2020). Given the phylogenetic positions of these species ( Supplementary Figure 1 ), including these sequences will facilitate better comparative analysis of the tRNA gene arrangement and structure in flowering plants.

A narrow range of nuclear tDNA numbers in angiosperms (500–600 tDNAs between five angiosperm genomes) had been previously reported (Michaud et al., 2011), and extending the coverage to 69 angiosperm genomes resulted in a broader range in the number of nuclear tDNAs that were detected (approximately 150–1,500 tDNAs; Figure 1A ). This tDNA range is comparable to that reported by Bermudez-Santana et al. (2010), although their range also included tRNA pseudogenes (432–1,290 tDNAs between seven land plant genomes). In addition, the green algae Volvox carteri and Chlamydomonas reinhardtii were reported to have 1,051 (including tRNA pseudogenes) and 256 tDNAs, respectively (Bermudez-Santana et al., 2010; Michaud et al., 2011). Therefore, nuclear genomes from the green lineage can have tDNAs as few as 150 or as many as 1,500. This variation in plant nuclear tDNA numbers is relatively small compared to other eukaryotes. Tetraodontiformes have approximately 700 tDNAs, while the related zebrafish, Danio rerio, has approximately 20,000 (Bermudez-Santana et al., 2010). Concurrently, in mammals, old-world monkeys and apes had 496–736 tDNAs, while cows and rats exceeded 100,000 tDNAs (Bermudez-Santana et al., 2010); a reannotation of the cow tRNAs showed that the majority of these putative tDNAs include tRNA-like sequences (Theologis et al., 2000; Tang et al., 2009). In nuclear eukaryotic genomes, the number of tDNAs can vary even within species of the same lineage or clade. Indeed, ANA, eudicot, and monocot genomes have varying numbers of nuclear tDNAs, and no lineage-specific pattern could be observed ( Figure 1A ).

The varying genome sizes in eukaryotes, including plants, could explain this variation in the number of tDNAs. While earlier studies suggested a strong correlation among plants, with Arabidopsis being an outlier (Bermudez-Santana et al., 2010; Michaud et al., 2011), our data showed a weak overall correlation in the 69 angiosperm genomes studied (R² = 0.41, p-value <0.0001; Figure 1C ), especially for the eudicot lineage, with an R-squared value of 0.29 (p-value = 0.0002). More recent studies have similarly reported a weak correlation among plants (Mohanta et al., 2020; Santos and Del-Bem, 2023). However, this was not the case for the monocot lineage, which exhibited a strong correlation (R² = 0.79, p-value <0.0001; Figure 1C ). A strong correlation between the monocot genome sizes and the number of tDNAs had been previously reported (Planta et al., 2022).

At least for the eudicot genomes, a likely explanation is related to the unique case of Arabidopsis (Michaud et al., 2011). A weak correlation between the number of tDNAs and genome size was initially shown in A. thaliana, with an R-squared value of 0.16. This correlation contrasted with the other analyzed plant genomes, which all had moderate to high R-squared values. Compared to four other angiosperms (Medicago truncatula, Populus trichocarpa, Oryza sativa, and Brachypodium distachyon) and one green alga (C. reinhardtii), A. thaliana had a higher number of tDNAs in each chromosome (Michaud et al., 2011). Except for A. thaliana, the other genomes had at most only two tDNAs per Mb of chromosome. Chromosomes 2–5 of A. thaliana had approximately four tDNAs per Mb, while Chromosome 1 had eight tDNAs per Mb of chromosome (Michaud et al., 2011). This unusually high number of tDNAs in Chromosome 1 of A. thaliana is largely due to the existence of two large tDNA clusters in this chromosome: tandem repeats of 27 tRNA^Pro and tandem repeats of 27 tRNA^Tyr–tRNA^Tyr–tRNA^Ser (Theologis et al., 2000). These clusters, indicative of gene duplications (Theologis et al., 2000; Bermudez-Santana et al., 2010), are likely the cause of the weak correlation between the tDNA number and genome size of A. thaliana. Indeed, removing the tRNA isotypes involved in the two identified clusters (tRNA^Pro, tRNA^Ser, and tRNA^Tyr) increased the R-squared value in A. thaliana from 0.16 to 0.70 (Michaud et al., 2011).

Similarly, the weak overall correlation found in the angiosperm genomes in this study might be explained by the prevalence of gene duplication events. This is likely the case, given that generally less than half of all tDNAs of each lineage were found to be unique ( Figure 3 ). This may also explain the observation that plants, alongside vertebrates, appear to have higher tDNA count and redundancy compared to other organisms (Santos and Del-Bem, 2023). However, this does not explain why the monocots showed a strong correlation between tDNA number and genome size (R² = 0.79; Figure 1C ), as opposed to weaker correlation observed in eudicot genomes (R² = 0.29). The key difference may lie in the existence of tDNA clusters, like the ones found in A. thaliana.

We considered tRNA genes to be clustered if at least three tDNAs were within 1 kb of each other. Using this criterion, 324 (in 40/44 genomes), 103 (in 13/20 genomes), and 2 (in 1/4 genomes) tDNA clusters were identified in eudicots, monocots, and ANA, respectively. Eudicots thus appear to have a stronger tendency toward gene duplication in the form of tDNA clustering compared to the other plant lineages, and this should explain the weaker correlation between tDNA numbers and genome sizes in eudicots compared to those in monocots. While ANA genomes appear to have a weak correlation like eudicots ( Figure 1C ), they had very few tDNA clusters. It is very likely that the linear regression model did not properly represent the correlation between ANA genome size and tRNA gene count due to the high p-value (0.7677). This may also be a result of our stricter criteria for tDNA clustering compared to other tRNA studies (Bermudez-Santana et al., 2010; Morgado and Vicente, 2019), which considered clusters as having at least two tDNAs within 1 kb of each other.

We identified tDNA clusters in Chromosome 1 of A. thaliana, similar to the two large clusters that were previously reported (Michaud et al., 2011) as follows: (i) consecutive tRNA^Pro clusters, adding up to 25 tandem repeats of tRNA^Pro, and (ii) consecutive tRNA^Tyr-tRNA^Ser clusters, comprising a long stretch of alternating tRNA^Tyr and tRNA^Ser genes. Unlike previously reported, these stretches of tRNA^Tyr and tRNA^Ser genes were not strictly tandem repeats of the triplet tRNA^Tyr–tRNA^Tyr–tRNA^Ser. The difference in the size and order of these clusters compared to those found by Theologis et al. (2000) is likely due to the updated genome assembly for A. thaliana. These tRNA^Pro and tRNA^Tyr–tRNA^Ser clusters were also found in other plant genomes. Most eudicots (34 out of 44 genomes, including A. thaliana), a few monocots (6 out of 20 genomes), and C. demersum were also found to have stretches of tRNA^Pro genes. On the other hand, a long stretch of alternating tRNA^Tyr and tRNA^Ser genes was also found in eight other eudicot genomes (Boechera stricta, Diptychocarpus strictus, Iberis amara, I. tinctoria, Lunaria annua, Lepidium sativum, Malcolmia maritima, and Myagrum perfoliatum). This tRNA^Tyr–tRNA^Ser tDNA cluster was not found in any other monocot or ANA genome. Another tDNA cluster detected is a tandem repeat of 28 tRNA^Ile found exclusively in Chromosome 2 of Z. mays. Among the clusters found in this study, this is the longest in size. Interestingly, this cluster is followed by three more clusters consisting purely of tRNA^Ile (5x tRNA^Ile, 3x tRNA^Ile, then 4x tRNA^Ile) within the same chromosome.

It remains to be seen whether these tDNA clusters serve any biological purpose. tDNA clusters are implicated in genome breakage resulting in genome rearrangement (Rienzi et al., 2009). They are also found to be involved in mobile genetic elements and horizontal gene transfer (Morgado and Vicente, 2019). tDNA clusters are likely dynamic and fragile genomic regions, and this inherent instability might be the reason for the evolution and prevalence of these tDNA clusters rather than being products of positive selection. Moreover, a study on the tDNA clusters of Arabidopsis shows that these clusters are predominantly methylated and transcriptionally repressed (Hummel et al., 2020). However, the case of tRNA^Pro clusters is intriguing given its frequency among the plant genomes studied.

Proline is found to have diverse roles in plants. They are involved in cell wall and plant growth (Kishor et al., 2015), but the more well-documented function of proline is related to plant stress. In response to different environmental stresses, e.g., drought or water loss, salt, metal, and pathogen attack, plants accumulate proline (Kishor et al., 2005; Verslues and Sharma, 2010; Patriarca et al., 2021; Vujanovic et al., 2022). Being an osmolyte, proline can maintain cellular metabolism and even reduce plant growth in stressful conditions (Maggio et al., 2002; Vujanovic et al., 2022). This physiological response of proline accumulation would involve tRNA^Pro activity and could thus be a reason behind the prevalence of tRNA^Pro clusters and duplications ( Supplementary File 8 ). While these clusters might be initially repressed by methylation (Hummel et al., 2020), the plant stress response could induce the removal of these epigenetic marks, thereby increasing global tRNA^Pro transcription levels. To confirm this link, future studies are encouraged to look into the expression profile of these clustered tDNAs in plants. The potential biological functions of these tDNA clusters themselves may also be investigated further by future studies.

Another interesting observation is the apparent lack of certain tRNA isotypes in the nuclear genome of H. annuus and S. alba, even though their organellar counterparts are present. After further investigation, we found that prior to filtering via EukHighConfidenceFilter, H. annuus and S. alba had 117 tRNA^Gly and 82 tRNA^Asp predicted genes, respectively. None of these first-pass tRNA genes had an isotype score that met the cutoff for EukHighConfidenceFilter, which was 95 by default for these two isotypes. The tRNAscan-SE developers emphasized to only change the cutoff values with great caution, as they have already been tested on different large eukaryotic genomes (Chan and Lowe, 2019); thus, throughout our analysis, we opted to keep all default cutoff values unchanged. However, the fact that some of the first-pass tRNA^Gly and tRNA^Asp genes had scores that were very close to the cutoff value (as close as 94.5) indicates the need to reevaluate these score cutoffs.

To transcribe plant tRNAs, RNA polymerase III (Pol III) is recruited. One of the requirements for its recruitment is a TATA-binding protein (TBP), and the presence of TATA-box motifs upstream of plant tRNA genes is implicated in the efficiency of tRNA transcription (Dieci et al., 2006; Michaud et al., 2011). However, the proportion of angiosperm tDNAs containing such a motif is strikingly low ( Table 3 ). Previous studies have similarly reported the lack of TATA-box motifs upstream of many eukaryotic tDNAs (Hamada et al., 2001; Giuliodori et al., 2003; Dieci et al., 2006) as well as the little effect caused by the removal of TATA-box motifs in the transcription of plant tRNA^Leu genes (Choisne et al., 1998). For many Pol III-transcribed genes, TBP can be recruited without a specific TATA-like sequence. For these TATA-less genes, recruiting Pol III is instead facilitated by TFIIIC, which binds the DNA via the A and B boxes and recruits TFIIIB, which has a TBP as one of its subunits. TFIIIB recruits Pol III (Choisne et al., 1998; Yukawa et al., 2000; Dieci et al., 2006). This suggests that while some plants prefer the TATA-mediated recruitment of TBP [e.g., A. thaliana (Choisne et al., 1998; Hamada et al., 2001)], it may not be preferred or deemed necessary by other organisms that lack conserved TATA-box motifs. Dieci et al. (2006) hinted that the difference between a TATA-box-dependent and a TATA-box-independent organism might be found in their respective transcription machinery. Notably, the intragenic A and B boxes bound by TFIIIC were found in all detected nuclear tRNA genes ( Supplementary Files 2 and 3 ). However, this can mainly be explained by the fact that the tRNA D- and T-loops are encoded within these boxes (Galli et al., 1981; Hofstetter et al., 1981; Turowski and Tollervey, 2016) and that the tRNAscan-SE program detects tRNA genes based on the presence of A and B box sequences (Lowe and Eddy, 1997).

The CAA motifs, on the other hand, were found in most angiosperm tDNAs between positions −1 and −50 bp ( Table 3 ). Removal of these motifs upstream of plant tDNAs decreased in vitro expression levels of these tRNAs (Choisne et al., 1998; Yukawa et al., 2000). While previous studies reported functional CAA motifs to be between −1 and −10 bp in plant tDNAs (Yukawa et al., 2000, 2011; Michaud et al., 2011), more CAA motifs were found when the scope was extended up to −50 bp ( Table 3 ). This suggests that transcription start sites (TSS) for many plant tDNAs may be further upstream than others.

The majority of angiosperm tDNAs contained at least one downstream stretch of T residues ( Table 4 ), which is expected as it is considered an essential signal used by Pol III for transcription termination (Braglia et al., 2005; Arimbasseri and Maraia, 2015). In eukaryotic tRNAs, this poly(T) signal is commonly found to be approximately four to five bases long (Braglia et al., 2005). Aside from stretches of four to five T residues, there is also an abundance of poly(T) stretches that are 6 to 10 bases long, and those with extreme lengths—19, 26, and 23 bases—were found in the ANA, eudicot, and monocot tDNAs, respectively. While a significant percentage of angiosperm tDNAs do not contain a downstream poly(T) signal ( Table 4 ), it is possible that increasing the coverage to 100 or more nucleotides downstream (instead of only 50) will locate more poly(T) signals, backup poly(T) signals, and other poly(T) signals of extreme and variable lengths.

Our results provide a comprehensive overview of the tRNA gene content, structure, and organization of nuclear and organellar angiosperm genomes, utilizing the recent abundance of genomic data enabled by next-generation sequencing technologies. This study can thus supplement further studies on plant tRNA gene function and regulation. The specific function of these tRNA gene clusters and an explanation for the differences in the abundance of several regulatory motifs [e.g., TATA-boxes, CAA motifs, and poly(T) stretches] are some points that may be explored in the future.