SARS-CoV-2 lineage dynamics in England from September to November 2021: high diversity of Delta sub-lineages and increased transmissibility of AY.4.2

Lineage diversity

In round 14 the lineage was determined for 481 of 764 positive samples. All lineages were Delta or a Delta sub-lineage with the four most prevalent lineages being AY.4 at 65.1% (60.7%, 69.2%, n=313), AY.43 at 6.0% (4.2%, 8.5%, n=29), B.1.617.2 (original Delta variant) at 5.2% (3.6%, 7.6%, n=25) and AY.4.2 at 4.6% (3.0%, 6.8%, n=22) (Figure 1-A, Supplementary Table 1). In round 15 the lineage was determined for 840 of 1399 positive samples. Again all samples were Delta or a Delta sub-lineage with the most prevalent lineages again being AY.4 at 57.6% (54.2%, 60.9%, n=484), B.1.617.2 at 12.8% (10.8%, 15.3%, n=108), AY.4.2 at 11.8% (9.8%, 14.1%, n=99) and AY.43 at 4.8% (3.5%, 6.4%, n=40).

• Download figure
• Open in new tab

Figure 1:

Proportion of positives by Delta sub-lineage. (A) The proportion of positives samples by round designated to the 8 lineages most prevalent over both rounds 14 and 15 (AY.39, AY.4, AY.4.2, AY.43, AY.44, AY.5, AY.6, B.1.617.2). (B-D) Proportion of positive samples by round and region with lineage designated as (B) B.1.617.2, (C) AY.4 and (D) AY.4.2.

The next four most prevalent lineages over both rounds combined were AY.5, AY.6, AY.39, and AY.44. However, even a single detection of a lineage corresponded to an average of 971 (157, 5852) swab-positive individuals at any time during round 14 and 1051 (175, 6264) swab-positive individuals at any time during round 15. During rounds 14 and 15 there were 33 and 31 unique lineages detected, respectively with 44 unique lineages detected overall. There was no apparent difference in genetic diversity between the two rounds as estimated by the Shannon diversity (p = 0.831) (Supplementary Table 2).

Distribution by region and age

During round 15 the proportion of B.1.617.2 was found to be highest in London at 22.1% (14.9%, 31.4%), being greater than the proportion in South East, East of England and Yorkshire and The Humber (Figure 1-B, Supplementary Table 3). Conversely, in round 14 and 15 the proportion of AY.4 was lowest in London at 48.1% (35.4%, 61.1%) and 44.2% (34.6%, 54.2%) respectively and was found to be higher in North West, West Midlands and Yorkshire and The Humber during both rounds (Figure 1-C, Supplementary Table 3). This reduced proportion of the nationally most prevalent lineage (AY.4) in London coincided with a higher level of genetic diversity in London. The Shannon diversity was highest in London during both rounds at 1.814 in round 14 and 1.809 in round 15 (p <0.001 and p = 0.002 respectively, reference = West Midlands, Supplementary Table 2). Higher levels of genetic diversity were also found during both rounds in the South East and South West, relative to the West Midlands (which showed the lowest levels of genetic diversity in round 14 and the second lowest in round 15). There were no regional differences in the proportion of AY.4.2 during round 15 (Figure 1-D, Supplementary Table 3). Regional differences during round 14 and regional differences for other lineages could not be investigated due to small sample sizes but numbers are provided in Supplementary Table 4.

Sub-regional analysis was performed in order to investigate the presence of clustering in each round for each lineage (see Methods). Despite being highly geographically dispersed (Figure 2) clustering was detected in round 14 for AY.4 (p = 0.037) and AY.4.2 (p = 0.029) (Supplementary Table 5). However, during round 15 clustering was no longer evident for both AY.4 (p = 0.706) and AY.4.2 (p = 0.067). The only lineage for which clustering was detected in round 15 was AY.6 (p = 0.003) which was found mainly in London and towards the South coast of England.

• Download figure
• Open in new tab

Figure 2:

Geographic distribution of all positive samples with a lineage designation (Green) with overlaid distribution of AY.4.2 (Pink, left), AY.4 (Purple, centre) and AY.6 (Orange, right) for both round 14 (top) and round 15 (bottom). The lineages shown had either a significant level of clustering in round 14 (AY.4 and AY.4.2) or round 15 (AY.6).

