High frequency of transition to transversion ratio in the stem region of RNA secondary structure of untranslated region of SARS-CoV-2

Stem and loop annotations of SARS-CoV-2 reference genome

The reference SARS-CoV-2 isolate Wuhan-Hu-1 genome consisting of 29,903 bases (NC_045512.2) has been used for annotations of functional regions of the SARS-CoV-2 (Wu et al., 2020). For the base substitution study, the untranslated regions, 5′-UTR (base position 1 to 265) and 3′-UTR (base position 29,675 to 29,903) were considered. Unlike the protein coding gene sequences, these UTRs are devoid of any translational selection on codon usage bias (Sharp & Li, 1986). The stable secondary structure reported by Miao et al. (2021), covered the 5′-UTR along with a portion of the nonstructural protein gene nsp1. This 5′-UTR secondary structure had been determined experimentally using radio-labeled transcript and RNase V1 enzymatic probing (Miao et al., 2021). The secondary structure reported by Huston et al. (2021) included a portion of the structural protein gene N, accessory protein gene ORF10, and the 3′-UTR. The well-defined stem-loop motifs (Table 1) of UTRs, except those paired with the bases of the neighboring coding regions were considered for the base substitution study. In subsequent sections in this article, these stem-loop structures are abbreviated as SL-I through SL-VII (Table 1).

