Comparative analysis of SARS-CoV-2 quasispecies in the upper and lower respiratory tract shows an ongoing evolution in the spike cleavage site.

Studies are needed to better understand the genomic evolution of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). This study aimed to describe viral quasispecies population of upper and lower respiratory tract by next-generation sequencing in patients admitted to intensive care unit. A deep sequencing of the S gene of SARS-CoV-2 from 109 clinical specimens, sampled from the upper respiratory tract (URT) and lower respiratory tract (LRT) of 77 patients was performed. A higher incidence of non-synonymous mutations and indels was observed in the LRT among minority variants. This might be explained by the ability of the virus to invade cells without interacting with ACE2 (e.g. exploiting macrophage phagocytosis). Minority variants are highly concentrated around the gene portion encoding for the Spike cleavage site, with a higher incidence in the URT; four mutations are highly recurring among samples and were found associated with the URT. Interestingly, 55.8% of minority variants detected in this locus were T>G and G>T transversions. Results from this study evidenced the presence of selective pressure and suggest that an evolutionary process is still ongoing in one of the crucial sites of spike protein associated with the spillover to humans.

Introduction

Since coronavirus disease 2019 (COVID-19) was initially reported in China on 30th December 2019 (Dong et al., 2020; Wu et al., 2020), SARS-CoV-2 has been spreading worldwide. As of 8th March 2022, there have been 447 million confirmed infections and more than 6 million deaths have been reported worldwide (Dong et al., 2020, https://ourworldindata.org/coronavirus). The origin of SARS-CoV-2 is still debated but a hypothesis suggests the Malayan pangolins (Manis javanica) to be the possible intermediate host for the virus and recombination signals between pangolin, bat and human coronavirus sequences have been identified (Lam et al., 2020; Wong et al., 2020; Wu et al., 2020; Xiao et al., 2020). In fact, the SARS-CoV-2 genome sequence showed a high percentage of genomic identity (around 96%) with BatCoV-RaTG13 virus as well as (around 88%) with two other SARS-like bat viruses (Bat-SL-CoV-ZC45 and Bat-SL-CoV-ZXC21) (Zhou et al., 2020). Instead, the comparison with the SARS-CoV genome sequences showed an overall lower identity, approximately 79.6% (Gralinski and Menachery 2020; Zhou et al., 2020). Similar identity scores were observed when comparison analyzes were focused on the Spike (S) protein (around 75%) (Gralinski and Menachery 2020; Zhou et al., 2020). The S protein is the main determinant of viral tropism and is responsible for receptor binding and membrane fusion (Belouzard et al., 2012). For this reason, amino acid changes on this protein might have effects on infectivity, viral pathogenesis as well as transmissibility. It was initially reported that mutation D614G in the S protein, which has emerged and has become dominant, might have induced an enhancement of viral replication and viral fitness (Shi et al., 2020). Monitoring of emerging mutations, especially in the S protein, has been performed extensively with the establishment of Virus Evolution Expert Working Group (VEWG) by the WHO (World Health Organization WHO, 2021). The great effort on this concern is also highlighted by the huge number of SARS-CoV-2 sequences submitted on public repositories (e.g. GISAID). Initial report suggested a little viral diversity for SARS-CoV-2 (Karamitros et al., 2020; Simmonds, 2020), however, since December 2020 a process of positive selection with presumed advantages such as increased transmission rates has been documented for a series of variants of concern (VOCs) such as alpha, delta and omicron (Harvey et al., 2021; Tao et al., 2021; Volz et al., 2021). These VOCs have demonstrated a significant public health impact, with changes in the virus transmissibility and reduce the efficacy of vaccines (Harvey et al., 2021; Tao et al., 2021). Quasispecies is believed to be a strategy of virus evolution (Domingo and Perales, 2019) and although the viral kinetics of SARS-CoV-2 infection from the upper respiratory tract (URT) to the lower respiratory tract (LRT) have been gradually clarified more work is still needed to explore the inter-host and intra-host variations of SARS-CoV-2. Overall, studies aimed at investigating intra-host evolution or the dynamic of SARS-CoV-2 quasispecies have been mainly focused on samples collected from upper respiratory tract (URT) (Al Khatib et al., 2020; Jary et al., 2020; Pérez-Lago et al., 2021; Shen et al., 2020; Siqueira et al., 2020; Sun et al., 2021; To et al., 2020). In fact, the dynamic of SARS-CoV-2 population in the lower respiratory tract (LRT) of patients showing severe acute respiratory infections (SARIs) is poorly investigated. It would be important to elucidate the role of specific mutations in the progression of SARS-CoV-2 from upper to lower respiratory tract or identify specific mutational patterns associated with severe infection. In the present study, high-depth next-generation sequencing (NGS) of the S gene was performed in a set of 109 respiratory samples from URT (n = 58) and LRT (n = 51) in order to: (i) evaluate the genetic diversity in two different body compartments; (ii) identify minority variants potentially associated with progression from upper to lower respiratory tract on paired samples from patients admitted to intensive care unit (ICU) with severe infection.

Material and methods

Patients and samples

A total of 109 clinical specimens from URT (nasopharyngeal swabs; NPS) and LRT (bronchoalveolar lavage; BAL or broncho aspirate; Brasp) were collected and analyzed from 77 COVID-19-positive patients (Appendix Table S1). Respiratory samples were prospectively collected from patients hospitalized at ICUwith severe to critical COVID-19 disease and from patients with mild symptoms not requiring hospitalization according to the WHO clinical management of COVID-19 guide (World Health Organization WHOa, 2020). Among patients admitted to ICU with severe to critical COVID-19 infection, whenever possible, paired URT and LRT samples were collected. All specimens were collected from late February 2020 to January 2021 at the Microbiology and Virology Department of Fondazione IRCCS Policlinico San Matteo in Pavia as Regional Reference laboratory for COVID-19 diagnosis. The presence of SARS-CoV-2 RNA in respiratory specimens was assessed by using specific real-time reverse transcriptase–polymerase chain reaction (RT-PCR) targeting the RNA-dependent RNA polymerase and E genes, following the WHO guidelines and the Corman et colleagues’ protocols (Corman et al., 2020; World Health Organization WHOb, 2020). Quantification cycle (Cq) values according to MIQE guidelines (Bustin et al., 2009), were used as a semiquantitative measure of SARS-CoV-2 viral load. The sequence investigation of patient samples was approved by the Ethics Committee of our institution (P_20200085574).

S gene amplification and sequencing

Total RNA was extracted using the QIAamp Viral RNA Mini Kit according to the manufacturer's instructions using a starting volume of 400 µL elute in a final volume of 60 µL.

The extracted RNA was subjected to a one-step RT-PCR using the SuperScript IV One-Step RT-PCR System (Thermo Fisher Scientific, USA). Two strategies were adopted: a “long PCR” for the amplification of the entire S gene (near 4000 bp) using the primer pairs SARS-2-S-F3 (tatcttggcaaaccacgcgaacaa) and SARS-2-S-R3 (acccttggagagtgctagttgccatctc) or using a semi-nested approach with the following two primers pairs: the first step with SARS-2-S-F3 and SARS-2-R6 (ttctgcaccaagtgacatagtgtaggca), followed by a second step with SARS-2-F6 (tcaggatgttaactgcacagaagtcc) and SARS-2-S-R3 (complete list of primers and their position are reported on Appendix Table S2). The thermal profile for the retro transcription was 55 °C for 10 min, followed by the “long PCR” with an initial denaturation/RT inactivation step at 98 °C for 2 min, the amplification for 42 cycles including the initial denaturation (98 °C for 10 s), the annealing step at 60 °C for 10 s, and the extension at 72 °C for 3 min; a final extension of 5 min at 72 °C. The semi-nested PCR was performed using the Platinum SuperFi DNA Polymerase with 5 µl of the first-step DNA and the following thermal program: 98 °C for 30 s, the amplification for 40 cycles including the initial denaturation (98 °C for 10 s), the annealing step at 60 °C for 10 s, and the extension at 72 °C for 2 min and 30 s. The proper PCR products were purified with AMPure Beads, with elution in TE buffer. Enriched DNA samples were used to prepare sequencing libraries with the Nextera XT kit. Sequencing was performed on an Illumina MiSeq machine, aiming for ∼1,000,000 250 bp paired-end reads per sample.

In silico analysis of sequences

Quality control of sequencing reads was performed using the program FastQC (www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads in each sample were filtered and trimmed for quality using fastp (Chen et al., 2018). In addition, 28 bases were cut from each end of all reads, to remove sequences generated from the semi-nested PCR primers. Filtered reads were aligned to the S gene of the Wuhan-hu-1 reference genome (NC_045512.2) (Wu et al., 2020) using bowtie2 (Langmead and Salzberg, 2012). Haplotype sequences for each sample were obtained from the alignment SAM files, using the software CliqueSNV (Knyazev et al., 2018). The hedgehog algorithm was used on the most abundant haplotype of each sample to classify the SARS-CoV-2 strains using only the S protein sequences (O'Toole et al., 2022). In parallel, alignment data was processed with samtools (Li et al., 2009) and bam-readcount (https://github.com/genome/bam-readcount) in order to calculate the number of occurrences of each nucleotide and indel in all positions of the reference. Only nucleotides and indels with at least 1% prevalence were considered in the following analyzes.

Python and R scripting (scripts are available from github link) was used to extract and classify all mutations using the following algorithm:

For each position in all samples, the nucleotide or indel with the highest prevalence was called the “majority variant”.

All other bases with at least 1% prevalence were called “minority variants” (MVs).

The correlation between the presence of MVs and the respiratory tract district of sampling was tested for all positions of the gene using the Fisher exact test (p < 0.05).

MVs in each sample were counted and classified by mutation type (synonymous, non-synonymous or indel), gene sub-domain (Huang et al., 2020), and mutation pattern (from which majority base to which low prevalence base). The differential distribution of the number of MVs between URT and LRT samples was tested using the Wilcoxon rank sum test (p < 0.05). The test was repeated for all classes determined, weighting each count on the total number of MVs of the sample.

This allowed us to test their association with LRT or URT (Fisher exact test) and to measure the gene variability sampled in this work. Sequencing reads are available on the SRA database under BioProject ID PRJNA686083. The scripts generated to perform this project are available on GitHub at https://github.com/SteMIDIfactory/DeepSpike.

Statistical analyzes

Comparisons of Cq, number of minority variants and haplotypes were performed with Mann-Witney test for continuous unpaired variables and in paired respiratory samples with the Wilcoxon rank sum test for continuous paired variables. Correlations between two quantitative variables were measured by the Spearman correlation test. Fisher's exact test for categorical variables was used for analysing mutation frequencies between groups of patients. Descriptive statistics and linear regression lines were performed using Graph Pad Prism software (version 8.3.0).

Results

Patients population

A total of 77 patients were included in the study. Fifty-five (71.4%) of them had severe presentations and were admitted to the ICU, while 22 (28.6%) had mild infections not requiring hospitalization. Among 55 ICU patients, for 28 (50.9%) of them paired URT and LRT samples were available, for 19 (34.5%) patients only LRT samples were analyzed (two of them with two and one with three serially collected BAL, respectively) and for 8 (14.6%) patients only URT samples were available. Among 28 patients with paired URT and LRT samples, the LRT sample was collected at the same time for a great majority of paired samples (range -4 to 9 days of difference). From the 22 patients with mild disease only URT samples were collected and analyzed.

Dataset description

A total of 109 samples were included in the NGS analyzes, of which 58 (53.2%) were collected from the URT and 51 (46.8%) from the LRT. The DNA of the S gene was enriched with PCR methods and deep short-read sequencing was performed. A total of 69279628 reads were obtained from sequencing, with an average of 635593 reads for sample (range 56214-3346036). Reads were mapped to the SARS-CoV-2 reference strain and the average depth obtained was 6923x (range 1301x-7954x). Reads mapping allowed to extract, classify, and count all genomic variants present in the samples with a prevalence of at least 1% of the sequencing depth (corresponding to ∼10x depth in the samples with the worst read yield).

The hedgehog algorithm was used in order to classify the SARS-CoV-2 strains using only the S protein sequences. The great majority of SARS-CoV-2 included in this study belonged to A_1 lineage harbouring mainly the D614G change (84/109; 77.1%). Thirteen (11.9%) strains belonged to B.1.177_1 lineage (A222V, D614G), 7 (6.4%) to alpha VOC (lineage A_9; del69-70, N501Y, A570D, D614G, P681H, T716I, S982A, D1118H), 4 (3.7%) strains to lineage B.1_14 (D614G, D839Y) and 1 (0.9%) strain to lineage B.1.177.52_1 (A222V, D614G, P1162R) (Supplementary Material Table S2)

Viral load and correlation between Cq and intra-host variability

The median viral load measured as Cq value observed in URT samples (23.8; range 13.0–36.0) was similar to those observed in the LRT samples (22.6; range 12.0-34.9; p = 0.12) (Fig. 1 A). Similarly, among 28 paired samples no difference in the median Cq value was observed (p = 0.83). In order to describe the vial load variability on paired samples, we performed a correlation analysis and a plot on difference in Cq value between paired URT and LRT samples. In 13/28 (46.4%) paired samples the Ct value was lower in URT (higher viral load) as compared to LRT samples with a median ΔCq of 2.8 (range 0.04–11.0 ΔCq) (pink circles on Fig.1B), while in 15/28 (53.6%) paired samples, Ct value was lower (higher viral load) in LRT as compared to URT samples with a median ΔCq of 4.7 (range 0.1 to 9.3 ΔCq) (light blue circles on Fig. 1B). In addition, an overall correlation was observed among paired samples (rho=0.74, Fig. 1C).

Fig. 1

(A) Distribution of Cq in the URT and LRT samples. (B, C) Differences and correlation plots for Cq values in 28 paired URT and LRT samples. The statistic Spearman's correlation coefficient and linear regression R2 value are also reported.

In general, viral replication has been associated with the intra-host diversification of viral population. For this reason, we compared the Cq values, as expression of the viral load, and the number of MVs and haplotypes observed. No evidence of correlation between Cq and number of MVs in both URT (σ=0.16) and LRT samples (σ = 0.11) was observed (Fig. 2 A). Similar findings were observed in the comparison of Cq and the number of haplotypes in URT samples (σ = 0.01) as well as in the LRT (σ = -0.003) (Fig. 2B).

Fig. 2

Comparison of number of haplotypes and minority variants in upper vs lower respiratory tract samples. (A) Correlation between the number of minority variants and viral load expressed in cycle of quantification (Cq). (B) Correlation between the number of haplotypes and viral load expressed in Cq.

Haplotype and minority variant counts

The number of MVs is slightly higher in LRT (median 13.5; range 3–99) than in URT samples (median 8; range 0–263), but with no statistical significance (p = 0.07, Fig. 3 A). Conversely, the number of haplotypes identified was significatively greater in LRT samples (median 2; range 1–13), compared to URT ones (median 1; range 1-9; p = 0.02, Fig. 3B). The ratio of non-synonymous (dN) to synonymous (dS) substitutions (dN/dS) was calculated. The median value observed in LRT samples (median 2.65, range 0–11) was greater than those observed in UTR samples (median 1.81, range 0–8; p = 0.02), suggesting a higher positive selective pressure in the lung environment (Fig. 3C).

Fig. 3

Distribution of the number of (A) minority variants, (B) haplotypes and (C) dN/dS ratio in the URT and LRT samples. Values are represented as a boxplot with all points inscribed. (D) Distribution of the weighted incidence of synonymous, non-synonymous, and insertions and deletions (Indel) in the URT and the LRT samples. (E) Distribution of the weighted incidence of frameshifting insertions and deletions in the URT and the LRT. Values are weighted by dividing them by the total number of minority variants in the sample.

With a more in-depth analysis, the weighted incidence of synonymous mutations was higher in the URT samples compared to LRT samples (p = 0.01, Fig. 3D), while, although not significant, a greater number of indels was observed in the LRT samples as compared to URT samples, (p = 0.10, Fig. 3D). Lastly, no difference was observed in the number of non-synonymous mutations. The analysis was repeated considering only deletions and frameshifting insertions; both of them had greater incidence in the LRT, without significant difference (p = 0.12 and p = 0.13, Fig. 3E).

The weighted incidence of MVs was calculated in all regions of the gene corresponding to the functional and structural domains of the protein. We observed that MVs in the N-terminal Domain (NTD) are more likely to occur in the LRT (p < 0.001). On the other hand, mutations in the Fusion Peptide (FP, p = 0.006) and in Subdomain 2 (SD2, p < 0.001) are more common in the URT (Fig. 4 A). Furthermore, MVs in the region coding for the protein subunit S1 are more common in the LRT, while in subunit S2 we measured a significantly higher abundance in the URT (Fig. 4B). This result is expected as it is the reflection of the values observed in the functional domains.

Fig. 4

(A) Distribution of the weighted incidence of minority variants in the S gene subdomains (NTD: N-Terminal Domain; RBD: Receptor-Binding Domain; SD1: Structural Domain 1; SD2: Structural Domain 2; FP: Fusion Peptide; HR1: Heptad Repeat 1; HR2: Heptad Repeat 2; TM: TransMembrane domain) in the URT and LRT samples. (B) Distribution of the weighted incidence of minority variants in the two S gene subunits in the URT and the LRT (S1: Subunit 1; S2: Subunit 2). (C) Distribution in the URT and in the LRT of the weighted incidence of minority variants, classified by mutation patterns. Values are weighted by dividing them by the total number of minority variants in the sample.

Incidence of mutation patterns were tested as well, both counting the total events, and weighting them on the total number of MVs of each sample. Fig. 4C shows the weighted incidence of each mutation pattern in both respiratory tract compartments, while Appendix Fig. S3 shows the absolute one. Both analyzes show a higher presence of A>C and T>G mutations in the MVs of URT samples, while A>T, T>A, and T>C have a higher incidence in the LRT.

Lastly, we tested the correlation in all nucleotide positions between the presence of MVs (binary value) and the two respiratory tract districts (shown in Fig. 5 A). Fig. 5B, instead, shows the incidence of MVs along the gene. Mutations are equally distributed on the sequence, with the exception of the area around the cleavage site between the two subunits of the gene. Here mutation sites are highly concentrated both in the upper and in the lower respiratory tract samples, especially in the former, as already seen in Fig. 5A. In addition, the presence of mutations is associated with the respiratory district in nine codons, five correlated with URT, four with LRT (see Table 1 ). Four of the MVs associated with URT are located around the cleavage site of the two subunits. In this area, we detected a high concentration of MV transversions between T andG: from position 2000 to position 2150, 248 out of 643 (38.6%) mutations are T>G and 11 (17.3%) are G>T (total = 359; 55.8%)..

Fig. 5

Graphical distribution of changes along S protein gene. (A) Number of samples containing minority variants in each position. Two separate histograms are used for URT and LRT samples, which are indicated upside down for image clarity. (B) Correlation of the presence of minority variants with URT and LRT in each position of the gene. Bar height represents the log10 (p-value) of the Fisher exact test. In the middle, a scheme of the gene subdomains and subunits is used as separator. NTD: N-Terminal Domain; RBD: Receptor-Binding Domain; SD1: Structural Domain 1; SD2: Structural Domain 2; FP: Fusion Peptide; HR1: Heptad Repeat 1; HR2: Heptad Repeat 2; TM: TransMembrane domain; S1: Subunit 1; S2: Subunit 2.

Table 1

Description of the 9 Minority Variants (MVs) positions associated with either URT or LRT compared to the reference sequence (NC_045512.2). Global frequency is referred to the frequency of the mutation in the same amino acid in GISAID global database. Data are accessible at www.cov.lanl.gov .

Codon position	Amino acid	Reference	Mutation	Mutation type	Associated with
212	71	C(S)	A(Y)/T(F)	NotSyn	LRT
2055	685	T(R)	G(R)	Syn	URT
2058	686	T(S)	G(R)	NotSyn	URT
2060	687	T(V)	G(G)	NotSyn	URT
2100	700	T(G)	G(G)	Syn	URT
3005	1002	A(Q)	T(L) / -	NotSyn/Del	LRT
3483	1161	A(S)	C(S)	Syn	URT
3485	1162	C(P)	T(L)/G(R)	NotSyn	LRT
3596	1199	A (D)	G(G)/T(G)	NotSyn	LRT

Discussion

The evolution of SARS-CoV-2 was initially quite slow, when compared to other RNA viruses (van Dorp et al., 2020). Yet, its rapid global spread has allowed to record thousands of mutations in public databases; some of those were favourable and have emerged worldwide (Long et al., 2020). The emergence of VOCs was favored by more than 400 million of infections worldwide (Parra-Lucares et al., 2022) with a more significant number of mutations observed in S sequences as compared to other genomic regions (Yusof et al., 2021). Generally, mutations in viral structural such as S glycoprotein can play a crucial role in their virulence by possibly determining changes in their cellular tropism and the generation of antibody escape variants as reported for Delta and Omicron variants (Andrews et al., 2022; Dejnirattisai et al. 2022; Mlcochova et al., 2021). The emergence of these variants has been promoted by the viral quasispecies evolution and the severity of SARS-CoV-2 infection is driven by progression from URT to LRT (Ke et al., 2020). In this perspective, our study has investigated the genetic diversity in SARS-CoV-2 quasispecies focusing on structural S protein sequences in two body compartments in order to: (i) evaluate the genetic diversity in two different body compartments; (ii) identify minority variants potentially associated with progression from upper to lower respiratory tract on paired samples from patients admitted to intensive care unit with severe infection.

Among SARS-CoV-2 evolution, G614 variant in S protein has become worldwide predominant since April 2020 and was associated with an increased fitness advantage (Korber at al., 2020). This finding was also confirmed by other studies that compared D614 and G614 variants and found that G614 was associated with increased replication in human lung epithelial cells (Shi et al., 2020). On the contrary, G614 variant was not associated with an increased disease severity and its role in pathogenesis has yet to be elucidated (Long et al., 2020; Shi et al., 2020). More than 60% of patients included in this study had severe infections developing severe pneumonia and requiring oxygen therapy. Mutations associated with these symptoms were explored by obtaining sequences from the LRT samples. All SARS-CoV-2 sequences generated in this study harboured the G614 variant, while the original D614 variant was found neither among majority nor among minority variants. Thus, we do not report any evidence of the G614 variant favouring the severe presentation. Moreover, no evidence of mutations on S gene associated with progression from upper to lower respiratory tract emerged from our analysis. This result was also observed in a series of paired samples and is consistent with the finding, previously reported by Rueca et al. in a lower number of patients (Rueca et al., 2020). Our data are in keeping with Wylezich et colleagues reporting no evidence of compartment-specific pattern of mutations between different respiratory compartments (Wylezich et al., 2021). This finding suggests us that disease severity could be mainly determined by host factor such as comorbidities, age and absence of pre-existing immunity (Al Khatib et al., 2020). In addition, deletions have been observed as MVs in a few samples but with a lower frequency than that observed in a recent publication (Andrés et al. 2020).

Overall, the S gene sequences originated in our study showed a greater variability (number of haplotypes) in LRT as compared to URT samples. This difference was unrelated to the viral load (measured as Ct) since comparable Ct values were observed in URT and LRT samples. Indeed, no correlation between viral load and viral diversity was observed and this is consistent with the finding previously reported by Siqueira et colleagues, who investigated quasispecies variation in cancer patients (Siqueira et al., 2020). The difference in viral population between URT and LRT could be explained by the hypothesis of an independent replication in the two respiratory districts also suggested by Wölfel et al. (2020). In addition, a higher positive selective pressure in the lung environment has been observed as compared to upper respiratory tract (dN/dS>1). Similar observation was reported by Sun et al suggesting that diversifying of the quasispecies mutants indicated potential independent virus replication in different tissues or organs (Sun et al. 2021). The great variability in the LRT samples resulted also in an increased number of frameshifting deletion and insertions. The presence of deleterious mutations could indicate a loss of function of the S protein in a fraction of the viral population. This subpopulation might be maintained thanks to replication and cell invasions events that do not imply the ACE2 receptor (e.g. within syncytia or in macrophages after phagocytosis)(Abassi et al., 2020). However, this and other theories on how such mutations can influence viral replication in lung tissues might be elucidated and extensively investigated with additional studies.

Nucleotides changes arise during the virus replication and persistence, in particular A>G. These were shown to be related to the host editing mechanisms such as the APOBEC and ADARs (Carpenter et al., 2009; Di Giorgio et al., 2020). The A>G transition is caused by the deamination from Adenosine to Inosine (A>I) generated by the ADARs. Thus, the significant rate of T>C observed in the LRT both by Di Giorgio et al. (2020) and in this study agrees with the hypothesis that T>C in SARS-CoV-2 could also be related to the ADARs mechanism32. Although not associated with host editing mechanism, the T>G pattern in MVs is also of particular interest in this study, as they are the most prevalent pattern observed in the entire dataset and they were found associated with the URT.

Finally, the cleavage site between S1 and S2 of the S protein corresponds to one of the two genetic sites in which Andersen and colleagues found crucial mutations associated with the spillover to humans (Andersen et al., 2020). For this reason, we can hypothesize the presence of an evolutionary selective pressure in this site and that the driving force of this evolution might reside in the URT. Modifications (mutations or deletion) in the S1/S2 junction site has been associated with virus attenuation in hamsters (Wang et al., 2021). Alternatively, these observations might be explained by an absence of negative selection on random occurring mutations: this is a cleavage site; thus, SNPs affect the protein structure much less than in other sites. This last hypothesis explains the high density of mutations in the cleavage site (S1/S2 junction) both in URT and LRT samples. However, it does not explain the higher incidence in the URT and the presence of associated mutations. Moreover, other studies found specific low frequency mutations around the cleavage site, i.e. deletions that were associated with milder symptoms (Andrés et al., 2020; Lau et al., 2020). Whereas, in our study mutations are mainly T>G and G>T transversions in this site (359/643 MVs). Since such mutations are usually rare changes in nucleic acids, this observation underlines further the presence of a selective pressure. Interestingly, the T>G pattern was also identified as an inexplicable intra-host mutational signature in HPV (Zhu et al., 2020).

In conclusion, the results of the present study indicate that severe SARS-CoV-2 infections are not associated with a specific mutational pattern. However, a great variability was observed on viral population in LRT also associated with a positive selective pressure. How the difference may impact immune response escape, tissue tropism and pathogenicity is still to be elucidated. Moreover, we observed the evidence of possible ongoing evolution in one of the gene loci that were crucial for the spillover to humans. This highlights the importance of genomic surveillance to predict and avoid vaccine escape mutants.

CRediT authorship contribution statement

Stefano Gaiarsa: Writing – original draft, Software, Formal analysis, Data curation. Federica Giardina: Writing – original draft, Formal analysis, Investigation. Gherard Batisti Biffignandi: Methodology, Software, Formal analysis, Data curation. Guglielmo Ferrari: Investigation. Aurora Piazza: Validation, Investigation, Data curation. Monica Tallarita: Investigation. Federica Novazzi: Investigation. Claudio Bandi: Writing – review & editing. Stefania Paolucci: Resources, Visualization. Francesca Rovida: Investigation, Resources. Giulia Campanini: Writing – review & editing. Antonio Piralla: Conceptualization, Methodology, Writing – original draft, Project administration, Funding acquisition. Fausto Baldanti: Supervision, Project administration, Writing – review & editing.

Declaration of Competing of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper