Identification of Delta Variant-Specific Mutations of SARS-CoV-2 Genomes Found in Uzbek Patients by Amplicon Sequencing

1. Introduction

Considerable research has been conducted on the SARS-Cov-2 genome since 2020. To date, 15,778,185 SARS-CoV-2 genome sequences have been shared via online data platforms such as GISAID [1], which helps track the new variants and mutations. The variants of SARS-CoV-2, the causative agent of COVID-19, continue to change during its travel by human-to-human transmission [2].

The World Health Organization (WHO) officially declared a pandemic on March 11, 2020 [3]. SARS-CoV-2 has scattered globally and is currently a significant problem. Europe and Latin America have suffered more than other nations, particularly in the early stages of the coronavirus epidemic. As of 12 July 2023, WHO had received reports of 767 972 961 confirmed COVID-19 cases, including 6 950 655 fatalities. A total of 13 462 024 421 vaccination doses have been given as of 8 July 2023. From 3 January 2020 to 12 July 2023, Uzbekistan reported 171 083 confirmed cases of COVID-19, with 1 016 fatalities. So far, 80 043 836 vaccination doses have been given as of 3 June 2023 [4].

The established nomenclature systems for naming and tracking SARS-CoV-2 genetic lineages through GISAID, Nextstrain, and Pango are currently being used by scientists in scientific research. The Technical Advisory Group on Virus Evolution was convened by WHO and recommended using letters of the Greek Alphabet, i.e., Alpha, Beta, Gamma, Delta, etc., of the virus strains to make it easier and more practical to be discussed by non-scientific audiences [3].