During round 15 the proportion of B.1.617.2 was higher in individuals ages 25-34 years old at 24.2% (12.8%, 41.0%) relative to those aged 35-44 years old at 8.0% (4.1%, 15.0%) (p = 0.026) (Supplementary table 6). The proportion of AY.4 was found to be lower in 5-12 year olds at 52.1% (44.6%, 59.5%) relative to 35-44 year olds in which the proportion of AY.4 was 65.0% (55.3%, 73.6%) (p= 0.042) in round 15, while it was not in round 14.There were no differences between age groups in the proportion of AY.4.2 during round 15. Differences between age groups during round 14 for AY.4.2 and other lineages could not be investigated due to small sample sizes but numbers are provided in Supplementary Table 7.

Detection of increasing sub-lineages

Logistic regression models were fitted to the proportion of each lineage detected in either round 14 or 15, allowing daily growth rates in proportion to be estimated (Figure 3, Supplementary Table 8). Of the 44 unique lineages detected, 6 were estimated to have growth rates different to zero. AY.4, AY.39, AY.98.1 and AY.111 were decreasing in proportion, whereas AY.4.2 and B.1.617.2 were increasing in proportion. The decrease in proportion of AY.4 corresponded to a daily growth rate of -0.009 (−0.015, -0.003). The increase in proportions of B.1.617.2 and AY.4.2 corresponded to growth rates of 0.029 (0.017, 0.041) and 0.028 (0.016, 0.041) respectively.

• Download figure
• Open in new tab

Figure 3:

(A) Estimated daily growth rate of the log odds of each lineage detected relative to all other lineages. Shown are both lineages with a growth rate not significantly different to zero (black) and those with a growth rate significantly different to zero (coloured). (B-G) Raw estimates of the daily proportions (points) with 95% confidence intervals (error bars) for lineages with a growth rate significantly different to zero: B.1.617.2 (B, pink), AY.4.2 (C, yellow), AY.4 (D, dark green), AY.111 (E, orange), AY.39 (F, purple), AY.98.1 (G, light green). Also shown is the best-fit Bayesian logistic regression model with central estimate (solid line) and 95% credible interval (shaded region).

Comparing estimates of the reproduction number R from round 14 to round 15 for AY.4 and AY.4.2 (see Methods) we estimate a multiplicative R advantage of 1.15 (1.08, 1.23), assuming no change in the generation time distribution.

Differences in cycle threshold values

There were quantitative differences between lineages in the N- and E-gene Ct values. The mean N- and E-gene Ct values were lowest for AY.6 though not materially lower than the values obtained for AY.4 (Figure 4, Supplementary Table 9). Mean N-gene Ct value was 22.14 (20.30, 23.99) for AY.6 compared to 23.98 (23.68, 24.28) for AY.4 (p = 0.054). Mean E-gene Ct value was 20.74 (18.90, 22.59) for AY.6 compared to 22.46 (22.16, 22.76) for AY.4 (p= 0.071). Mean N- and E-gene Ct values were found to be comparable to AY.4 for both AY.4.2 and AY.5. Relative to AY.4, mean N- and E-gene Ct values for AY.43, AY.44, AY.39 and B.1.617.2 were all higher.

• Download figure
• Open in new tab

Figure 4:

Estimated mean N-gene (A) and E-gene (B) Ct values for the 8 lineages most prevalent over rounds 14 and 15 (AY.39, AY.4, AY.4.2, AY.43, AY.44, AY.5, AY.6 and B.1.617.2) as calculated using Gaussian regression. Point estimates (points) and 95% confidence intervals (lines) are shown for estimates obtained using data from both rounds (blue), data from just round 14 (green) and data from just round 15 (purple)

Differences in symptomatology

The proportion of individuals exhibiting the most predictive COVID-19 symptoms (loss or change of sense of taste, loss or change of sense of smell, new persistent cough, fever) in the month prior to swabbing was lower (p=0.026) in those infected with AY.4.2 at 38.1% (29.9%, 47.1%) relative to those infected with AY.4 at 49.2% (45.7%, 52.8%) (Figure 5-A, Supplementary Table 10). This difference was not explained by patterns in age or N-gene Ct value (Figure 5-B, Supplementary Table 11). In all logistic models for symptom status adjusting for round of the study, the point estimate of the effect of AY.4.2 infection was broadly similar, but, due to increased variability in the estimate, did not reach statistical significance. In particular, the odds ratio of AY.4.2-infected participants exhibiting the most predictive COVID-19 symptoms was 0.69 (0.46, 1.03) lower than for AY.4 infected participants (p= 0.069) when including round as a covariate (Model 3, Supplementary Table 11).

• Download figure
• Open in new tab

Figure 5:

(A) Proportion of positive individuals reporting any symptoms or reporting one of the four most predictive COVID-19 symptoms (loss or change of sense of taste, loss or change of sense of smell, new persistent cough, fever) in the last month by lineage of infection, for the 8 lineages most prevalent during rounds 14 and 15 (AY.39, AY.4, AY.4.2, AY.43, AY.44, AY.5, AY.6 and B.1.617.2). (B) Odds ratios of reporting the most predictive COVID-19 symptoms in the last months for multivariable logistic regression models including lineage (AY.4.2 with reference AY.4, red), age (relative to change of 10 years in age, blue), round of study (round 15 with reference round 14, green) and N-gene Ct value (relative to change in Ct value of 5, purple).

In addition, 68.6% (59.8%, 76.3%) of those infected with AY.4.2 reported any symptoms in the month prior to swabbing compared to 75.4% (72.2%, 78.3%) for those infected with AY.4 (p= 0.119). There were no differences evident in symptom reporting between AY.4 infected individuals and the other 6 most prevalent lineages (B.1.617.2, AY.5, AY.6, AY.43, AY.44 and AY.39).

Phylogeographic analysis

A relaxed molecular clock model was fit to the data and used to estimate a time-resolved phylogenetic tree (Figure 6). AY.4.2 was found to populate two closely related clades that emerged in June/July 2021. AY.43, AY.5 and AY.6 were also observed to have distinct clade groupings having emerged around June/July 2021 as well. The mutation rates inferred at the tree’s tips showed a large degree of variation in all of the 8 most prevalent lineages. The mean mutation rate for AY.4.2 was found to be 0.57 (<0.01, 1.10) × 10⁻⁴ lower than the mean mutation rate of AY.4 (p = 0.050) (Figure 6, Supplementary Table 12). The mean mutation rate inferred for samples collected in round 15 was found to be 1.00 (0.70, 1.40) × 10⁻⁴ lower than the mean mutation rate for samples collected in round 14 (p < 0.001).

• Download figure
• Open in new tab

Figure 6:

(A) Time-resolved phylogenetic tree of all positive samples obtained for which the lineage designated was Delta or a Delta sub-lineage. (B) Distribution of mutation rates inferred at each phylogenetic tree tip for the 8 lineages most prevalent in round 14 and round 15 for samples obtained in round 14 and 15. (C) Distribution of mutation rates inferred at each phylogenetic tree tip for samples collected in round 14 and round 15. (D) Difference in the mean mutation rates between lineages and rounds as inferred from a Gaussian regression model.

A mugenic model was run on the time-resolved phylogenetic tree to estimate the relative virus migration rates between regions, a measure of inter-region transmission (Supplementary Table 13). Overall levels of inter-region transmission were lowest for the North East during round 14 and 15. The highest overall level of inter-region transmission was observed for the North West during round 14 and 15, but looking at individual rounds there were higher levels for Yorkshire and The Humber in round 14 and for the South East in round 15. High rates of transmission during round 14 and 15 were found between the North West and Yorkshire and The Humber, the West Midlands and the South East, and also between the South East and London.

Other lineages

Though we have focused on AY.4.2 we have detected a diverse set of Delta sub-lineages, with even a single detection corresponding to approximately 1000 swab-positive infections in the community at one time during the study period. The short time over which AY.4.2 went from being an undeclared lineage to a variant under investigation shows how crucial it is to have careful surveillance of all lineages irrespective of frequency. For 38 of the 44 detected lineages, it was unable to be determined whether the proportion was increasing or decreasing.

Between rounds 14 and 15 a reduction in the mean mutation rate of the virus was detected suggesting a reduction in the rate of evolution. However, despite this slowdown evolution is still occurring and we observed an increase in the proportion of B.1.617.2, an indicator that the number of undeclared B.1.617.2 sub-lineages was increasing, suggesting even further diversity of Delta sub-lineages that have yet to be given a unique lineage designation.

