SARS-CoV-2 Displays a Suboptimal Codon Usage Bias for Efficient Translation in Human Cells Diverted by Hijacking the tRNA Epitranscriptome

1. Introduction: The Critical Role of Codon Bias in Translation Efficiency

The efficiency of genetic code translation into amino acid is of utmost importance for rapidly dividing organisms. This tedious process is carried out by the ribosome and necessitates accurate decoding of the genetic code in messenger RNA (mRNA), through the selection of transfer RNAs (tRNA). The standard genetic code is degenerate: with the exception of tryptophan (Trp) and methionine (Met), most amino acids are encoded by two or more synonymous codons, which are used unequally within a genome (Figure 1A) [1]. This non-random distribution of synonymous codons, known as codon usage bias (CUB) [2], evolved alongside the expansion of the genetic code, shaped by evolution and natural selection [3,4]. CUB affects translation efficiency based on the availability of cognate codon-specific tRNAs, the effector molecules of translation, to ensure incorporation of the correct amino acids during polypeptide synthesis through complementary codon–anticodon base pairing [5,6]. CUB varies in a given genome. Highly expressed genes typically display stronger CUB, aligning with abundant tRNAs for optimal translation [7,8]. Conversely, rare codons may be selected on purpose, particularly at the beginning of coding sequences in eukaryotes. Indeed low CUB in the 5′ region decreases the translation elongation rate and reduces the likelihood of ribosomal traffic jams as translation progresses towards the 3′ end [9]. Additionally, reducing the elongation rate can promote the recruitment of chaperons facilitating co-translational protein folding [10,11]. Notably, CUB influences gene expression beyond translation, affecting transcription efficiency [12], termination [13], chromatin structure, and splicing [14]. As CUB varies depending on the organism [15,16], it is crucial for investigating the adaptation of infectious agents in their hosts, as well as virus evolution and pathogenesis. This report focuses on these aspects using SARS-CoV-2 as a case study.

