Genomic Epidemiology of SARS-CoV-2 variants in unvaccinated children found to be similar to that of adults with COVID-19 vaccinations

Study design and ethical considerations and approval

This was a retrospective genomic epidemiological study conducted using SARS-CoV-2 genome sequences from Pakistan, spanning May 2021 to October 2022.

This study was conducted in accordance with the ethical standards set forth in the Declaration of Helsinki. Ethical approval was obtained from the Ethics Review Committee of Aga Khan University, which granted a waiver of written and verbal informed consent for this retrospective analysis (Approval No. 2021-6316-18156 and 2022-6871-22433). The waiver was justified because this was a retrospective analysis using secondary data, where re-contacting participants to obtain consent was neither feasible nor practical. Such an approach is consistent with ethical guidance for studies involving anonymized or previously collected data, where the risks to participants are minimal. All data was de-identified prior to analysis to ensure confidentiality.

Data sources and sample selection

We retrieved all available SARS-CoV-2 genome sequences (FASTA format) and associated metadata for Pakistan from the GISAID database (https://www.gisaid.org). Filters were applied to include only sequences with “high coverage” and “complete collection date.” From this dataset, 569 sequences were selected for downstream analysis based on the availability of vaccination status in the metadata.

Additionally, data for SARS-CoV-2–positive samples from individuals <18 years of age tested at Aga Khan University Hospital (AKUH), Karachi, were included for pediatric-focused analysis (n = 143). The specimens were selected as previously described by Nasir et al. (Nasir et al., 2023), by identifying the first 10 positive samples per day. Inclusion criteria included a PCR cycle threshold (Ct) value ≤35 using the Cobas SARS-CoV-2 assay. Samples with Ct >35 or duplicates from the same individual were excluded.

Classification by vaccination status

Vaccination status was assigned based on metadata and age. Individuals ≥18 years were considered vaccinated or unvaccinated based on reported status. Since COVID-19 vaccines were not available to individuals <18 years in Pakistan until late 2021 (GoP, 2021), all pediatric cases (<18 years) during the study period were considered unvaccinated unless otherwise specified.

Phylogenetic and genomic analysis

Phylogenetic analysis was conducted using the Nextstrain Augur pipeline and visualized via the Auspice platform (https://auspice.us). Sequences were aligned to the SARS-CoV-2 reference genome (NC_045512.2) using MAFFT. A maximum-likelihood phylogenetic tree was constructed using IQ-TREE2, employing the General Time Reversible (GTR) substitution model. Ancestral sequence inference and molecular dating were performed within Augur, and Nextstrain clade assignments were included. Divergence values were extracted from the JSON output of Augur for downstream statistical comparisons. Visualization filters (e.g., by age group, vaccination status, and time) were applied in Auspice for subgroup analysis.

Statistical analysis

Statistical analyses were performed using IBM SPSS v24 and GraphPad Prism v8.0.1. The Kolmogorov–Smirnov test was used to assess normality of divergence data. Chi-square, unpaired samples t-test and Mann Whitney U tests were used to analyze for significant differences. For comparisons involving more than two groups (e.g., Ct values across variants), the Kruskal–Wallis test followed by post-hoc pairwise comparisons was used. A p-value < 0.05 was considered statistically significant.

Study cohort description

We wanted to understand the effect of COVID-19 vaccination on the transmission of SARS-CoV-2 variants in Pakistan by examining the genomic epidemiology of strains circulating in 2021 and 2022. The 569 SARS-CoV-2 genomes we analyzed were those available in GISAID for samples collected between May 2021 and October 2022 and which had metadata available including COVID-19 vaccination status (Supplementary table 1). Overall, Omicron and Delta variants were dominant at 45.87% and 45.17% respectively ( Table 1). 54.1% of the strains had been collected in 2021 of which the dominant variant was Delta (82.8%). The remaining (45.9%) strains were from 2022 when Omicron was dominant (92%). 47.1% genomes were collected from females. The greatest number of sequences (71%) were from strains collected from children (<18 years, p < 0.001).

Of the COVID-19 cases studied for vaccination status, there were 29.3% (n = 167) from vaccinated individuals and 70.7% (n = 402) from unvaccinated individuals ( Table 2). Most (69.5%) of the vaccinated persons were in the 18-50 years old category (p < 0.001). Most of genomes studied from unvaccinated COVID-19 cases were from children aged <18 years.

The SARS-CoV-2 genome submission were from across Pakistan. Of the 569 samples, most (38.7%) were from Islamabad (Islamabad Capital Territory, ICT), followed by Sindh (32.7%), Punjab (10.5%), Kyber Pakhtunkhwa (8.8%), Baluchistan (2.3%) and the rest from Gilgit-Baltistan and Azad Jammu and Kashmir ( Figure 1).

Figure 1. Province-wise distribution of samples.

Graph shows geographical distribution of SARS-CoV-2 genomes downloaded from GISAID for phylogenetic analysis. Samples were from the provinces of Pakistan; ICT – Islamabad Capital Territory, KPK - Khyber Pakhtunkhwa. *includes samples from Gilgit-Baltistan and Azad Jammu Kashmir provinces as well as those for which location was unavailable.

SARS-CoV-2 phylogenetic diversity

The phylogenetic diversity of SARS-CoV-2 genomes sequences during the period of study showed the divergence of genomes between the Delta variants in 2021 and Omicron variants in 2022 was evident ( Figure 2A). The most common Delta lineage was B.1.167.2 (53.7%; n = 137) followed by AY.108 (24%; n = 60 (Supplementary Table 2)). In 2021 there were 1.3% Wildtype, 3.9% Alpha and 3.2 % Beta strains. In 2022, within the Omicron lineage BA.5.2 was prevalent (34.5%, n = 87). Importantly, isolates from both unvaccinated and vaccinated groups were found amongst the SARS-CoV-2 lineages identified. We further compared the occurrence of variants between vaccinated and unvaccinated groups. Strain divergence is visualized through a scatterplot and a similar pattern was observed in both vaccinated and unvaccinated cohorts ( Figure 2B). Genetic differences between vaccinated and unvaccinated individuals did not show any statistically significant difference (p = 0.75, unpaired samples t-test), (Supplementary Table 3).

Figure 2. SARS-CoV-2 phylogenetic divergence based on vaccination and time of sample collection from May 2021 to October 2022.

569 SARS-CoV-2 genomes were analyzed using the Augur pipeline. (A) The graph shows the diversity of SARS-CoV-2 genomes collected from unvaccinated (dark blue dots) or unvaccinated persons (light green dots) respectively. The X-axis shows samples along the timeline of the pandemic from 2019 to 2022, while the Y axis shows the two main branches of the phylogenetic tree. (B) The scatterplot shows genetic divergence between SARS-CoV-2 strains from either unvaccinated (dark blue) or vaccinated (light green) individuals. The X-axis of the graph shows divergence values of each sample in the cohort, and Y-axis shows the groups. The data was downloaded from GISAID (gisaid.org).

Age-wise distribution of SARS-CoV-2 variants

Next we compared the occurrence of variants across different age groups. Delta and Omicron were the most abundant variants across all age groups ( Figure 3A) and age-wise comparison of the strains showed a distribution of the prevalent SARS-CoV-2 variants (Beta, Alpha, Delta and Omicron) across older and younger age groups, as divided between those below 18 years and those aged 18 years and above ( Figure 3B, supplementary table 3). We divided the dataset for cases less than 18 years and more than 18 years. The statistical test between the two age categories did not give a significant difference in terms of divergence. Genomic divergence analysis of variants with age also displayed the spread of those <18 years and >18 years old cases to be similar ( Figure 3B). Also, unpaired samples T-test on the genomic divergence values did not give a significant statistical difference between the two groups (p-value = 0.6038, supplementary table 3).

Figure 3. Age-wise comparison of SARS-CoV-2 variants.

We compared the occurrence of SARS-CoV-2 variants based on the age of cases. (A) The graph depicts Next clade-based phylogeny of strains from individuals according to their age. The scatter plot shows the clade on the X-axis and age (years) on the Y-axis as represented by the dots for each individual.

(B) Strains were plotted to show genetic divergence as per age of individuals. SARS-CoV-2 divergence between samples is shown on the X axis whilst the Y-axis shows age. The dotted horizontal line marks at age 18 years to give a visual comparison of individuals aged less than 18 years and those above 18 years.

Comparison of SARS-CoV-2 variants in children from Karachi with nationwide trends

There has been limited information on COVID-19 in children aged <18 years and we now focused on this age group by comparing trends based on genomic sequences available countrywide and compared with those from our facility-based testing in Karachi. Of the SARS-CoV-2 genomes analyzed from across Pakistan, 402 were from those aged <18 years with 218 isolated in 2021 and 184 from 2022. In 2021, 89.9% of strains from this younger age group comprised Delta variants, and in 2022, Omicron was the most prevalent (90.2%) ( Figure 4A). We then focused on SARS-CoV-2 isolates received at the Aga Khan University Hospital as more information including, viral loads (Ct values of diagnostic PCR testing) were available. Of the 143 isolates collected from the under 18 year-olds with COVID-19 at our institute, 103 were from 2021 and 40 from 2022. Of these, 70% were from outpatient cases and 30% were from inpatients COVID-19 patients. We observed the dominant variants to be Delta in 2021 (57.3%) and Omicron in 2022 (100%) ( Figure 4B). Unfortunately, data on clinical disease severity of all individuals was not available.

Figure 4. Study of SARS-CoV-2 variants in children in the study cohort.

We first identified Wildtype and variants Alpha, Beta, Delta and Omicron amongst the strains studied. (A) The graph shows the frequency of variants in 402 genomes sequenced from children (<18 years) within the study cohort. The graph depicts the frequency of isolates from 2021 (n = 218) and 2022 (n = 184). (B) Identification of VOC from children diagnosed at AKUH. The graph depicts the frequency of 143 variants identified (2021 n = 103; 2022 n = 40). (C) The SARS-CoV-2 diagnostic PCR cycle threshold value was compared between the variants identified in those aged under 18 years at AKUH. Kruskal-Wallis test shows the medians vary significantly among the 5 groups, P = 0.0002. Mann-Whitney U test showed significant differences between Wildtype and Delta p=0.0007; Wildtype and Omicron p = 0.0001; Alpha and Omicron p-0.0449; Beta and Omicron p-0.0249; *CT1 denotes cycle threshold of orf1ab gene target. *p <= 0.05 was considered significant.

Further, we examined the SARS-CoV-2 viral loads of 143 AKUH specimens from those aged <18 years and found SARS-CoV-2 viral loads (as assessed by Ct values), to differ significantly between variants (Kruskal Wallis p value = 0.0002). There was a higher Ct in wild type compared with VOC Delta (p = 0.0007) and Omicron (p = 0.0001) ( Figure 4C). This indicated higher viral loads in cases involving VOCs as compared with Wildtype strains. There was no significant difference between Wildtype and Alpha, p = 0.6653; Wildtype and Beta, p = 0.3239; Alpha and Beta, p = 0.6238; Alpha and Delta p = 0.0957; Beta and Delta, p = 0.1350; Delta and Omicron, p = 0.2408. Alpha variant viral loads (median CT = 27) were lower than in the case of Omicron (median CT = 21) (p = 0.0449) infections; and the same trend was observed for Beta variants (p = 0.0249). Of note, samples containing Omicron variants had the lowest median Ct value (higher viral load) in children.

Limitations

This study has several limitations. Vaccination status was not available for many adult cases in GISAID, leading to underrepresentation of this group. Additionally, clinical data such as symptom severity, hospitalization, and comorbidities were unavailable, limiting our ability to correlate viral variants with disease outcomes. Finally, the retrospective design limits causal inferences.

Implications

Our findings suggest that in settings with widespread prior exposure, variant distribution and transmission dynamics may not differ significantly by vaccination status - a key consideration for future surveillance and vaccine deployment strategies. The high seroprevalence, young population structure, and mixed vaccine types in Pakistan offer a unique lens through which to understand SARS-CoV-2 evolution in diverse immunological landscapes.

Continued genomic surveillance, particularly in pediatric and rural populations, remains vital to detect emerging variants and guide context-specific public health responses. Moreover, insights from populations with high natural immunity and mixed vaccine platforms can inform global preparedness for future respiratory virus outbreaks.