Further, though we capture the dynamics within England, SARS-CoV-2 is a global problem and new variants of concern can arise anywhere in the world and then spread through international travel. Higher proportions of B.1.617.2 were detected in London as well as higher levels of diversity; this likely reflects the role London continues to play in the introduction of international variants [33]. Within England, the North West region played a major role in the dissemination of the virus, having the greatest inferred rate of inter-region transmission.

Analysis of N- and E-gene Ct values found decreased levels in AY.4 and AY.4.2, which is unsurprising given both have successfully disseminated across the country, but AY.5 and AY.6 were also found to have similarly low Ct values suggesting similar viral loads; the mean N- and E-gene Ct value appeared slightly lower for AY.6 compared to AY.4. Clustering was also detected in round 15 for AY.6; careful consideration of AY.6 should be given in the future in case the current lack of growth so far reported [11] has only been due to its geographic isolation.

Limitations

We have presented the inferred dynamics between Delta sub-lineages in England between 9 September and 5 November 2021. Our sample’s main strength over those obtained from routine surveillance is the random nature of the testing program leading to a relatively unbiased set of positive samples. However, as the sample sizes we obtain are relatively small compared with routine national surveillance our estimates have lower precision.

Lineages were only successfully determined for ∼61% of positive samples, with the ability to determine a lineage heavily influenced by a sample’s Ct value; this has potentially led to biases with lineages with lower Ct values more heavily represented in the dataset. Detecting distinct sub-lineages is a high-dimensional problem, with often many common mutations being shared between distinct lineages with only a small number of distinguishing mutations. This is exacerbated when all the lineages are highly related, as in the current nature of the pandemic in England where all samples are descendants of Delta (B.1.617.2), and can lead to incorrect designations [23]. Further, only sub-lineages that have been defined are able to be assigned to a sample. During the emergence of a new sub-lineage there is a phase of ambiguity when numbers are small and it is unclear if the mutations present warrant the declaration of a new sub-lineage. This can be seen in the detection of AY.4.2 and AY.43; both lineages had been circulating for months by October 2021 [11] but were not yet declared sub-lineages by pangoLEARN [24] in early October 2021, and so did not appear in the publicly available technical briefings [25].

Viral genome sequencing

The methods of the REACT-1 study have been described elsewhere [37]. REACT-1 is a repeat cross-sectional study whereby in each round a random subset of the English population (selected from the National Health Service general practitioners’ patient list) is invited to obtain a self-administered swab test (parent/guardian administered for 5-12 year olds). These tests are then sent to a laboratory to undergo rt-PCR testing for the presence of SARS-CoV-2. A round of the study covers a ∼2- to 3-week period and has occurred approximately monthly since May 2020 with between 100,000 and 185,000 individuals taking part in each round. Since round 8 in January 2021 all positive samples with a low enough N-gene Ct value (the threshold was 34 in rounds 14 and 15 presented here) and sufficient volume have been sent for genome sequencing. Amplification of the extracted RNA was performed using the ARTIC protocol [38] (version 4 primers), with sequence libraries prepared using CoronaHiT [39]; sequencing was performed on the Illumina NextSeq 500 platform. Raw sequences were analysed using the bioinformatic pipeline [40] and then uploaded to CLIMB [41]. Lineages were assigned using PangoLEARN [24] (database version 2021-11-04), a machine learning-based assignment algorithm, using Pango nomenclature [9]. For some sequences of low overall quality, a lineage designation was not possible and so they were not included in the analyses. Samples with less than 50% of bases covered were further excluded from the analysis.

Phylogeographic model

For all sequences from REACT-1 rounds 11 (15 April - 3 May 2021), 12 (20 May - 7 June 2021), 13 (24 June - 12 July 2021), 14 (9 September - 27 September 2021) and 15 (19 October - 5 November 2021), in which the lineage designated was Delta or a Delta sub-lineage, a maximum likelihood phylogenetic tree was constructed using a HKY model implemented in IQ-TREE [42]. An uncorrelated relaxed clock model implemented in TreeTime [43], assuming a normal distribution of rates with mean 0.0008 substitutions per site per year and a single coalescent rate for the time scale, was then fit to the maximum likelihood phylogenetic tree producing a time-resolved phylogenetic tree. The mutation rates at the tree’s tips were extracted from the model and a Gaussian regression model was fit to the samples obtained during round 14 and 15 for the 8 most prevalent lineages (AY.39, AY.4, AY.4.2, AY.43, AY.44, AY.5, AY.6, B.617.2) including lineage and round as covariates. A mugenic model (implemented in TreeTime [43]) was run on the time-resolved phylogenetic tree, treating the region in which each sample was isolated as a discrete state. This allowed estimates of the migration rates between regions to be calculated (assumed to be symmetric).