Out of these seven stem-loop motifs given in Table 1, the first six motifs SL-I through SL-VI were from 5′-UTR, and the last motif SL-VII was from 3′-UTR. The stem-loop motifs in terms of dot-bracket notations given in Table 1 were visualized using RNAfold Forna software (http://rna.tbi.univie.ac.at/forna/) (Gruber et al., 2008; Kerpedjiev, Hammer & Hofacker, 2015; Lorenz et al., 2011; Mathews et al., 2004). The two stem-loop structures SL-II and SL-III were pseudoknotted. In SL-II, base positions 51, 52, and 53 in the loop regions were paired with positions 37, 38, and 39. Base positions 67, 68, 69, and 70 of the SL-III were paired with base positions 76, 77, 78, and 79. The secondary structures of the 5′-UTR and 3′-UTR considered in this study are given in Figs. S1 and S2, respectively. For the base substitution analysis, bases were categorized into two groups (i) paired or belonging to stem region and (ii) unpaired or belonging to loop regions. Out of the 204 positions analyzed in the 5′-UTR, 149 were categorized under the stem region, and the remaining 55 were categorized under the loop region. Out of the 156 positions analyzed in the 3′-UTR, 68 were categorized under the stem region, and remaining 88 were categorized under the loop region. In total, the percentage of bases considered under stem and loop were found to be 60.0 and 40.0, respectively.

Retrieval and sequence alignment of UTRs of SARS-CoV-2 genome

In the present investigation, 46,076 high-coverage SARS-CoV-2 genome sequences were extracted on 24th July 2020 from the GISAID database (https://www.gisaid.org/) (Shu & McCauley, 2017). These genome sequences were sampled from patients, drawn from 95 countries, across the globe. These genome sequences represent the early stage of adaptation phase of SARS-CoV-2 pandemic in the human population. The downloaded genome sequences were processed to filter out sequences displaying size mismatch with the reference sequence NC_045512.2, including those with ambiguous nucleotides other than A/T/G/C. The final filtered set of 42,725 strains was retained and used to create a local BLAST database. Using 5′-UTR and 3′-UTR sequences of the reference genome (NC_045512.2) as query sequences, alignments of the two (5′- and 3′) untranslated regions were extracted from the local BLAST database, for base substitution analysis. In total, 4,049 sequences of the 5′-UTR and 2,811 strains of the 3′-UTR were available for the analysis. Alignments of the 5′-UTR and 3′-UTR sequences are given in Table S1. In addition to the above strains, dominant variants of Delta and Omicron strains reported in the GISAID database up to 24th September 2023, were also analyzed. After preprocessing, 8,227 sequences of the 3′-UTR and 7,020 sequences of the 5′-UTR, were available for base substitution analysis (Table S2).

Identification of base substitutions in the stem and loop motifs of RNA secondary structure of SARS-CoV-2

To identify inter-species mutations in an alignment of homologous sequences of a few closely related species, researchers have often used methodologies based upon reconstructing a phylogenetic tree, and changes from ancestral sequences at various tree nodes (Wu & Maeda, 1987). Taking advantage of the large volume of SARS-CoV-2 genome sequences available in the public domain, a simple approach was employed in this intra-species base substitution study. A consensus sequence considering the most frequent nucleotide at each position in the aligned sequences was generated. This consensus sequence was then compared with each sequence to identify base substitutions. Identification of base substitutions was carried out using the consensus sequence as shown in the hypothetical example (Fig. 3). The mutation frequencies were further normalized by dividing the total count of a given mutation, by the total number of nucleotides, in which the mutation occurred. For example, if C → U substitution count was found to be 1 in a sequence and the count of base C in that sequence was 10, then the normalized mutation frequency was calculated as 1/10 = 0.10. This consensus sequence-based method for estimating intra-species base substitutions is reported to be quite effective in mutation studies in bacteria genome tRNA gene secondary structures (Sen et al., 2022), estimating dN/dS for protein-coding genes (Aziz et al., 2022) and polymorphism analysis in intergenic regions (Beura et al., 2023). In subsequent sections, these normalized mutation frequency values are referred to as mutation frequencies. Each substitution in a sequence was further mapped to the stem-loop structure, and classified into loop and stem regions. If a pair of substitutions were observed in a paired position in the stem region, then they were considered compensatory; otherwise, they were designated non-compensatory substitutions. It is significant at this juncture, that the compensatory substitutions can be considered relatively older and more stable than the non-compensatory substitutions (Higgs, 2000). Therefore, the non-compensatory substitutions were found to be more frequent than the compensatory ones. Among all the stem-loop structures, only the non-compensatory substitutions were analyzed in this investigation Python scripts were written for identifying substitutions in the alignment of the SARS-CoV-2 secondary structures and their categorizations. The Python script, along with the executable and supporting stem-loop motif sequence files are available online for researchers in GitHub (https://github.com/MDash-NITAP/SLanalysis.git).

Higher mutation rate in the loop than the stem in the RNA secondary structure of SARS-CoV-2

Assuming the deleterious effect of mutations in the stem region, a comparative analysis between the rate of mutation in the stem and the loop region in the seven stem-loop motif sequences was carried out (Table 1). The unique mutations found in the stem-loop motif SL-I, SL-III, and SL-IV are shown in Fig. 4, and the remaining motif mutations are given in Table S3. In total, 360 base positions have been considered in the 5′-UTR and 3′-UTR, for the base substitution analysis. Of these base positions, 217 were categorized under stem region, in which mutations were observed in 58 positions. On the other hand, 143 positions were categorized under the loop region, in which mutations were observed in 45 positions. Therefore, per position, the rate of mutation in the loop regions was found to be 0.315, whereas the rate was only 0.267 in the stem region. The observed difference in mutation counts between the stem and loop regions was found to be statistically significant (p- < 0.01). In addition, the most frequent base substitution C → U in the stem region was compared with the loop region (Fig. 4). The rate of base substitution C → U in the stem region (0.340), was found to be less than half of the rate in the loop region (0.690). This proportionately lower mutation rate in stem is in concordance with the notion, that the mutations in the stems destabilize RNA secondary structures and, therefore, they are counter-selected.

Comparative analysis of transition and transversion in the stem and loop regions of RNA secondary structure of SARS-CoV-2

Twelve mutation frequencies considering all the stem-loop structures of the 5′-UTR and 3′-UTR were calculated and presented in Table 1. The mutation frequency values are shown in Fig. 5. Among the four transitions, C → U exhibited the highest frequency (0.364), which was more than the sum of the remaining three transitions U → C (0.132) G → A (0.126), and A → G (0.080). G → U (0.299) was the most frequent among the transversion mutations. This transversion value was even more than the transitions U → C, G → A, and A → G. This higher transversion frequency of G → U is consistent with earlier reported mutation patterns across different functional regions in SARS-CoV-2 genome. C → A (0.061) transversion was the next, followed by similar frequencies of A → U (0.050) and U → A (0.047), and G → C (0.046). Whereas, A → C (0.040), U → G (0.038), and C → G (0.024) were among the least frequent transversions. In total, transition and transversion frequencies were 0.704 and 0.605, respectively, resulting in a ti/tv ratio of 1.164, which was in accordance with the expected mutation pattern of the whole genome of the virus.

The transition to transversion ratio in the stem is higher than that in the loop region of RNA secondary structure of SARS-CoV-2

Considering the differential impact of base substitutions in the stem and loop regions, base substitutions were calculated separately in the stem and loop regions (Fig. 5). Though the size of the stem region was larger than the loop region, the number of substitutions in the stem region was lower than the loop region. The base substitution rate in the stem region was 0.28, whereas the same was 0.36 in the loop region. This lower rate of mutations in the stem region suggested that the stem region is more conserved in comparison to the loop region. In general, transition mutations were more frequent than the transversions in stem regions as well as in loop regions (Fig. 5). Among the transitions in the stem region, C → U was the most frequent substitution with a frequency of 0.294, followed by U → C, G → A, and A → G with frequencies of 0.197, 0.121, and 0.056 respectively. G → U substitution was the most frequent transversion, with a frequency of 0.224, that was comparable with the C → U transition. Among the substitutions in the loop region, the C → U transition had the highest frequency (0.484), followed by the transversion G → U (0.448). The two other transitions G → A and A → G were third and fourth in order of frequencies with values of 0.138 and 0.109, respectively. Interestingly, the transition U → C frequency was very low (0.044) compared to the other three transitions in the loop region. The other transversion frequency values were within the range of 0.103 (G → C) to 0.022 (U → G).

Considering these mutation frequencies, transition to transversion ratio (ti/tv) in the stem and loop regions was calculated (Fig. 6). In the stem region, total transition and transversion frequency were found to be 0.667 and 0.453, respectively, resulting ti/tv values of 1.472. In the loop region, total transition and transversion frequency were 0.775 and 0.888, respectively, resulting in ti/tv value of 0.872. The higher ti/tv value in the stem region, compared to the loop regions, was an outcome of the lower frequency of the transversion in the stem regions than the loop regions. This result suggested a deleterious effect of transversions in the stem region implicating that transversions might influence secondary structure of RNA.

In order to obtain statistical support for the higher ti/tv ratio in the stem regions, a Spearman rank correlation study between the twelve substitution frequencies in the stem and loop regions was done. When all the 12 mutation frequencies were considered, the Spearman rank correlation coefficient (ρ) was 0.510, suggesting that the order of the frequency values in stem and loop were similar. Yet, the frequency of U → C transition in the stem and loop region ranked the third and ninth respectively, suggesting U → C transition was fairly accommodated in the stem region without distorting the secondary structure. In contrast, the frequency of G → C and C → G transversions in the stem region exhibited the lowest values, whereas in the loop region, they displayed the fifth and the eighth highest values, suggesting that the avoidance of G → C and C → G transversions was stronger in the stem region in comparison to the loop region. However, higher frequencies of G → U transversion were observed both in stem and loop regions.

Analysis of base substitutions in the RNA secondary structure of Delta and Omicron variants of SARS-CoV-2

The base substitution analysis was extended to the stem and loop motifs of the UTRs of Delta and Omicron variants of SARS-CoV-2 (Fig. S4 ). In general, mutation frequency in the loop region was found to be with higher (2.537) in comparison to the stem region (2.131), indicating the differential role of secondary structure on base substitution. However, transition to transversion ratio (ti/tv) values in the stem and loop region were found to be similar. In the stem region, ti value was 1.135 and tv value was 0.996, resulting ti/tv value of 1.140. In the loop region, the ti value was 1.348 and, the tv value was 1.189 resulting ti/tv of 1.133. For understanding the similar ti/tv ratios in the stem and the loop regions, a further Spearman rank correlation study between the twelve substitution frequencies was carried out. When all the 12 mutation frequencies were considered, Spearman rank correlation coefficient (ρ) was found to be 0.79 as expected. The C → T and G → A transitions and G → T transversion were with top three ranks in both stem and loop region. However, the notable differences in rank values of the base substitutions in stem and loop regions were as follows. The frequency of transitions A → G, and transversions G → C and C → G were with higher rank in the loop region compared to the stem region. In contrast, T → A and A → C transversions were with lower rank in the loop region compared to the stem region (Fig. S4).