SARS-Cov-2 belongs to the β-coronavirus genus lineage B, which has large enveloped, positive-sense, single-stranded RNAs (30 kb) and encodes four structural proteins as well as other accessory or non-structural proteins (including the viral pp1a-pp1ab replicase, the 3C-like protease (3CLpro), the papain-like protease (PLpro), and the RNA-dependent RNA polymerase (RdRp) [5,6].

SARS-CoV2 whole genome sequencing is critical for discovering meaningful information on viral lineages, variants of interest, and variants of concern [7]. The Global Initiative on Sharing All Influenza Data (GISAID) [1], the National Center for Biotechnology Information (NCBI) [8], and the Virus Pathogen Database and Analysis Resource (ViPR) [9] have full and partial genome sequences of SARS-CoV-2 obtained from various countries and laboratories.

The worldwide phylogenetic tree of SARS-Cov-2 was constructed using sequenced virus genomes, and the sequences are divided into several clades. Our previous study examined the earliest clades and their characteristics [10]. Many viral variants were later discovered and uploaded to global databases. The B.1.617.2 variant was detected in India in December 2020 [11], and the World Health Organization classified it as VOC Delta on May 11, 2021, according to its increased transmissibility [12]. Delta variants with increased transmissibility have been linked to critical S-protein mutations such as D614G, L452R, P681R, and T478K [13]. The spike protein facilitates virus attachment to host cell-surface receptors and virus-cell membrane fusion. It is also the primary target for neutralizing antibodies produced after infection [14].

For the first time, an infectious epidemic was found in Uzbekistan in March 2020. The government decided to conduct public testing for COVID-19. This incident signaled the beginning of the pandemic in Uzbekistan, necessitating accurate identification and monitoring of coronavirus genotypes that were rapidly spreading throughout the population [10,15].

We, therefore, started the whole genome sequencing of COVID-19 samples with the primary objective of identifying virus genotypes spread in our territory as well as analyzing genomic diversity, types of mutations, and the emergence of new variations of SARS-CoV-2 [10,15]. In this study, we describe the following whole-genome sequence data from individuals in the Republic of Uzbekistan who were infected with COVID-19. We analyzed 134 mutations with nonsynonymous and synonymous types and successfully assembled the 17 high-quality sample genome sequences for coronavirus genotypes. In the global phylogenetic tree, comparative analysis utilizing SARS-CoV-2 genomes that are known as the Uzbekistan sample genomes into the GK subclade. The sequenced genomic data of coronavirus genotypes were used to develop the edible vaccine in Uzbekistan against coronavirus threats [16].

2. Materials and methods

2.1. Sample collection

Using nasopharyngeal and oropharyngeal swabs sticks (Huachenyang Technology, Shenzhen, China), samples were obtained from 100 symptomatic patients with high fever and intermittent cough and immediately placed 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. Following SARS-CoV-2 laboratory testing, biological samples were taken randomly from PCR-positive patients.

The Ethics Committee of the Ministry of Health of the Republic of Uzbekistan (#6/20-1582) approved the research investigation. 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. There were no minors 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 samples' identities were anonymized in this investigation.

A MagMax Viral/Pathogen kit was used to extract total RNA from the collected samples (Life Technologies Corporation, Austin, Texas, USA). Using a CFX Connect Real-Time PCR System, each sample was analyzed for the presence of SARS-CoV-2 (threshold - 0,050) (BioRad, Hercules, California, USA). RT-PCR was carried out in 40 ul final volumes. For cDNA synthesis at 35°C for 20 minutes, a first denaturation at 95°C for 5 minutes was followed by 50 cycles of 94°C for 15 seconds and 64°C for 20 seconds (set fluorescence measurement for Fam, Hex, Rox, and Cy5 channels at 64°C) and a SARS-CoV-2/SARS-CoV multiplex real-time PCR assay targeting the nucleocapsid (N), envelope (E), and spike region [10]. Among all evaluated patients, 48 PCR-positive samples (28 females and 20 men) were chosen at random for further research.

2.2. SARS-CoV-2 sequencing

As in our previous study, SARS-CoV-2 sequencing and data analysis were carried out [10]. Variant callers were filtered to remove variants with reading depths of less than 1000 and ion torrent quality scores of less than 400 to keep reliable variants only (Table 2 and Table S1). The filtered variants were used for sample clustering with Maximum Likelihood Tree in Molecular Evolutionary Genetics Analysis (MEGA, https://www.megasoftware.net) software. The consensus for each SARS-CoV-2 genome sequence was then submitted to the GISAID (www.gisaid.org) under the accession numbers of EPI_ISL_3668627 to EPI_ISL_3673673 (available for registered users) and NCBI under the accession numbers GI:2085183815 to GI:2085183892 (or MZ892621.1 to MZ892627.1; NCBI database; Table 1).

3. Results

3.1. Sample selection for sequencing

Over a hundred symptomatic patients with high fever and intermittent cough were chosen randomly from a private COVID-19 testing lab in Tashkent for this investigation. 48 PCR-positive samples were chosen for sequencing. Only 17 high-quality sequences were selected to be submitted to worldwide databases (five males and twelve women attended, with an average age of 34; Table 1).

Because of insufficient sequencing coverage, 31 SARS-CoV-2 samples were excluded from further analysis, resulting in many gaps in the consensus sequence. The remaining 17 samples (five males and twelve women, Table 1) had an average of 779 739 mapped reads per sample. The target reads were 99.61%, and the mean read depth was 3814. The average coverage consistency in chosen samples was 87.8%; (Table 2).

3.2. Identifying unique and reliable mutations in all high-quality sequenced samples

Using the maximum likelihood method, a phylogenetic tree was generated in the MEGA X program based on the 17 viral sequences of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in samples collected from COVID-19 patients in Tashkent city. This maximum likelihood tree was generated using 134 mutations discovered in the 17 sequences and was rooted by the Wuhan reference strain NC 045512.2 (see Table S1, 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 2). Most viral genomes contained between 39 and 44 mutations when compared with the NC_045512.2 reference genome, with mutations in samples 39 and 40 (39 mutations), samples 10, 33, 34 and 43 (40 mutations), samples 36, 37, 41, 44 and 45 (41 mutations), sample 35 (43 mutations) and samples 42 and 46 (44 mutations). The mutations mentioned above were calculated together with shared and unique mutations. Samples 38 and 46 do not carry unique mutations (Figure 2).

The most common nucleotide substitution detected (see Table S1) was from cytosine to thymine (51/134 mutations), followed by guanine to thymine (28/134 mutations) and thymine to cytosine and guanine to adenine (12/134 mutations), and was adenine to guanine (10/134 mutations). All mutations were homozygous. We observed 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 gene encoding the S protein. One unique (c.100C>T and c.463C>T) and one shared (c.245T>C) missense mutation was identified in the M (matrix) 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 found in N (nucleocapsid) region. No mutations were found in E (envelope) region (Tables S1 and S2).

Furthermore, one frameshift, 41 missense mutations, 27 synonymous mutations, and two upstream gene mutations were found in the ORF1ab region. One synonymous mutation and ten missense mutations were detected in the ORF3a region. The ORF6 region showed one missense mutation, while the ORF7a region exhibited six missenses and one synonymous mutation, and the ORF8a region exhibited one missense, one upstream gene, and one conservative insertion-deletion mutation. Finally, one synonymous mutation (c.57C>T) was found in the region of ORF10. 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 (Table S1).

3.3. The phylogenetic tree was built from viral sequences

We wanted to examine the significant alterations in all sequences to understand what distinguished our cases from others worldwide. MEGA X was used to create a phylogenetic tree based on the 17 viral sequences using the greatest likelihood technique. (Figure 1). Variant caller (v.5.12.0.4) identified 134 mutations from 17 SARS-CoV-2 viral genome sequences in COVID-19 patient samples (Table S1).

Most of the mutations shared 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 (Table S1).

These variants have grouped with USA, Indian, and England variants. All are considered as 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).