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.