Using the technique from our earlier research (7), samples were taken from 100 symptomatic patients with a high temperature and an intermittent cough. The samples were then immediately put in a viral transport medium. Patients suspected of having COVID-19 infection were referred to the diagnostics laboratories of the private BiogenMed COVID-19 testing laboratory in Tashkent, Republic of Uzbekistan. After SARS-CoV-2 testing in the laboratory, specimens of biological origin were randomly collected from PCR-positive individuals.

The scientific inquiry was approved by the Ethics Committee of the Ministry of Health of the Republic of Uzbekistan (#6/20-1582). All experiments were conducted in conformity with the applicable standards and legislation. Samples were renumbered and de-identified so that no one, not even researchers, knew who the patients were. Only anonymized data, including age and biological sex, were preserved for reporting purposes. All patients in the sample collection provided verbal consent for voluntary participation. We stated to all participants that the collected samples would be used for a sequencing project exclusively, without disclosing their identities or disturbing them in the future. In this study, verbal consent was preferred over written consent since patients were hesitant to sign any written document due to their concern about COVID-19 infection at the time and a lack of knowledge of the genome sequencing investigation. No minors were involved in the sample collection. There was no need to do so because the sequencing experiment was non-invasive, participants were not subjected to further downstream clinical procedures, and the identity of samples was anonymized in this investigation.

RNA extraction, Real-Time PCR, and sample evaluation were carried out after our previous work (7). Among all evaluated patients, 48 PCR-positive samples (28 females and 20 men) were randomly chosen for further research.

SARS-CoV-2 sequencing and data analysis were carried out as in our previous study (8). We used coverage analysis (v5.12.0.0) and variant caller (v5.12.0.4) while estimating the probability of the data belongings. Variant callers were filtered to retain only dependable variants and eliminate variants with reading depths of less than 1000 and ion torrent quality scores of less than 400 (Table 2 and Supplementary Table 1). We utilized the Maximum Likelihood Tree in the Molecular Evolutionary Genetics Analysis (MEGA)¹ program to cluster samples for the filtered variants. The consensus for each SARS-CoV-2 genome sequence was then submitted to the NCBI under accession numbers GI:2085183815 to GI:2085183892 (or MZ892621.1 to MZ892627.1; NCBI database; Table 1) and to GISAID² under accession numbers EPI_ISL_3668627 to EPI_ISL_3673673 (available for registered users).

More than a hundred patients suffering from high fever and sporadic cough symptoms were randomly selected in this work from a Tashkent-based private COVID-19 testing laboratory. Forty-eight PCR-positive samples were chosen for sequencing. Only 17 high-quality sequences were selected to be submitted to worldwide databases (twelve women and five men with an average age of 34, Table 1).

The results in this work show twice the higher infectivity level of women after random analysis. Twelve of seventeen patients were women among those subjected to SARS-CoV-2 infection during random selection. Thus, we expect greater sensitivity of women among Uzbekistani patients than men to SARS-CoV-2 infection. We did not observe the tendency of older people to be more sensitive to coronavirus infection. Among the studied patients, only four out of seventeen people were >50 years old after random selection. The Delta variant is notably more transmissible, estimated to have a 50–80% increase in transmissibility compared to the Alpha variant (13). In addition, the Delta variant is believed to have higher transmissibility compared to earlier strains of SARS-CoV-2, particularly among children and adolescents. According to a decision analytic model analyzing 106,866 confirmed COVID-19 infections, susceptibility to the Delta variant was significantly higher among the 10 to 15-year-old age group, showing a 1.92-fold increase compared to the pre-Delta variant (14). This corresponds with our results: seven patients were <30 years old after random selection.

There are numerous gaps in the consensus sequence as a result of 31 SARS-CoV-2 samples being disqualified from additional analysis due to inadequate sequencing coverage. The average number of mapped reads for the remaining 17 samples (12 women and five men) was 779 739 (Table 1). The target reads were 99.61%, and the mean read depth was 3814. The average coverage consistency in the chosen samples was 87.8% (Table 2).

Using the maximum likelihood method, the MEGA X program generated a phylogenetic tree based on the 17 viral sequences of SARS-CoV-2 taken from samples of COVID-19 patients in Tashkent city. The Wuhan reference strain NC 045512.2 was the root of this maximum likelihood tree, which was constructed using 134 mutations found in the 17 sequences (Supplementary Table 1 and Figure 1). Bootstrap values are illustrated.

Variant Caller produced all evaluated mutations in these selected samples. The number of mutations ranged from 36 (samples 38,47, and 36) to 45 (sample 48), altogether unique and shared mutations (Figure 1). Compared to the NC_045512.2 reference genome, most viral genomes represented between 39 and 44 mutations. These included thirty-nine mutations in samples 39 and 40, forty mutations in samples 10, 33, 34, and 43, forty-one mutations in samples 37, 41, 44, and 45, and forty-three mutations in only samples 35 and forty-four mutations in samples 42 and 46. The mutations mentioned above were calculated together with shared and unique mutations. Samples 38 and 46 do not carry unique mutations (Figure 1).

The total number of shared mutations and unique mutations is calculated for each sample.

The most commonly occurring nucleotide substitution found was cytosine to thymine (51/134 mutations), then guanine to thymine (28/134 mutations), thymine to cytosine and guanine to adenine (12/134 mutations), and then adenine to guanine (10/134 mutations). The Supplementary Table 1 shows the breakdown of these mutations. Every mutation had homozygosity. In the gene encoding the S protein we identified one shared disruptive inframe deletion (c.467_472delAGTTCA), two unique (c.*4300A>G, c.*4352G>T), two shared (c.*4308G>T, c.*4358G>T) downstream gene variant mutations, four unique (c.646C>T, c.2533G>T, c.3761G>T and c.3790G>T) and nine shared missense variant mutations, four unique (c.936C>T, c.1995C>T, c.2259G>A and c.2388T>C) and one shared (c.3183C>T) synonymous variant mutations. In the M (matrix) region one unique (c.100C>T and c.463C>T) and one shared (c.245T>C) missense mutation was observed; in N (nucleocapsid) region one shared (c.-3delA) upstream gene variant mutation, three unique (c.200C>T, c.518C>T and c.1085C>T) and six shared (c.188A>G, c.608G>T, c.643G>T, c.1129G>T, c.1152G>C and c.1154G>A) missense mutations and four unique (c.105G>A, c.123G>A, c.987G>A and c.1080C>T) synonymous mutations were identified. No mutations were found in the E (envelope) region (Supplementary Tables 1, 2). The nucleotide substitution C>T was established to enhance the viral adaptation in hosts (15). The reason was linked with the occurrence of amino acids of a hydrophobic nature that could contribute to better cell penetration. This change was reported as evidence of host-dependence and approach to combat coronavirus (16, 17).

In addition, our analysis revealed one frameshift mutation, 41 missense mutations, 27 synonymous mutations, and two upstream gene mutations within the ORF1ab region. Within the ORF3a region, one synonymous mutation and ten missense mutations were identified. In the ORF6 region, one missense mutation was observed, while the ORF7a region displayed six missense mutations and one synonymous mutation. Additionally, the ORF8a region exhibited one missense mutation, one upstream gene mutation, and one conservative insertion-deletion mutation. Finally, one synonymous mutation (c.57C>T) was found in the ORF10 region. We found 134 mutations in total, 65 shared and 69 unique, representing one frameshift mutation, one conservative and disruptive inframe deletion, four upstream region mutations, four downstream region mutations, 39 synonymous mutations, and 84 missense mutations (Supplementary Table 1).

We wanted to examine the significant alterations in all sequences to understand what distinguished our cases from others worldwide. The software MEGA X was employed to construct a phylogenetic tree utilizing the maximum likelihood method, incorporating the analysis of 17 viral sequences (Figure 1). The variant caller identified one hundred thirty-four mutations from the abovementioned sequences (v.5.12.0.4) (Supplementary Table 1).

Most of the shared mutations with more sequences, for example, eighteen mutations (C3037T, G15451A, C16466T, C21618G, GAGTTCA22028G, T22917G, C22995A, A23403G, C23604G, G24410A, C25469T, AGATTTC28247A, T26767C, TA28270T, A28461G, G28881T, G29402T, and G29742T) with 17 sequences and four mutations (G210T, C241T, C14408T, T27638C, C27752T) with 16 samples, seven mutations (C6402T, C8986T, G9053T, A11201G, A11332G, C19220T and G28916T) with 12 sequences and four mutations (C14925T, A21137G, A24110C and A25439C) with six sequences, accordingly (Supplementary Table 1).

These variants have been grouped into the USA, Indian, and England variants. All are considered Delta strains (B.1.617.2) in the phylogenetic tree provided by genomedetective.com. Interestingly, Beta variant samples are located in the Neighbor Cluster in the tree, whereas Alpha, Gamma, and Omicron variant samples were found in the further branches (Figure 3).

This phylogenetic tree was generated in www.genomedetective.com using our sequenced samples.

*B.1.1.7 I20 belongs to the Alpha variant. B.1.351 20H belongs to the Beta variant. P.1 20J belongs to the Gamma variant. B.1.617.2 21A, 21I, 21J belong to the Delta variant. BA.1 21K, BA.2 21L, BA. 1/BA.2, BA.4 22A, and BA.5 22B belong to Omicron variant of SARS-CoV-2.

Our findings show that the whole-genome sequences obtained from the 17 symptomatic COVID-19 patients represent significant nucleotide diversity (Table 1 and Supplementary Tables 1, 2). Based on the similarity of mutations observed in the Delta strain, the 17 sequenced samples have been clustered into a predominant clade of SARS-CoV-2 within the GISAID public database, designated as clade G (or GK subclade). The G clade has been known for S-protein mutations such as D614G, L452R, P681R, and T478K [11]. These mutations were also found in all 17 samples using the viralvar.org website. However, some other mutations in this region were found in each sample. For example, 34 S-protein mutations were identified in sample 39, while 23 mutations were found in sample 47 (Supplementary Table 3).

In accordance with the findings of Cherian et al. (18), the recently identified lineages B.1.617.1 and B.1.617.2 were prominently represented within the phylogenetic tree. These lineages were characterized by signature mutations, including L452R, T478K, E484Q, D614G, and P681R, within the spike protein, notably within the receptor-binding domain (RBD) (18). Those mutations, except E484Q, were also shared with all 17 samples in our result. In addition, mutation T19R was shared with all samples, whereas G142D (in 15 samples), D950N (in 12 samples), I850L (in 6 samples), T95I (in 5 samples), A222V and F797C (in 2 samples) were found (Supplementary Table 3).

The structural analysis of the receptor-binding domain (RBD) mutations, specifically L452R, T478K, and E484Q, indicated a potential enhancement in ACE2 binding affinity. Conversely, the presence of P681R within the furin cleavage site suggested a possible acceleration in the rate of S1-S2 cleavage, thereby facilitating heightened transmissibility. Moreover, the RBD mutations L452R and E484Q were found to diminish the binding affinity to specific monoclonal antibodies (mAbs), potentially impacting their neutralization efficacy (18). Interestingly, the E484Q mutation did not appear in our sequenced samples.

Most of the unique amino acid changes were found in the S-spike region in Sample 39 (26 mutations; F106L, G107V, T108L, T109L, L110stop, D111I, S112R, K113R, T114P, Q115S, S116P, L117Y, I119L, V120L, N121I, N122T, A123L, T124L, N125M, V126L, V127L, I128L, V130S, C131V, E132N, and Q134N) and in Sample 47 (16 mutations; A260V, G261W, A262C, A264S, Y265L, Y266L, V267C, Y269L, L270S, A263C, L277I, Q271S, P272T, R273stop, T274D, and L276S). In contrast, five AA changes were found in Sample 33 (Supplementary Table 3). The virus genome must be sequenced to identify SARS-CoV-2 strains and investigate local and worldwide dissemination. In addition, the comprehensive sequencing of the entire viral genome associated with infection is a valuable tool in elucidating outbreak dynamics. This includes investigating variations in pandemic size over time, spatial dissemination patterns, and transmission routes, as emphasized in the WHO COVID Report of 2021. Indeed, this investigation revealed that many diseases originated in the United States, India, and England (GK subclade) (Figure 2). The global distribution of the virus was caused by worldwide travel.

The disparities observed in virus structure among the samples illustrate the evolution of viral morphology over time in response to varying environmental conditions. The genomic sequence data derived from these 17 samples, subsequently deposited in global databases such as NCBI and GISAID, hold significant promise for public health and research entities in elucidating patterns of disease dissemination within the region. Moreover, these findings are anticipated to inform the development of diagnostic assays, therapeutic interventions, and vaccine candidates. Several institutions, including our center, have initiated efforts toward developing a national vaccine based on the SARS-CoV-2 genome sequences identified in symptomatic patients from Uzbekistan. The sequence information presented herein enriches the COVID sequence database (GISAID) with updated sequence data and mutational profiles specific to our geographical region. This contribution holds significance for forthcoming molecular epidemiological investigations and evolutionary phylogenetic studies conducted by health and scientific organizations.

The structural map of SARS-CoV-2 proteins was developed based on the one by Jamison et al. (19). The proteins’ 3D structures were developed by modeling their sequences in Swissmodel.³ The ones with no crystals obtained yet were not modeled in the database; thus, we did not include their models.

The mutations observed in this work belonged to both structural and non-structural proteins (Figure 3). Among non-structural proteins, we did not determine changes in NSP1, NSP7, NSP8, NSP9, NSP10, and NSP11. The number of mutations observed in this work coincides with proteins’ molecular parameters; the highest number of mutations were determined in NSP3 and spike protein. NSP3, the largest coronavirus protein, is one of the essential proteins involved in replication/transcription processes (20). NSP3 is one of the well-established proteases in SARS; NSP3 is responsible for papain-like functions (21) and NSP5 is known as 3CL (chymotrypsin-like main) protease (22).

Compared to protein length, the highest numbers of mutations in the genome were observed in N protein (9 mutations) and ORF3a (10 mutations). Even though the highest numbers (12 and 13) of mutations were found in NSP3 and spike proteins, respectively, their frequencies were intermediate relative to protein length. Compared to previous variants, the delta variant has higher transmissibility, and the efficacy of vaccines was lower against it (9). We needed to carry out whole genome reads to understand the genome changes better and fight those new variants. Our results demonstrated twice higher missense mutations than synonymous ones.

The high pathogenic effects of coronavirus during the pandemic were linked to many features. The affinity receptor-binding domain (RBD) of spike protein to human angiotensin-converting enzyme 2 (ACE2) was one of the most cited points contributing to the deadly effects of the virus (23). Mutations in the RBD in various virus strains determined its impact on ACE2, leading to reduced or enhanced binding affinity (24). In this work, we observed six mutations (plus one deletion) in all studied 17 patients (Table 1). Only four changes were determined in one patient. Two mutations were determined in 2 and 6 patients (Table 3).

The changes D614G, L452R, P681R, and T478K were earlier detected in other works (25, 26). Among these changes, the ones causing the changes in electrostatic interactions on the surface can be considered the most significant. Surface charge is regarded as one of the points that could contribute to spike protein binding efficacy based on electrostatic interaction that might facilitate conformational changes (27, 28). Among the patients studied in this work, we observed several changes that did not affect the affinity of the RBD of the spike protein to ACE2. Most of the changes either belonged to sites other than RBD or were not included in RBD-ACE2 interactions. Three of six mutations detected in all patients in this work were changes to either Arg or Lys. The significant change belonged to pLeu452Arg, where the Leu at 452 position changed to Arg. Aggarwal et al. (29) established enhanced binding affinity resulting from the pLeu452Arg mutation (29). This mutation was earlier found in Indian patients (18, 28). The location of Leu at the 452 position is described in Figure 4, together with an indication of spike protein RBD and ACE2.

Another significant amino acid change was T478K, the one observed in all patients, leading to the substitution of hydrophobic threonine (T) to charged lysine (K), which was established to enhance spike protein surface electrostatic potential that might affect spike RBD affinity with human ACE-2 (30).

Studying the consequences of patients’ health conditions if infected with various mutated/unmutated SARS-CoV-2 would result in a deeper understanding of the outcomes of multiple mutations. However, we did not study any clinical consequences due to anonymized samples ruled out by the ethical committee. In such cases, using computer software is another reliable approach (31). We will use computational analysis, such as molecular docking or molecular dynamics simulations, to predict the consequences. Because of the large number of mutations in this work, we planned to conduct computational calculations as separate research for several mutations. For now, our study is limited to the outcomes and discussions obtained.