Besides the three codons that signal translation termination ((UAA, UAG or UGA), 61 codons in mRNA are decoded through sequence complementarity with tRNA anticodon. Human cells display over 270 isodecoder genes (tRNAs with the same anticodon but different body sequence) among more than 610 annotated tRNA genes. Yet there are only 49 isoacceptor families (tRNAs with distinct anticodons but incorporating the same amino acid) to decode the 21 amino acids specified by the genetic code. Thus, the efficient deciphering of the 61 amino acid codons implies that some tRNAs can recognize more than one codon.

This conundrum led to the wobble hypothesis, introduced by Francis Crick in 1966 [17], proposing that only the first two bases of the codon pair precisely with corresponding bases in the anticodon, while the third position allows for flexibility or “wobble”. Accordingly, 30–40% of all codon recognition in a given organism is achieved through tRNA wobble recognition [18]. The modified wobble hypothesis of 1991 [19,20,21] expanded on the original hypothesis by including the role of certain base modifications occurring in or near the tRNA anticodon loop (Figure 1B). These modifications, ranging from simple methyl groups, to more complex structures such as sugars, affect either the folding of the tRNA [22] or its ability to bind to codons [23]. Progressively, the impact of these modifications in almost all steps of RNA metabolism has emerged [24], and their potential consequences for translation fidelity are being characterized [21]. Currently, these modifications are recognized as major architects of the anticodon structure, capable of preventing, favoring or expanding wobble base pairing, depending on their nature and position on the tRNA [25,26].

Recently, modifications in the tRNAs anticodon-loop, particularly at the U₃₄ position of tRNA have emerged as a central pillar in controlling codon bias interpretation by the translational machinery, impacting translation fidelity [27]. Position 34 in the anticodon loop of tRNAs has indeed been identified as a hot spot for base modifications. The addition of a methoxycarbonylmethyl group at position 5 and a thiol group at position 2 (mcm⁵s²) of U₃₄ in tRNA^Lys_UUU, tRNA^Gln_UUG, tRNA^Glu_UUC and tRNA^Arg_UCU, by the stepwise action of the Elongator (ELP1-6), ALKBH8, CTU1/2 multi-enzyme complex [28] enhances base pairing with A-ending codons that are infrequent in the human genome (Lys AAA; Gln CAA; Glu GAA; Arg AGA codons) (Figure 1C). The first two enzymatic steps (Elongator and ALKBH8) generate mcm⁵ modifications at the U₃₄ position of two additional tRNA (tRNA^Gly_UCC and tRNA^Arg_UCU) whereas Elongator alone can attach carbamoylmethyl (ncm⁵) to U₃₄ of six different tRNA (tRNA^Ala_UGC, tRNA^Thr_UGU, tRNA^Pro_UGG, tRNA^Ser_UGA, tRNA^Val_UAC, tRNA^Leu_UmAA). Additionally, KIAA1456 [29] (a human Trm9 homolog) has been recently shown to generate mcm⁵U directly from uridine [30]. In humans, additional modifications at position 34 include pseudouridine (Ψ), inosine (I), methycytidine (m⁵C) and queuosine (Q). Without these tRNA modifications, a codon-specific slowdown of translation occurs, impacting the overall translation efficiency of mRNAs enriched in corresponding codons [31,32,33].

1.2. Codon Bias and tRNA Pool

Translation efficiency cannot be explained solely by CUB, since preferentially used codons are not necessarily translated faster than non-preferred ones [34]. The overall control of translation efficiency also depends on the relationship between codon usage and the concentration of cognate tRNA. In any given organism, the set of tRNA molecules, known as the tRNA pool, is fitted to the corresponding genome’s codon bias to ensure optimal translation. This mutual alignment, measurable by the tRNA adaptation index (tAI, defined below) [35,36], defines codon optimality, indicating how efficiently a codon is translated by the ribosome. Indeed, adapting CUB to the most abundant tRNAs, decreases the time required for selecting the correct tRNAs, thereby reducing the likelihood of binding non-cognate tRNAs. In most prokaryotes and eukaryotes, there is a close correlation between tRNA level and the efficiencies of each codon–anticodon pairing. Consistent with translational selection, most optimal codons also have abundant corresponding tRNA copies in the human genome, although some still require wobble tRNAs. Yet, a close correlation between codon usage and tRNA abundance is not consistently observed, especially in multicellular organisms with higher tRNA gene redundancy [9]. In mammals, direct adaptation to anticodon pools does not fully explain synonymous codon usage, which is also shaped by mutational biases and genetic drift [35], such as GC-biased gene conversion [37]. Despite weak correlations between synonymous CUB and tRNA pools [38], the global correlation between codon and anticodon pools in mammals is strong and stable across different cell types and cell states [39,40]. Obviously, the extent of base modification in tRNAs also influences this balance especially for codons decoded through tRNA wobble recognition.

1.3. Viral Manipulation of Host Translational Machinery through tRNA Modifications

Viruses are obligatory parasites that rely entirely on the host’s translation machinery to translate their genome [41]. Numerous analyses comparing virus and host cell genomes in terms of codon usage have revealed that they use synonymous codons at different frequencies. Since discrepancy between the viral codon usage and the availability of the corresponding tRNAs in the host tRNA pool induces ribosome pausing [27], the central question we would like to bring to light is how viruses overcome the inadequacy between their own CUB and the suboptimal tRNA pool composition of the cells they infect, to efficiently translate their genome. In this context, viruses may need to maneuver host tRNAs to decode their skewed codons and optimize translation [42,43,44,45]. Moreover, because the tRNA pool in each cell type contains a subset of all the isodecoder and isoacceptor tRNAs needed for correct amino acids incorporation through complementary codon–anticodon base pairing, the differential expression of tRNA genes across tissues and individuals [46,47] may also influence viral tropism.

2. Results and Discussion

2.1. Analysis of SARS-CoV-2 Codon Bias

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the causative agent of the recent devastating coronavirus disease 2019 (COVID-19) pandemic [48,49], which has infected over 600 million people and caused more than 6 million deaths worldwide (https://covid19.who.int). Coronaviruses belong to the order “Nidovirales”, with those infecting humans falling into two genera: alphacoronaviruses (HCoV-229E and HCoVNL63), and betacoronaviruses (HCoV-HKU1, HCoV-OC43, Middle East respiratory syndrome coronavirus (MERS-CoV), the severe acute respiratory syndrome coronavirus (SARS-CoV1 and SARS-CoV-2). Since SARS-CoV-2 belongs to the same genus as SARS-CoV-1 and MERS-CoV, they share several structural characteristics [50] briefly outlined below [51,52].

SARS-CoV-2 is an enveloped virus with a positive-sense, single-stranded RNA genome of ~30 kb. Upon cell entry, two overlapping open reading frames (ORFs) ORF1a and ORF1b are translated from the positive strand genomic RNA (Figure 2A). The translation of ORF1b involves a -1 ribosomal frameshift enabling read-through of the ORF1a stop codon. ORF1a and ORF1b encode continuous polypeptides that are processed into 16 nonstructural proteins (Nsps) [53,54]. The viral genome is then used by the viral RNA-dependent RNA polymerase, Nsp12, to produce negative-strand RNA intermediates that serve as templates for the synthesis of positive-strand genomic RNA and subgenomic RNAs [55,56]. The subgenomic transcripts contain a common 5' leader sequence fused to different segments from the 3′ end of the viral genome [57], along with a 5′ cap structure [58] and a 3′ poly(A) tail [59]. They encode four conserved structural proteins: spike protein (S), envelope protein (E), membrane or matrix protein (M), nucleocapsid protein (N), along with several accessory proteins. By homology with SARS-CoV1, SARS-CoV-2 is predicted to translate 9 accessory proteins [60]. Nevertheless, the current annotation of SARS-CoV-2 (Reference Sequence: NC_045512.2) includes only 6 accessory proteins (3a, 6, 7a, 7b, 8, and 10), and not all of these ORFs have been experimentally reproducibly confirmed [61,62]. Using ribosome profiling techniques, the accurate quantification of canonical viral ORFs expression was established, along with 23 novel unannotated viral ORFs [63]. These include several in-frame internal ORFs lying within existing ORFs, resulting in N-terminally truncated products, and internal out-of-frame ORFs producing novel polypeptides.

2.1.1. Codon Bias Analysis: A Tool to Shed Light on Virus History, Origins and Evolution

It is well established that mutation pressure and natural selection are the primary factors shaping the codon usage of an organism [64]. Codon usage can also be influenced by nucleotide composition, synonymous substitution rate, gene length, expression level, and transfer RNA (tRNA) abundance [65,66]. Due to their reliance on host tRNA, viruses may evolve their codon usage to optimize or deoptimize translation in relation to their host’s codon usage [67,68]. Therefore, exploring the codon usage of viral genes is critical for uncovering viral evolutionary history [69], understanding virus-host interactions, and identifying the evolutionary forces shaping viral genomes [70,71]. Such information can also help characterize newly emerging viruses and trace their propagation across different host species.

Codon analysis toolbox

The level of bias can be investigated according to widely used bioinformatic tools:

Codon Adaptation Index (CAI) - The CAI estimates the degree of adaptation between a gene segment or an entire virus genome, and a potential host [72,73,74]. It compares the codon usage in the viral sequence to a reference set of highly expressed genes from the host which are assumed to use the most optimal codons for the host’s translational machinery. The CAI value ranges from zero to one, with one indicating that a gene uses the most frequently synonymous codon for each encoded amino acid, and values close to zero indicating the use of less common synonymous codons. Higher CAI values between different genes on genome segments, indicate a better adaptation to the host’s cell translational machinery. However, the CAI index is not strictly speaking a measure of CUB as codon usage is inherently multivariate and requires complementary approaches for comprehensive analysis.

Plotting the effective number of codons (Nc plot) – Nc plots graph the effective number of codons used (Nc) against the G+C frequency at the third base position of the codon (GC3s). This quantifies how far a gene’s codon usage deviates from equal usage of synonymous codons [75]. Neutrality plots, which plots GC1+2 (mean G+C frequency at the 1st and 2nd position) against GC3 for each gene, are commonly used alongside Nc plots to estimate the respective contribution of mutation pressure and natural selection in shaping the CUB [76,77,78]. The slope of the curve in neutrality plots indicates the percentage contribution of mutual pressure to the overall codon bias.

Relative Synonymous Codon Usage (RSCU) - RSCU determines the intrinsic preference of a given cell or virus gene for synonymous codons by calculating the ratio of observed codon frequency to the expected frequency, assuming equal usage of all synonymous codons for the same amino acid [79,80]. The codon bias is considered positive for RSCU > 1.6 and negative for RSCU < 0.6, and unbiased for values in between.

tRNA Adaptation Index (tAI) - The tAI quantifies translational efficiency in a given context by considering the availability of tRNAs for each codon, factoring in wobble base pairing efficiency (Wi) [81]. The normalized Wi values indicate the nominal speed at which the ribosome translates a codon relative to tRNA concentration, gene copies and pairing efficiency. The tAI of a gene reflects the adaptation of its coding sequence to the intracellular tRNA pool, providing a measure of translational efficiency complementary to CAI, Nc and RSCU [82]. If tRNA levels are not available, tRNA gene copy number can be used instead [36,83] and retrieved from dedicated database GtRNAdb [84] (http://gtrnadb.ucsc.edu).

2.1.2. SARS-CoV-2 Adaptation to Various Species

Several reports based on state-of-the-art bioinformatics suggest that SARS-CoV-2 originated from an ancestral coronavirus in bats, specifically related to the bat RaTG13 coronavirus [85], that likely passed through intermediate hosts such as pangolin (Manis javanica) before crossing species barriers again to infect humans [86,87,88]. We performed CAI [73] comparisons of the Wuhan reference strain NC_045512 of SARS-CoV-2 (weighted by the size of each ORF) against codon usage table (CUT) of several species from the HIVE-CUTs database (https://hive.biochemistry.gwu.edu/cuts/). Our analysis confirmed a high degree of SARS-CoV-2 codon usage adaptation to both pangolin and bat (Figure 2B), supporting their proposed role as animal sources for this emerging virus. Interestingly, our analysis also revealed an even higher adaptation index to snakes (Naja atra among others) and to the marsupial wombat (Vombatus ursinus), suggesting these animals may also be relevant in understanding the virus's transmission. Interestingly, the hypothesis that snake could serve as a potential intermediate host of SARS-CoV-2 between bats and humans has been debated [89,90,91] and remains plausible considering that snakes are common wildlife meal of Chinese people and are ordinarily sold at the Wuhan seafood and animal market [92]. The wombat seems less likely since it only lives in Australia. However, the recent trend among wealthy Chinese people to hunt wombats in Australian lodges [93,94] makes it a tiny but potential virus transmission niche. Global CAI analysis of SARS-CoV-2 ORFs shows oscillations around the 0.69 mark (Figure 2A), indicating that none of the multiple ORFs are fully adapted to their human host codon preferences, with some regions showing poor adaptation with CAI values below 0.6, compared to CAIs of 0.869, 0.882, and 0.819 for the well-adapted highly expressed β-globin, β-myosin and β-tubulin humans genes with protein abundance above 10,000 ppm, and CAI values below 0.78 for the poorly expressed RHA, RIG-I and Kallmann Syndrome human genes (protein abundances below 300 ppm) [95].

2.1.3. Nc and Neutrality Plots

If GC3s is the only determinant factor shaping the codon usage, the Nc values would align with a dotted theoretical curve representing random codon usage [96] (Figure 3A and 3B). If G+C compositional constraint alone influences the codon usage the gene of interest, then the GC3s and Nc correlated spots would lie on or below the expected curve. Nc values range from 20 to 61. A value of 20 indicates a maximum level of codon bias, whereas a value of 61 indicates a complete lack of bias. In general, if the Nc value is ≤ 35, the coding sequence is considered to have significant codon usage bias. The Nc-GC3s plot has been widely used to determine whether codon usages of given genes are affected by mutation only (corresponding points would lie around the expected curve) or also by other factors such as selection (corresponding points would depart away from, considerably below the expected curve). When the Nc and GC3s values were plotted for the individual SARS-CoV-2 ORFs (Figure 3A, right), all points, except for ORFs 7a and 10, lied below the expected random curve, indicating that G+C compositional constraints might significantly influence SARS-CoV-2 codon usage. The deviation of almost all SARS-CoV-2 ORFs from the random curve towards the lower Nc-values suggests that not only mutation but also other factors, such as translational selection, are likely to be involved in determining the selective constraints on codon bias in SARS-CoV-2 genes towards a preferred and limited set of codons. As a reference to human genome, translational selection is much more pronounced in housekeeping genes [97] such as Globin, Myosin and Tubulin than in poorly expressed humans genes (DHX9, DDX58, IFN-beta, or KAL) as outlined in the corresponding Nc plot (Figure 3A, center). Accordingly, the weighted average for all SARS-CoV-2 ORFs (black dot with Nc=45.1 and GC3s=0.28 in Figure 3A) markedly differs from other RNA viruses such as the Flaviviridae Zika virus ORFs (Figure 3A, left), that also have a bias towards A-ending codons [98,99,100], showing a much lower GC content and overall Nc value for SARS-CoV-2, stressing again the unique codon-wise characteristics of SARS-CoV-2. The same analysis was performed for Orf1a and Spike gene segments of various Coronavirus (Figure 3B) including SARS-CoV-2 (Wuhan initial reference isolate and Omicron recent isolate), SARS-CoV-1, MERC-CoV and various hCoV lineages (229E, OC43 and NL63). We observed that all points representing different strains were lower than the standard curve for both genes, Orf1a showing more dispersion on the Nc axis than Spike. Additionally, with the exception of SARS-CoV-2 Orf1a (Wuhan and Omicron), the coronavirus strains were not clustered together, highlighting again that mutational pressure combined with other factors may have contributed to the codon usage bias of SARS-CoV-2. Spike Nc plot shows more clustering between strains revealing less bias between strains for this gene segment, with a potential higher contribution of mutational pressure to the Spike codon bias. This was further confirmed with the complementary neutrality plot analyses (Figure 3C) that revealed the greater contribution of mutation pressure to Spike than to Orf1a CUBs (Spike:32% >> Orf1a:23%), while the relative constraints on GC3 (natural selection) being the main evolutionary force driving CUB is conversely higher for Orf1a (Orf1a:77% >> Spike:68%).

2.2. SARS-CoV-2 Genome is Enriched in U₃₄-Sensitive Codons

With the lack of information on the level of virus-mediated alteration of U₃₄ tRNA modification in mind, we compared U₃₄ sensitive codon frequencies in coronaviruses (HCOV-OC43, SARS-CoV-1, MERS-CoV and SARS-CoV-2) and human genomes using the RSCU indicator [15]. A striking contrast emerged when these frequencies were arranged in clusters (Figure 4). Codons preferred by SARS-CoV-2 are (i) barely used in human cells and (ii) predominantly include codons requiring U₃₄ modifications on their cognate tRNAs for efficient decoding (with the exception of Gly (GGA)). Notably, SARS-CoV-2 exhibits a greater divergence from human codon frequencies compared with SARS-CoV-1 and MERS-CoV. Although the latter still shares some codon usage preferences with humans, SARS-CoV-2 contains a clear enrichment in U₃₄-sensitive codons. This pattern was consistent across both nonstructural (Orf1a and Orf1b) or structural viral genes (Spike). U₃₄-sensitive codon usage between SARS-CoV-2 and humans revealed a significant viral preference for U₃₄-sensitive codons with up to two-fold enrichment for some codons such as Gln^CAA, Arg^AGA, and Leu^TTA. (not shown).

2.2.1. Comparison of Coronavirus Translation Adaptation (tAI)

Previous CAI measurements have the disadvantage of relying on highly expressed reference host genes. The translation adaptation index (tAI) offers a more nuanced approach as it can be based on either intracellular tRNA levels (when available) or tRNA gene copy numbers in the host genome. Here, we compared tAI of various ORFs encoding non-structural (ORF1ab) and structural proteins (S, E, M and N) among coronaviruses infecting humans. Our analysis revealed SARS-CoV-2 as the least adapted to the translational machinery of its human host, with all its ORFs having tAI below 0.34, excepted for the Matrix protein (M) (Figure 5A). This tAI level obtained for SARS-CoV-2 is significantly lower than that observed for human ORFs encoding highly abundant proteins (such as β-myosin, β-globin, or β-tubulin). Calculations were performed using human tRNA gene copy number retrieved from the genomic tRNA database (http://gtrnadb.ucsc.edu) known to mirror the global tRNA abundancy of a given organism but lacking information related to tissue specificity. However, it seems from the recent of the Nedialkova’s group [101] that tRNA transcript levels may fluctuate without affecting significantly tRNA anticodon pools abundancy and availability that primarily dictate decoding pace. Accordingly, tRNA anticodon pool was shown to remain stable across cell types (human primary cells (cardiomyocytes (CM), neuronal precursor cells (NPC) and mature neurons) differentiated from iPSC cells) ensuring consistent decoding rates throughout development, independently of cell identity. This stability across cell types of tRNA pool was further authenticated by plotting normalized tRNA gene copy number (GCN) against the normalized experimental data of tRNA anticodon expression level in these human primary cells (Figure 5B), showing a quasi-constant correlation between GCN and experimental tRNA levels in the four cell types. This observation validates the use of human tRNA GCN in the calculation of aforementioned tAI used in Figure 5A.

2.2.2. The Enigma of Spike Protein’s Furin Cleavage Site

Among the SARS-CoV-2 ORFs, the region encoding the spike protein has been extensively studied as it mediates attachment to the host cell by binding to the ACE2 membrane protein and facilitates viral fusion to the host cell membrane following efficient cleavage by furin proteases [102,103]. Early examination of the Wuhan SARS-CoV-2 isolate early revealed an unusual furin-like cleavage site at the S1/S2 junction of spike ORF [104,105]. This site, created by insertion of a 4-amino acid sequence (PRRA), corresponds to the insertion of 12 nucleotides (...U CCU CGG CGG GC...) absent from all other sarbecoviruses, including the closely related bat sarbecovirus RaTG13 with ~96% genome sequence identity [48]^, [106] (Figure 6A). In SARS-CoV-2, the furin site insertion lies in a region where codons are moderately adapted to the human host, as depicted by the CAI curve. Most codons in this region are commonly used in the virus genome (high RSCU), except for the arginine dicodon (CGG CGG (R R)). Interestingly, the corresponding unique CGG codon is less preferred than AGA in SARS-CoV-2 (SARS-CoV-2 RSCU_CGG=0.17, SARS-CoV-2 RSCU_AGA=2.63,)), which is not the case in the human genome (human RSCU_AGA =1.29, human RSCU_CGG =1.21) (Figure 6C). During the pandemic, mutations impacting the furin site have been rare, suggesting a strong purifying selection in humans [107,108]. However, non-arginine residues in the PRRAR motif remain permissive to potential optimization during human viral evolution in different variants of concern and interest (Figure 6B). This complex interplay between codon usage and furin site evolution warrants further investigations to unravel its significance in SARS-CoV2 origins.

2.3. Suitability of the SARS-CoV2 Highly Biased Codon Composition for Viral Translation in Target Tissues

The translation efficiency of the viral genome is heavily influenced by codon optimality, which is determined by the balance between the viral codon usage bias and the availability of a suitable tRNA pool in target cells. However, recent reports highlighted the significant variability in tRNA gene expression across human tissues [46,109,110,111]. This variability suggests the necessity to reevaluate codon optimality by considering tissue-specific codon usage in compliance with virus tropism [112]. Whereas the evaluation of the effect of tissue-specific codon optimality on viral protein synthesis may remain experimentally elusive, a recent in silico study by Hernandez-Alias et al. [113] analyzing the relative codon usage landscape over 500 human-infecting viruses alongside with tissue-specific tRNA expression profiles from 23 human tissues from The Cancer Genome Atlas (TCGA), has suggested that tRNA repertoires could determine tissue-specific translational efficiency [114]. They proposed that the SARS-CoV-2 proteome was well-adapted to tissues like the upper respiratory airways, lung alveoli [115,116] and gastrointestinal tract [117], seemingly matching the observed SARS-CoV-2 tropism. However, the translation appropriateness of these tissues was also matching with other viruses such as Flavivirus or Alphavirus that exhibit tropisms that do not share one shred of similarity with SARS-CoV-2 tropism. Several flaws of this study, including the use of tRNA abundance data derived from cancer cells in which profound deregulation of tRNA expression occurs [118,119,120,121], restrict its significance. Interestingly, recent studies showed that highly expressed genes in human lung primary tissue have a codon composition perfectly aligned with SARS-CoV-2 codon usage, suggesting that the virus might have optimized its codon bias to take advantage of lung cells [122]. The recent Nedialkova’s report [101] that illustrated the broad stability of tRNA anticodon pools across human cell types (shown above in Figure 5B), reinforces even more the prohibition of cancer cells in approaches aiming tRNA pool dynamics.

Taken together, these observations underscore the crucial need to further investigate how the codon composition of the viral ORFeome influences the translation rate of host genes and promotes viral translation.

2.3.1. tRNA Modifications: Do RNA Viruses Have a Wobble?

In analyzing virus-host interactions, it is crucial to consider the modification of host cell tRNAs promoted by the virus. While the functional importance of the virus mRNA epitranscriptome has been extensively reported in human cells [123,124] and more recently in the SARS-CoV-2 RNA genome and subgenomic transcripts [125,126,127,128], our understanding of virus-induced modifications in host tRNAs is still limited although increasingly relevant [129]. As outlined earlier, modifications of the anticodon loop of tRNA regulate translation rate and fidelity, contributing to the translational adaptation [130]. For instance, changes in tRNA modifications significantly impact translation efficiency in response to physiological stresses [131,132]. Hypomodification, occurring especially at position U₃₄ of specific tRNAs, namely tRNA^Lys_UUU; tRNA^Gln_UUG; tRNA^Glu_UUC, slows down translation, disrupts protein homeostasis and reduces cellular fitness. These mechanisms are crucial for maintaining cellular function and viability during stress until normal conditions are restored. While these mechanisms may have a critical influence on viral translation, this aspect of virus/host relationships remains largely unexplored.

2.3.2. SARS-CoV-2 Codon Bias Dynamics during the Pandemic

Besides providing crucial information on viral genetic adaptation to host [133,134], evaluating nucleotide composition and CUB in viral genome also provide further insights into the mutational patterns of viruses over time and can be crucial for vaccines and antiviral therapy development [135]. A recent study analyzing over 3.5 million SARS-CoV-2 sequences revealed unique mutational trends with consistent nucleotide and codon frequencies [136]. This study also highlighted variations over time, including synonymous mutations, silent at the amino acid level, and nonsynonymous mutations inducing amino acid changes, impacting protein levels. It also revealed an unexpectedly high proportion of nonsynonymous mutations in the Spike gene when compared with glycoprotein genes from other RNA viruses.

To deepen this analysis, we examined the sequences from the various SARS-CoV2 clade (Alpha, Beta/Mu, Delta, Lambda/Gamma and Omicron) extracted from the Fumagalli’s paper [136] and generated their respective RSCU profile, monthly from December 2019 to July 2022. The RSCU pattern will here give us a snapshot of the average CUB for each month (Figure 7A). Focusing on SARS-CoV-2 Orf1ab and Spike genes, we analyzed the overall codon bias patterns by cluster analysis in order to compare the temporal codon usage variations of each clade. Notably, all seven clades were present and detectable at the very beginning of the pandemic albeit with variable abundancy. The Alpha clade was the first to strongly emerge and spread worldwide, while Beta/Mu and Lambda/Gamma remained at low levels of diffusion. Regarding their respective CUBs, all clades exhibited very similar patterns during the first 3 months of 2020. From April to September 2020, the RSCU pattern of Alpha clade began to fluctuate markedly before reaching a long period of stability from October 2020 to September 2021, coinciding with the peak of Alpha clade expansion worldwide. Surprisingly, the CUB pattern during this high transmission period was reversed compared to the initial pattern. The pattern of a stabilized, reversed CUB was observed not only in the Alpha clade but also in the Delta and Omicron clades, when they both reached their highest levels of transmission.

The clear correlation between CUB inversion and peak of diffusion can be interpreted in two opposite ways: a) the high level of expansion has selected sub-variants with a stable CUB, or b) a highly adapted sub-variant with optimized CUB preceded the peak, driving the exponential expansion. Although both scenarios might be involved, the persistence of inverted CUB beyond peak diffusion suggests natural selection of an optimized CUB likely played a role (favoring option b).

Comparing the RSCU of both genes at peak expansion for each clade reveals that Spike exhibits a much stronger preference for codons ending by A/U, especially those recognized by tRNAs modified at U₃₄ position (e.g. Arg_AGA, Lys_AAA, Glu_GAA and Gln_CAA) (Figure 7B). This necessitates analyzing RSCU patterns for each gene individually at expansion peaks. Interestingly, each clade shows different CUB optimization pattern during peak diffusion, indicating that multiple codon patterns can support an optimal viral translation (Figure 7C and 7D). For instance, Omicron has evolved a codon bias diametrically opposed to the CUB of other clades, while Alpha and Delta Spike sequences share a similar CUB pattern (Figure 7C). Conversely, the Orf1ab gene in Alpha and Gamma/Lambda clades show the closest patterns (Figure 7D). Furthermore, within a given clade Orf1ab and Spike genes have developed distinct RSCU pattern. These varying preferences for A/U ending codons or codons requiring U₃₄-modified tRNAs, suggest that efficient translation might depend on subtle but specific factors like translation speed and protein folding parameters.

2.4. Hypothesis: Virus-Driven Manipulation of the tRNA Pool and tRNA Modifications Forces Translation of SARS-CoV-2 Genes

Bioinformatic approaches like the tAI and the CAI measures, reasonably predict gene expression but can be improved. More precise estimations of amino acid-loaded tRNAs (“ready-to-translate” tRNAs) availability would be more realistic than using concentration of tRNA molecules or its estimate from tRNA gene copy number. The real availability of functional tRNAs, influenced by tissue-specific pools, is required for high translation efficiency. Another important mechanism for efficient protein synthesis is the channeling effect, which involves the direct transfer of aminoacyl-tRNA/tRNA between aminoacyl-tRNA synthetases (ARS), elongation factors, and ribosomes. This process is crucial because it regenerates and concentrates frequently used tRNAs near specific translation sites [137,138]. Additionally, the global CUB measures (CAI, tAI, RSCU) do not consider the order of favorable and unfavorable codons along the transcript, which can create fast or slow translation segments [42,139]. Environment-dependent dynamic variations in the tRNA pool and tRNA demand, should also be integrated in future models of translation efficiency. Using tRNAs composition from cancer cells can distort prediction of virus/host-tissue translational compatibility due to altered tRNA levels [118,119,120,121,140] and tRNA modifying enzymes in pathological conditions [141]. U₃₄ modifications, in particular, rapidly respond to metabolic changes [142] such as methionine metabolism, carbon balance, or phosphate homeostasis [143,144,145].

Translational reprogramming is clearly appearing as a key element in cell adaptation to changing environments [146] and may aid virus adaptation to hosts. We have here highlighted striking discrepancies between SARS-CoV-2 preference for U₃₄-sensitive codons and availability of cognate U₃₄-modified tRNAs (ie: mcm⁵s²) in target cells. Since SARS-CoV-2 hasn't fully adjusted its codon bias to match human target cells, its recent rapid expansion may be due to its ability to manipulate U₃₄ modifications, optimizing translation and facilitating its infection cycle. By reducing the need for precise codon usage adaptation, this ability could allow the virus to infect a broader range of hosts [147].

2.5. Evidences Supporting the Ability of SARS-CoV-2 to Exploit the tRNA Epitranscriptome in Order to Favor Viral Translation

2.5.1. Potential Manipulation of tRNAs by SARS-CoV-2

Direct evidence supporting this hypothesis for SARS-CoV-2 is currently limited in the literature, except for a recent study by Tao Pan’s group [148], which reported the presence of cellular tRNA in SARS-CoV-2 particles. This study suggested that each virion contains at least four different tRNA molecules. Notably, among the eight tRNAs preferentially enriched in SARS-CoV-2 particles, 75% require U₃₄ modification for efficient decoding, including tRNA^Lys_UUU and tRNA^Glu_UUC, both bearing the mcm⁵s²U₃₄ modification. This observation suggests preferential packaging of critical tRNAs complementing skewed SARS-CoV-2 codons, reinforcing the assumption that U₃₄-sensitive codons in viral genes require adaptation of the host tRNA pool for efficient viral translation. Other enzymes involved in tRNA modifications might also be involved in this virus-mediated translational control, such as the tRNA methyltransferase TRMT1 that generates the m²₂G₂₆ mark. Indeed, recently we [149] and others [150,151] have shown that TRMT1 was specifically proteolyzed by the SARS-CoV-2 Nsp5 main protease leading to a decrease of the m²₂G₂₆ modification on tRNAs in infected cells negatively impacting viral replication. This manipulation, suggests a role of m²₂G₂₆ tRNA modification patterns in cellular pathogenesis and biology of SARS-CoV-2 infection.

This situation is not unique to SARS-CoV2. Recent studies have shown that during the Alphavirus CHIKV infection, deregulated expression of KIAA1456, an ALKBH8 homolog able to generate the mcm⁵U₃₄ modification of tRNAs, consequently reprogramed codon optimality and favored viral RNA translation [30]. This mechanism was also shown to occur during the Flavivirus DENV infection, which like CHIKV exhibits a strong enrichment in U₃₄-sensitive codons [30,129]. In parallel, a recent preprint [152] proposed that DENV is also able to exploit the host tRNA by a different mechanism involving the ALKBH1 RNA modifier and the virally-encoded NS5 protein (with dual RNA methyltransferase (MTase) and RNA-dependent RNA polymerase (RdRp) enzyme activities). The cellular and viral enzymes can both mediate f⁵Cm-modification of the cytoplasmic tRNA-Leu(CAA) at the wobble position C_34, and collaboratively promote pro-viral translational remodeling during DENV infection.

2.5.2. Experimental Data Revealing SARS-CoV-2 Induced tRNA Epitranscriptome Modulations

To assess the validity of our hypothesis, we investigated the dynamic changes in tRNA modifications within SARS-CoV-2-infected cells. We first explored by LC-MS/MS the behavior of epitranscriptomic marks in the tRNA subpopulation extracted from SARS-CoV-2 infected VeroE6 cells (Figure 8A). Using this approach, we tracked 21 different tRNA post-transcriptional modifications and discovered that modifications at position 34, including ncm⁵U, mcm⁵U and mcm⁵s²U, were noticeably increased. We extended our analysis to human Caco2 cells, a more suitable cellular model for exploring SARS-CoV-2 infections given the known COVID-19 gastrointestinal manifestations. Our focus on U₃₄ modifications once again highlighted the early and rapidly changing nature of the three marks we examined (Figure 8B). By simultaneously tracking tRNA levels, we also found that two of the four U₃₄-modified tRNAs (tGlu^UUG and tGln^UUC) were upregulated at later time points (Figure 8C). Overall, these preliminary results indicate the possible role of tRNA modifications in SARS-CoV-2 infection and reinforce the idea that SARS-CoV-2 can manipulate the host's tRNA transcriptome. However, it remains unclear whether these changes can genuinely benefit viral infection.

2.5.3. SARS-CoV-2 Infection is Impaired when the tRNA U34 Modification Pathway is Disrupted

In mammals, enzymes responsible for the chemical U₃₄ modification include the Elongator complex (Elp1–6), Alkylation repair homolog 8 (Alkbh8), and thiouridylases (Ctu1/Ctu2) [153]. If the induction of ncm⁵U, mcm⁵U and mcm⁵s²U marks indeed benefits viral translation, virus infection should be closely tied to the accurate activity of the Elongator complex. Using primary fibroblasts from patients with Familial Dysautonomia (FD) (Figure 9A), which lack Elp1 expression (=IKBKAP-/-) and consequently have reduced levels of all three U₃₄ tRNA modifications (Figure 9B) [154], we observed significantly lower levels of SARS-CoV-2 infection in FD cells compared to control fibroblasts (wt) from healthy individuals (Figure 9C). These preliminary data emphasize the critical role played by U₃₄ modifications of host tRNAs in SARS-CoV-2 lifecycle and provide the first evidence of SARS-CoV-2 ability to rewire the tRNA epitranscriptome to facilitate translation of its heavily codon-biased genome.

3. Concluding Remarks

3.1. Is altering tRNA Epitranscriptome a Common Viral Strategy?

Viruses are highly dependent on the host cell's translation machinery, including host tRNAs, for the efficient translation of their genetic material [42,155,156], and this process is heavily influenced by chemical modifications of tRNAs, catalyzed by various tRNA-modifying enzymes, particularly within the tRNA anticodon loop region [27,157]. We here showcased SARS-CoV-2's ability to directly target this crucial step by upregulating U₃₄ tRNA modifications to facilitate the translation of its genome, which is enriched in U₃₄-sensitive codons. In addition to Chikungunya and Dengue for which the ability to interfere with the host tRNA epitranscriptome has been suggested [30,129,152,158], Zika can have also evolved the same tactic to overcome its high degree of preference for U₃₄-sensitive codons [98]. Besides tRNA modification, translation efficiency favoring viral translation can also be manipulated by altering tRNA levels [159], as we briefly illustrated in the case of SARS-CoV-2. Retrovirus such as HIV [160] or DNA viruses, such as SV40 [161], EBV [162], Adenovirus [163] or HSV-1 [164], are able to manipulate tRNA levels by stimulating Pol III transcription of tRNA genes.

3.2. Future Priority Investigation Areas

New approaches to explore the dynamics of tRNA epitranscriptome during viral infections are imperative. The direct quantification of tRNA pools with new experimental approaches like the recently developed mim-tRNAseq, which simultaneous measures tRNA abundance and modifications [165,166] could provide insight alongside with Ribo-seq analysis of viral and cellular translatomes. Expansion of the toolbox for quantitative recording and understanding the chemical biology of the tRNA epitranscriptome is clearly needed. It will include emerging technologies for mass spectrometry-based [167,168], nanopore-based tRNA modification mapping [169,170] and analysis of ribosome-bound tRNAs [171]. In this respect, MLC-seq [172] (mass spectrometry ladder complementation sequencing), a recent groundbreaking mass spectrometry approach, offers a potential solution to these challenges by providing quantitative, and site-specific mapping of RNA modifications, revealing the truly complete informational content of tRNA. The recent development of DORQ-seq [173], a combination of cDNA hybridization and deep sequencing, will also deliver a detailed tRNA composition matrix from femtomolar amounts of total tRNA. Advances in nanopore sequencing are anticipated with optimized basecalling models that could allow enhanced detection of RNA modifications and mapping [174].

The essential role of tRNA modifications in translation regulation is clearly appearing as a new crucial host process hijacked by RNA viruses to adapt host translation to their sub-optimal codon bias. It is therefore decisive to investigate the underlying mechanisms involved, particularly those targeting U₃₄ tRNA modifications. By genetically inactivating U₃₄ enzymes using CRISPR/Cas9 or RNA interference, we can assess their impact on SARS-CoV-2 translation and replication. Additionally, exploring how individual viral proteins influence the abundance and activity of U₃₄ enzymes can provide valuable insights into the mechanisms of this new level of virus-host interaction.

Technological and scientific advances in RNA modification, highlighted their role in viral RNA structure, localization, splicing, stability, and translation [175,176]. Understanding viral- or host-induced alterations of RNA marks is essential for understanding gene regulation, identify essential marks for virus cycle and designing appropriate drugs. Nucleoside-derived inhibitors targeting SARS-CoV-2 nsp14 (N7-guanine)-methyltransferase have recently shown promise in crippling the stability of viral RNA [177]. Exploring virus–tRNA epitranscriptome interactions could open promising new avenues for therapeutic intervention.

4. Material and Methods

4.1. Bio-Informatics - Codon Analysis

Codon adaptation index (CAI) and Relative Synonymous Codon Usage (RSCU) calculations were done using CAIcal web available tools (http://genomes.urv.es/CAIcal) [73]. Nc data were obtained using CAIcal and plotted with Rstudio-ggplot2. Codon frequency were calculated with Codon Utilization Tool (CUT) website (http://pare.sunycnse.com/cut/index.jsp) [178]. tRNA adaptation index (tAI) was calculated using the stAIcalc software [179] (http://www.cs.tau.ac.il/&sim;tamirtul/stAIcalc/stAIcalc.html). Cluster analysis were done with Genesis 1.8.1 [180] and Cluster 3.0 [181] and visualization were made with Java Treeview [182]. Protein abundance were derived from PaxDb database version 4.1 (https://pax-db.org/) [183].

To verify the accuracy of tAI interpretation, we compared the consistency of tRNA expression across different human primary cell types by visually examining the relationship between tRNA gene copy number (GCN) and experimentally measured tRNA levels in various human primary cell types. We used a previously published dataset from the Nedialkova group [101] including addition information concerning modifications at position 34 occurring on some specific tRNAs.

4.2. SARS-CoV-2 Sequences

Accession numbers of the main SARS-CoV-2 sequences used in this report are listed in annex 1 and downloaded from NCBI database with the exception of SARS-CoV-2 FRA that was from the european virus archive (http://www.european-virus-archive.com) and corresponded to the exact 2020 isolate from Paris-Ile-de-France we used in our infection experiments. Analysis of clade-related RSCU evolution of SARS-CoV-2 during the pandemic, was performed using sequences recovered from the dataset generated in the Fumagalli’s paper [136].

4.3. Experimental Data

4.3.1. Cells and Viruses

Patient primary fibroblasts were from Coriell Institute: GMO1652 derived from non-FD control (Skin fibroblast (arm) from 11 years old Caucasian female); GMO4959 derived from FD patient (Skin fibroblast (arm) from 10 years old Caucasian female). The SARS-CoV-2 was a French Ile de France isolate (www.european-virus-archive.com/virus/sars-cov-2-isolate-betacovfranceidf03722020). Viral stocks were generated by amplification on VeroE6 cells (epithelial kidney of an African green monkey, ATCC CRL-1586). The supernatant was collected, filtered through a 0.45 µm membrane, and tittered using a TCID50 assay. Caco2 cells (epithelial colon adenocarcinoma, ATCC HTB-37) were used for tRNA modifications quantification upon SARS-CoV-2 infection.

4.3.2. Quantification of tRNA Modifications by Mass Spectrometry (LC-MS/MS)

RNA preparations enriched in tRNAs were obtained using mirVana™ miRNA Isolation Kit (Thermo). RNA samples were then digested by Nuclease P1 and treated by Alkaline phosphatase. Samples were then injected into LC-MS/MS. The nucleosides were separated by reverse phase ultra-performance liquid chromatography on a C18 column with online mass spectrometry detection using Agilent 6490 triple-quadrupole LC mass spectrometer in multiple reactions monitoring (MRM) positive electrospray ionization (ESI) mode. Quantification was performed by comparing with the standard curve obtained from pure nucleoside standards running with the same batch of samples.

4.3.3. Assays for Viral Replication

For infections, the cells were previously transduced with a Lentiviral vector expressing ACE2 using the lentiviral construct RRL.sin.cPPT.SFFV/Ace2.WPRE (MT136) kindly provided by Caroline Goujon (Addgene plasmid # 145842) [176]. Seventy-two hours after transduction, accurate ACE2 expression was controlled on western blot probed with anti-ACE2 antibody (Human ACE-2 Antibody, AF933, R&D systems). ACE2-positive cells (70-80% confluence) were then infected with SARS-CoV-2 diluted to achieve the desired MOI. After 24 hours in culture, the cells were lysed with the Luna cell ready lysis module (New England Biolabs). The amplification reaction was run on LightcyclerR 480 thermocycler (Roche Diagnostics) using the Luna Universal One-Step RT-qPCR kit (New England Biolabs), and SARS_For: 5’-ACAGGTACGTTAATAGTTAATAGCGT; SARS_Rev: 5’-ATATTGCAGCAGTACGCACACA; GAPDH_For: 5’-GCTCACCGGCATGGCCTTTCGCGT and GAPDH_Rev: 5’-TGGAGGAGTGGGTGTCGCTGTTGA primers. Each qPCR was performed in triplicate, and the means and standard deviations were calculated. Relative quantification of data obtained from RT-qPCR was used to determine changes in SARS-CoV-2 gene expression across multiple samples after normalization to the internal reference GAPDH gene. Individual tRNA quantification was performed by RTqPCR using tRNA gene specific primers and miR103a primers for internal normalization purposesi (see Appendix B).