3 Results

In this study, both the tissue-specificity and expression volume restrictions were applied to facilitate identification of the most prominent interactions between host miRNAs and the viral genomes. Firstly, the scope of the study was limited to the lungs and colon tissues as the two main sites of SARS-CoV-2 replication in the human body [15]. Secondly, only those miRNA species were analyzed that represent the top 95% of the miRNAome for their tissue type [37], with the following rationale. In COVID-19, each infected cell expresses from 1,500 to 15,000 copies of the virus [38] and from 100,000 to 200,000 miRNA molecules [39]. To visualize the kinetics of their interactions, one should take into account that each miRNA species, in addition to the virus, binds to 25-50 host mRNAs species, each expressed as 50 to 100 copies per cell, thus, and therefore, being supplied with more than 1,250 endogenous target molecules [21]. Because of that, only highly represented miRNA species are likely to exert a significant impact on the evolution of the viral genomes. For lung tissue, there were 32 miRNAs (Figure 2 of Supplementary materials), and for the intestines, 40 miRNAs (Figure 3) that fulfilled the abovementioned criteria. Twenty-one of these miRNAs were found in both of these tissues, and 19 were expressed in either one or another tissue. In lung tissue, lung-specific miRNAs accounted for 81% of all miRNAs, while in the colon only 56% of expressed miRNAs were colon-specific. In lungs, the most prevalent miRNA was miR-143-3p, followed by miR-21-5p, miR-22-3p and miR-30a-5p. In colon, miRNA expression profile had a predominance of let-7b-5p, miR-92a-3p and miR-200c-3p. From the GISAD database, a total of 3,755 SARS-CoV-2 sequences were retrieved. All of them were tagged by location within the same metropolitan area (Berlin, Germany) and a time stamp between the 1st January 2020 and the 20thDecember 2022 (Figure 1), see “On choosing a city of the study” section of Supplementary materials for more details). All selected sequences were divided into two groups:
“Ancestral variants”, which included “Alpha” (B.1.1.7-like, 759 variants), “Beta” (B.1.351-like, 4 variants), Delta (B.1.617.2-like, 841 variants; AY.4.2like, 12 variants) as well as 7 other strains (Table 1 of Supplementary materials). All of these sequences were collected in January 2022 and earlier; “Ancestral variants” included “Alpha” (B.1.1.7-like, 759 variants), “Beta” (B.1.351-like, 4 variants), Delta (B.1.617.2-like, 841 variants; AY.4.2like, 12 variants) and 7 other strains (Table 1 of Supplementary materials) and belonged to the period until January 2022;
“Omicrons”, which included 8 varieties of “Omicron” strain. All of these were collected in December 2021 or later.
Figure 1: SARS-CoV-2 variants obtained from GISAID portal according to Berlin, Germany geotag.
To assess the evolutionary pressure exerted by colon and lung-specific miR-
NAs, the expression-normalized amounts of binding motifs numbers within 29,903 nucleotides of viral SARS-CoV-2 RNA were analyzed (see “Materials and methods” for more details). In this way, the regulatory contribution of each miRNA was taken separately into account. As compared to “Ancestral variants”, Omicron variants contain strikingly lesser amounts of binding sites for either colon miRNAs (Mann-Whitney U test, p ≤ 4.44 · 10−254, Figure 2, A) or lung miRNAs (Mann-Whitney U test, p ≤ 1.07 · 10−146, Figure 1 of Supplementary materials). By reshuffling the list of miRNAs with a predefined expression repertoire specific to the highly expressed miRNAs, the proportion of random shuffles for which the drop in the number of binding motifs was greater than the drop detected for the highly expressed miRNAs was calculated (the detailed description of the procedure can be found in “Materials and Methods” section). This proportion turned to be lower than 5% for colon miRNAs (p ≤ 2.50 · 10−2), but not for lung miRNAs (p ≥ 7.79 · 10−2). When the significance of the miRNA influence was tested for other tissue compartments, including the prostate, the breast, the bladder, and the liver, the decrease in the number of observed motifs was not significant for either of these tissues (Table 2 of Supplementary materials). Further study of the colon-specific dropdown showed that it was mainly observed within Omicron brunch (BA.2-like) (see Figure 2, A). However, even after exclusion of Omicron (BA.2-like) from “Omicrons” group, the indicated trend remained strong (Mann-Whitney U test, p ≤ 2.57 · 10−69). Moreover, the weighted number of miRNA binding motifs decreased steadily according to time stamps (Spearman correlation ≤ −0.34, p ≤ 7.37 · 10−80, Figure 2, B), thus, pointing at continuous evolution of “Omicron” sequences within the intestinal niche.
The difference was calculated between the number of binding motifs in the
“Wuhan” and the subsequent “Omicrons” (Figure 3). It turned out that the loss of the binding sites was primarily due to diminished binding to hsa-let-7b-5p, hsa-let7a-5p and hsa-miR-93-5p. Notably, for hsa-miR-194-5p and hsa-miR-29a-3p, the amounts of binding sites in the “Omicrons” have increased. The pairwise alignment of SARS-CoV-2 RNAs with reference “Wuhan” variant allowed us to compare the distributions of the miRNA seed regions within each variant of the virus. Further, the binding distributions for individual miRNA were calculated. It was found that the contribution of each miRNA in those distributions was proportional to its expression in colon tissues. Finally, the distributions were averaged over “Ancestral variants” and different groups of “Omicrons”, either including or excluding Omicron (BA.2-like) variants, and analyzed. As seen from Figure 4, the decrease in amounts of the complementary seeds observed in Omicron BA.2-like variants was primarily due to a loss of hsa-let-7b5p binding motifs located within NSP4. Indeed, C9861T nucleotide mutations occurred mostly in Omicron BA.2-like sequences (872 variants of 887), while the 9861st position in all “Ancestral variants” sequences retained the cytosine. Interestingly, that transition was silent as it has not lead to amino acid mutation: both CTC and CTT triplets encode for lysine.
Notably, a comparison of “Ancestral variants” and other “Omicrons” (Figure
3 of Supplementary materials) showed that the overall dropdown of miRNA binding sites was not due to any specific mutation within a particular miRNA binding motif. On the contrary, a gradual decrease in the number of binding motifs was observed over time (Figure 2, B). To access the contribution of virus RNA coding regions to the observed dropdown, the weighted (by miRNA expression) amounts of miRNA binding motifs per coding region were calculated for “Ancestral variants” and for “Omicrons” with or without Omicron (BA.2-like) sequences separately. By comparing these amounts within “Ancestral variants” with “Omicrons” using Mann-Whitney U test (Table 1), it turned out that Spike and NSP15 regions were changed the most. In addition, the sequences classified as Omicron (BA.2-like) displayed an additional region with accumulated changes, the one encoding for NSP4, while in other “Omicrons” the third most divergent location was the region of NSP12. Figure 2: Weighted number (by miRNA expression in colon tissues) of binding motifs per 100 miRNA and 29,903 nucleotides A: grouped by “Ancestral variants” and “Omicrons”, B: distributed over time. Spearman correlation between the weighted number and number of days since the start of pandemics is ≤ −0.34, (p ≤ 7.37 · 10−80). In this analysis, Omicron (BA.2 like) variants were excluded. Figure 3: Amounts of the binding motifs in analyzed SARS-CoV-2 variants and reference “Wuhan” variant differ. Here, miRNAs are sorted by the levels of their expression in the intestine. The share of miRNome represented by particular miRNA could be found adjacent to miRNA ID. Figure 4: Distributions of miRNA binding motifs along to “Wuhan” RNA sequences and averaged by A: “Ancestral variants”,B: difference between “Ancestral variants” and Omicron (BA.2-like) variants, and C: Omicron (BA.2like) variants.