Statistical analyses

The 95% confidence intervals for lineage proportions were calculated using the Wilson method [44] assuming a Binomial distribution. This method is preferred when the number of positives is low but is still valid when this is not the case [45]. Higher accuracy in confidence interval estimates for when the number of positives is low was chosen so that lower bounds on case numbers for rarer lineages were as accurate as possible.

Estimates of the true number of swab-positive infections in England during round 14 and round 15 for lineages in which only one sample was detected in a round were calculated by multiplying the estimated proportion of the lineage for each round, the weighted prevalence estimated for each round [46], and the population size of England [47].

The significance of differences in proportions of particular lineages by age group and region was calculated using Fisher’s exact test with a binary outcome variable (lineage of interest or not). Differences with a p-value less than 0.05 were considered statistically significant. Analysis was only completed for a lineage in a round if there were more than 90 samples (AY.4 round 14, AY.4 round 15, AY.4.2 round 15, B.1.617.2 round 15), so that there were, on average, more than 10 samples per parameter (9 regions in England).

Shannon diversity was calculated using all data for round 14 and round 15, and for each region for round 14 and round 15 [48]. The significance of any differences in Shannon diversity between round 14 and 15 (for all data) and between regions in each round was assessed using the Hutcheson T-test [49] and its associated p-value.

The relative growth rate of a lineage compared to all other lineages was estimated using a Bayesian logistic regression model fit to the binary outcome variable (lineage of interest or not) over time. The probability that the growth rate was greater than zero was calculated from the model’s posterior. Lineages were deemed to be different to zero if the posterior probability that the growth rate was greater than zero was greater than 0.975 or less than 0.025, similar to a p-value threshold of 0.05.

The growth rates of AY.4.2 and AY.4 infected individuals were estimated by fitting an exponential model to the daily weighted prevalence using all REACT-1 data (all negatives and all AY.4/AY.4.2 associated positives) for rounds 14 and 15 assuming a Binomial likelihood. Growth rates were then converted to estimates of the reproduction number R assuming a gamma-distributed generation time with the shape parameter, n =2.29, and rate parameter, b =0.36 [50] through the equation [51]. The multiplicative R advantage of AY.4.2 over AY.4 was estimated using the entire posterior distribution of R_AY.4.2/R_AY.4 with the median and 95% credible interval reported.

For each lineage with more than 1 sample in a round the presence of clustering was assessed. The pairwise distance matrix between all n samples that were designated to a specific lineage was calculated and from this a mean pairwise distance was calculated for the lineage. Next, 10,000 random combinations of n positive individuals (n positive individuals chosen each time without replacement), for which any lineage was determined, were selected and for each combination the distance matrix and mean distance was calculated. The proportion of the 10,000 estimated mean distances below the lineage-specific mean distance was then calculated. Clustering was deemed to be significant if this proportion was less than 0.05.

For the 8 most prevalent lineages across rounds 14 and 15 Gaussian regression was performed to estimate the mean N- and E-gene Ct values for each lineage and p-values used to assess the significance of any difference to the reference lineage (AY.4). Models were run on all data (rounds 14 and 15 combined) and then run on data from each individual round as a sensitivity analysis.

The proportion of individuals reporting any symptoms in the month prior to swabbing and any of the most predictive COVID-19 symptoms in the month prior to swabbing was calculated for the 8 most prevalent lineages across rounds 14 and 15. P-values were estimated for each lineage relative to AY.4 by performing logistic regression with the symptom status as a binary variable (any symptoms vs no symptoms, and separately most predictive COVID-19 symptoms vs none of the most predictive COVID-19 symptoms). The sensitivity of the results that AY.4.2 is less likely to exhibit the most predictive COVID-19 symptoms, relative to AY.4, was assessed by fitting further logistic regression models including age, round of study and N-gene Ct value as covariates (E-gene was also investigated but was no different to using N-gene and so this was not included).