Genomic epidemiology of seasonal influenza circulation in China during prolonged border closure from 2020 to 2021.

China experienced a resurgence of seasonal influenza activity throughout 2021 despite intermittent control measures and prolonged international border closure. We show genomic evidence for multiple A(H3N2), A(H1N1), and B/Victoria transmission lineages circulating over 3 years, with the 2021 resurgence mainly driven by two B/Victoria clades. Phylodynamic analysis revealed unsampled ancestry prior to widespread outbreaks in December 2020, showing that influenza lineages can circulate cryptically under non-pharmaceutical interventions enacted against COVID-19. Novel haemagglutinin gene mutations and altered age profiles of infected individuals were observed, and Jiangxi province was identified as a major source for nationwide outbreaks. Following major holiday periods, fluctuations in the effective reproduction number were observed, underscoring the importance of influenza vaccination prior to holiday periods or travel. Extensive heterogeneity in seasonal influenza circulation patterns in China determined by historical strain circulation indicates that a better understanding of demographic patterns is needed for improving effective controls.

Seasonal influenza epidemics cause a high level of disease and economic burden each year (Li et al. 2018; Lafond et al. 2021). In non-pandemic years, they exhibit strong seasonal cycles in temperate regions and semi-annual epidemics or year-round activity in tropical and subtropical regions (Azziz Baumgartner et al. 2012; Tamerius et al. 2013; Li et al. 2019). In mainland China, seasonal influenza patterns are complex and vary by latitude, geographic location, and virus type (Yang et al. 2018; Diamond et al. 2022). Similar to other temperate regions, winter epidemics are short in northern China, while southern China experiences semi-annual or year-round circulation with varied intensity (Shu et al. 2010). Over the last decade in China, A(H1N1) and two type B influenza viruses, B/Victoria and B/Yamagata, were sparsely detected during summer epidemics, in contrast to A(H3N2) viruses that circulated nationwide during winter and in lower latitudes in summer (Yang et al. 2018; Diamond et al. 2022). Notably, following a relatively mild 2017–2018 influenza season in China, severe outbreaks occurred during the 2018–2019 season alongside other countries in the Northern Hemisphere season, and while the 2019 Southern Hemisphere was severe (Barr et al. 2019; Mook et al. 2020), the normal 2019 summer epidemic was absent in China (Huang et al. 2020).

The epidemiology of seasonal influenza is characterized by source–sink dynamics with rapid global migration of all four lineages (Rambaut et al. 2008; Russell et al. 2008; Bahl et al. 2011; Bedford et al. 2015). Densely populated regions in Asia and Africa can sustain virus lineages for longer periods and act as source populations for winter epidemics in temperate regions by generating novel antigenic variants (Russell et al. 2008; Bedford et al. 2015; Dhanasekaran et al. 2022). Regional epidemics are caused by repeat introductions and co-circulation of multiple lineages with heterogeneity in prevalence. A(H3N2) viruses often dominate most influenza seasons and infect all age groups owing to faster antigenic evolution, whereas A(H1N1), B/Victoria, and B/Yamagata lineages exhibit slower antigenic evolution although they can prevail in some years (Bedford et al. 2015).

Despite a global decline in seasonal influenza activity during the Coronavirus disease 2019 (COVID-19) pandemic, influenza cases were reported to the World Health Organization (WHO) Global Influenza Surveillance and Response System (GISRS) throughout 2020. Genetic data submitted to Global Initiative on Sharing All Influenza Data (GISAID) indicates that A(H3N2), A(H1N1), and B/Victoria viruses circulated predominantly in South Asia, West Africa, and China, respectively, while B/Yamagata lineage viruses appear to have been eliminated (Woolthuis, Wallinga, and van Boven 2017; Koutsakos et al. 2021; Dhanasekaran et al. 2022; WHO 2022). As travel restrictions gradually eased in 2021, A(H3N2) viruses became increasingly prevalent across much of the Northern Hemisphere (Bolton et al. 2021), but in China, where international borders remained closed, the incidence of B/Victoria increased since late 2020 and throughout 2021 (CNIC 2021; Dhanasekaran et al. 2022).

China maintained an elimination strategy throughout 2020 and 2021 to kerb the spread of COVID-19. After the initial lockdown was lifted in April 2020, Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) outbreaks were locally suppressed through rapid isolation, large-scale testing, and other non-pharmaceutical interventions (NPIs) (Kraemer et al. 2020), enabling the majority of mainland China to maintain a baseline of light social measures and public health guidelines such as testing, mask mandates, and periodic interprovincial border closures while international borders remained closed (Zhang et al. 2021). Here, we investigate the evolution and epidemiology of B/Victoria in China amid heightened testing, international border closures, and NPIs associated with the COVID-19 pandemic during 2020–2021.

Methods

Dataset curation and phylogenetic analysis

We retrieved all available data on confirmed influenza B/Victoria cases in mainland China from January 2011 to November 2021 from the WHO GISRS (Hay and McCauley 2018). We also downloaded all available human seasonal influenza virus genomes from China (n = 2,299) and a representative subset of global sequences (n = 955), including all B/Victoria virus vaccine reference strains, from GISAID (accessed 17 February 2022) for the same period. After removing duplicate isolates reporting the same strain name, we aligned each gene segment using MAFFT independently (v.7.22) (Katoh and Standley 2013) and trimmed 5ʹ and 3ʹ coding regions in AliView (v.1.27) (Larsson 2014). We generated maximum likelihood (ML) phylogenetic trees using IQ-TREE (v.2) (Nguyen et al. 2015) with the best fit nucleotide substitution model (transvertion model (TVM) + F + R3) identified by ModelFinder (Kalyaanamoorthy et al. 2017). Branch support was assessed by an ultrafast bootstrap algorithm with 1,000 replicates (Hoang et al. 2018). Clades were labelled according to WHO nomenclature. We tested the temporal signature necessary for time-scaled analysis by regression of root-to-tip genetic distances from the ML tree against sample collection dates using TempEst (v.1.5.3) (Rambaut et al. 2016). Outliers that substantially diverged from the mean were assumed to contain excessive sequencing artefacts and were removed from further analysis. ML dating of the final haemagglutinin (HA) gene dataset was performed using TreeTime (Sagulenko, Puller, and Neher 2018) including preliminary molecular clock analyses for each segment.

Reassortment testing

To test for reassortment among major B/Victoria clades, we constructed ML trees for each gene segment using IQ-TREE (v.2) (Nguyen et al. 2015). All available B/Victoria sequences from mainland China and representative global sequences were included. By colouring each major clade, the phylogenetic topologies of each segment were compared to find potential reassortment events. To visualize reassortment, a maximum of 250 sequences were randomly selected from each major clade to plot a tanglegram using the ‘baltic3’ Python 3 library (https://github.com/evogytis/baltic).

Phylodynamic analysis

Whole genome sequences were concatenated for each lineage V1A.3-China (n = 140), 3a1 (n = 430), and 3a2 (n = 243). The final concatenated dataset was again tested for reassortment using RDP5 (Martin and Rybicki 2000) to apply seven methods: RDP, GENECONV, BootScan, MaxChi, Chimaera, SiScan, and 3Seq. Time-scaled phylogenetic analysis was conducted using BEAST (v.10.1.4) (Drummond and Rambaut 2007) with a general time reversible (GTR) + G4 nucleotide substitution model, an uncorrelated, relaxed molecular clock with a lognormal rate distribution prior and a non-parametric Gaussian Markov random field (GMRF) Bayesian Skyride coalescent tree model. Separately, we implemented a Bayesian Skyline coalescent tree model to estimate changes in the effective population size (N) for clades 3a1, 3a2, and V1A.3-China (Drummond et al. 2005). We also conducted an asymmetric discrete-trait phylogeographic analysis with Bayesian Stochastic Search Variable Selection in BEAST (v.10.1.4) (Drummond and Rambaut 2007) to reconstruct spatial diffusion and classify evidence as either ‘definitive’ (Bayes factor (BF) > 100) or ‘sufficient’ (100 > BF > 3). Chinese provinces were grouped into six geographic regions summarized in Supplementary Data S3. Finally, a birth–death serial skyline (BDSKY) model (Stadler et al. 2013) was implemented in BEAST (v.2.6.3) (Bouckaert et al. 2019) to estimate changes in the effective reproductive rate (R) over time for each lineage. R was estimated using a fixed tree topology estimated from the GMRF Bayesian Skyride coalescent model, where only the internal node heights were estimated with a strict clock model set to 1.35 × 10–3 nucleotide substitutions per site per year (nt/subs/site/year) and 1.65 × 10–3 nt/subs/site/year for 3a1 and 3a2 clades, respectively. The prior for the origin of BDSKY was used with the lower bounds set to 1 January 2019 (clade 3a1) and 1 January 2020 (clade 3a2) and the initial values to 2.2 (clade 3a1) and 0.96 (clade 3a2). To minimize the impact of regional sampling bias on phylogeographic analysis, we repeated the analysis for 3a1 and 3a2 using a subset containing <10 randomly sampled sequences per province per month (Supplementary Data S3 and Supplementary Fig. S8). For each model, Markov chain Monte Carlo runs were conducted at least twice each for 50–100 million steps, sampling every 1,000 to 20,000 steps, with 10 per cent discarded as burn-in. Sufficient sampling (effective sample size(ESS) > 200) from the posterior was inspected by Tracer (v.1.7.1) (Rambaut et al. 2018). Tree visualizations and circular plots were created with R packages ‘ggtree’ (Yu 2020) and ‘circlize’ (Gu et al. 2014). Map data in this study were downloaded from a public source (https://www.highcharts.com/).

Results

Persistence and resurgence of seasonal influenza activity in China

According to the Chinese National Influenza Centre, all seasonal influenza virus subtypes, including B/Yamagata, were detected in the first week of January 2020 when seasonal influenza activity peaked (CNIC 2020). The rapid decline thereafter coincided with lockdowns beginning in Wuhan and extending across Hubei and all other provinces (Lau et al. 2020). On average, NPIs associated with COVID-19 control reduced influenza activity by 79 per cent across China in comparison to past seasons (2011–2019) (Feng et al. 2021). Although little to no influenza activity was reported from April to November 2020 (Huang et al. 2020), cases in China began to rebound in December 2020 and fluctuated throughout 2021 (Supplementary Fig. S1). Remarkably, 98.7 per cent of confirmed influenza cases (n = 11,082) in China from 5 October 2020 to 5 September 2021 were B/Victoria viruses, 0.4 per cent were A(H1N1), 0.2 per cent were A(H3N2), and 0.2 per cent were B/Yamagata viruses (Huang et al. 2021). However, the reported B/Yamagata cases have not been confirmed by sequencing or virus isolation, consistent with the February 2022 WHO report (WHO 2022) where no confirmed B/Yamagata sequences were detected globally since March 2020. Based on available sequences in GISAID (accessed 7 November 2021), over 1,000 B/Victoria genome sequences were reported nationwide from 2020 to 2021, while one A(H1N1) sequence was reported in Fujian in December 2020, nine A(H3N2) sequences were reported in Yunnan province from May to August 2020, and zero B/Yamagata sequences were reported after January 2020.

V1A was the dominant B/Victoria lineage globally during 2016–2017, from which V1A.3 and subclades 3a1 and 3a2 subsequently emerged. From early 2019 to mid-2020, B/Victoria HA sequences in China belonged predominantly to clade V1A.3, while 3a1 and 3a2 were detected sporadically during September and November 2019, respectively (Fig. 1). Genomes sequenced since late 2020 mainly belonged to clades 3a1 and 3a2 with ∼99.5 per cent nucleotide identity (Supplementary Fig. S2), with only a few V1A and V1A.3 viruses detected. Clade 3a1 was almost exclusively detected in China and was derived directly from 3a collected prior to COVID-19, while clade 3a2 formed two monophyletic lineages, one of which circulated in China and the other globally (Fig. 1B, D). Similarly, the few A(H3N2), A(H1N1), and B/Victoria V1A and V1A.3 viruses detected since mid-2020 were most closely related to viruses detected during the 2019–2020 season in China (Supplementary Fig. S3). These results suggest persistent low-level domestic circulation of multiple influenza lineages within China despite highly stringent control measures in early 2020.

Figure 1.

Genomic surveillance and evolution of influenza B/Victoria in China since 2019. (A) Influenza B/Victoria sequences from China submitted to GISAID (Shu and McCauley 2017). (B) Proportional clade distribution by month inferred from (A). (C) Regression of root-to-tip genetic distances against sample collection dates. Grey points denote other clades. (D) Time-scaled ML phylogenetic tree of the HA gene of influenza B/Victoria viruses, with viruses from China coloured by clade. Tips with no colour represent global B/Victoria viruses. Amino acid mutations of 3a1 and 3a2 are annotated on ancestral branches. B/Victoria vaccine strains B/Colorado/6/2017 (clade V1A.1), B/Washington/2/2019 (clade V1A.3), and B/Austria/1,359,417/2021 (clade 3a2) are indicated by ‘X’.

A root-to-tip regression of HA distances and sampling dates suggests an increased evolutionary rate for the ancestral branch leading to clade 3a2 viruses (Fig. 1C). Clades 3a1 and 3a2 feature amino acid substitutions indicative of antigenic drift in the HA1 surface protein, including V177I and V220M in clade 3a1 and P144L and K203R in clade 3a2. Notably, clade 3a2 viruses circulating in China contain an additional H122Q mutation that is absent from 3a2 viruses outside China, including the vaccine strain B/Austria/1,359,417/2021, selected in September 2021 (Fig. 1D) (WHO 2021). While HA1 site 122 lies within the highly variable 120-loop (amino acid positions 116–137) (Wang et al. 2008; Shen et al. 2009), position 122 is highly conserved in both influenza B virus lineages: >97 per cent of B/Victoria viruses globally since 2000 contained histidine (H) at site 122, while >97 per cent of B/Yamagata viruses contain glutamine (Q) (Supplementary Data S1).

Furthermore, analysis of all eight gene segments showed that 3a2 viruses circulating outside China were distinct reassortants, having acquired the neuraminidase (NA) and internal genes PB2, PA, MP, and NS from V1A viruses in 2019 (Fig. 2). In contrast, all 3a1 and 3a2 viruses sampled before March 2020 were derived from 3a viruses (shown in light pink), as were the majority of viruses collected in late 2020 and 2021 in China, hinting that no reassortment occurred to drive their circulation there (Supplementary Fig. S4).

Figure 2.

Tanglegram of B/Victoria virus reassortment. Coloured lines connect each virus across all eight genes, showing incongruence between and within major clades. Arrows indicate where clade 3a2 viruses in global circulation acquired NA, PB2, PA, MP, and NS genes from clade V1A viruses.

Temporal phylodynamics of seasonal influenza in China

Root-to-tip regression showed the concatenated whole genomes of circulating clades in China, 3a1 and 3a2, a well as V1A.3-China (the largest monophyletic clade of V1A.3 in China during 2019–2021, as shown in Fig. 1D) contained sufficient temporal signal (r2 > 0.3) for reliable estimation of phylodynamic parameters (Supplementary Data S2 and Supplementary Fig. S5); however, 3a2 sequences collected in 2019 (n = 2) were excluded due to poor temporal signal or incorrect metadata. The mean evolutionary rates of 3a1, 3a2, and V1A.3-China estimated across the genome were 0.5–1.8 × 10–3 nucleotide substitutions per site per year (nt/subs/site/year) (Supplementary Data S2), similar to the long-term evolutionary rates estimated individually for the internal segments (0.14–1.7 × 10–3 nt/subs/site/year) but lower than mean estimates for the HA and NA segments (∼2 × 10–3 nt/subs/site/year) (Dhanasekaran et al. 2015).

Estimation of the time to most recent common ancestor (tMRCA) points to 3a1 circulation in China since mid-July 2019 (mean 19 July 2019; 95 per cent highest posterior density interval (HPD) is 21 June 2019 to 14 August 2019) (Supplementary Fig. S6). Clade 3a1 collected after COVID-19 lockdowns had a mean tMRCA of 21 April 2020 (HPD, 12 March 2020 to 29 May 2020) and derived from 3a1 collected in China during 2019 (Fig. 3 and Supplementary Fig. S6). Clade 3a2 collected after COVID-19 lockdowns had a mean tMRCA of 18 September 2020 (HPD, 20 August 2020 to 14 October 2020) (Fig. 3) and appeared to derive from Chinese samples in 2019; however, uncertainty remains due to high similarity to 3a2 HA collected globally during late 2019 and early 2020.

Figure 3.

Phylodynamic analysis of influenza B/Victoria clades 3a1 and 3a2 in China since April 2020. Maximum clade credibility trees coloured by geographic region. Bars at nodes denote tMRCA. Curves under the trees show change in effective population size (N) with 95 per cent confidence intervals.

Temporal phylogenetics of clades 3a1 and 3a2 also point to regional persistence and intermittent spatial expansion (Fig. 3). We therefore characterized fluctuations in genetic diversity and transmission rates over time using a non-parametric Bayesian Skyline coalescent model to infer effective population sizes (N) (Drummond et al. 2005) and a Bayesian birth–death skyline model to estimate the reproductive rate (R) (Stadler et al. 2013). We found that N of 3a1 increased substantially from June 2020 to December 2020 and fluctuated slightly in 2021, whereas the N of 3a2 steadily increased from October 2020 to August 2021 (Fig. 3).

Transmission rate estimates showed that R of clade 3a1 increased to R = 1.5 (HPD, 0.8–2.3) after the initial COVID-19 lockdown ended (in May 2020), slightly increased again (R = 1.8, HPD, 1.5–2.2) following National Day holidays, and fell below 1 in 2021 (Fig. 4). R of clade 3a2 started from a high of R = 4.0 (HPD, 0.2–9.1) at the end of September 2020, increased slightly following National Day holidays similar to 3a1 viruses, and declined to R = 1.6 (HPD, 0.5–2.5) by November and December 2020. A brief rebound to R = 2.1 (HPD, 1.5–3.0) was observed around Spring Festival (during February–March) and again during summer vacation (July 2021) at R = 1.6 (HPD, 1.0–2.4). For V1A.3-China, N and R peaked in early 2019, gradually decreased, and nearly died out after 2020 (Fig. 4 and Supplementary Fig. S6).

Figure 4.

Effective reproductive rate (R) of B/Victoria 3a1, 3a2, and V1A.3-China lineages. Curves coloured by clades (3a1, 3a2, and V1A.3-China) represent the mean effective reproductive rate (R) with 95 per cent confidence intervals based on influenza B/Victoria positive cases per week reported from China to GISRS (Hay and McCauley 2018). Vertical dashed red line indicates when the nationwide COVID-19 lockdown was lifted in April 2020 (Kraemer et al. 2020). Vertical dashed lines mark the start of National Day (1–7 October 2020), Spring Festival (11–17 February 2021), and summer vacation (June–August 2021) holidays in mainland China.

Phylogeography of B/Victoria in China

Whole genome phylogenetic trees suggest that 3a1 and 3a2 circulated in eastern China prior to widespread detection in December 2020 and March 2021, respectively (Fig. 3). To a lesser extent, clades 3a1 and 3a2 were regionally maintained in the southwest and northwest, respectively, for several months. Of the 30 possible migration routes between the six geographic regions, we identified four routes with definitive evidence and five with sufficient support for B/Victoria clade 3a1 dissemination (Fig. 5 and Supplementary Data S4). We identified seven definitive and five sufficient migration routes for clade 3a2, and three definitive and three sufficient migration routes for the V1A.3-China lineage (Fig. 5 and Supplementary Data S4). Eastern China acted as a major source for 3a1, 3a2, and V1A.3-China epidemics detected nationwide during 2019–2021. Within eastern China, Jiangxi province was identified as a major epicentre and source for recent influenza B/Victoria virus dissemination (Supplementary Data S3). From the east, four migration routes (to central, south, north, and southwest regions) were observed both prior to and during COVID-19-related control measures (Fig. 5 and Supplementary Data S4). Interregional migration most frequently occurred from east to central (38, 11, and 12 migration events for 3a1, 3a2, and V1A.3-China, respectively), east to south (13, 8, and 17, for the three lineages), and east to southwest (10, 6, and 13, respectively) (Fig. 5). However, new migration routes from northwest to southwest (once for 3a1 and 5 times for 3a2) and southwest to central regions (5 times for 3a1 and once for 3a2) were observed for clades 3a1 and 3a2 (Fig. 5). Northern China had the least number of confirmed cases, serving primarily as a regional sink. Notably, a 3a1 lineage from the southwest (Sichuan, Chongqing, Yunnan, and Guizhou provinces) and a 3a2 lineage from the northwest (Gansu and adjacent provinces) acted as important secondary sources for inter-seasonal epidemics throughout 2021 (Fig. 5 and Supplementary Data S4). Overall, while 3a1 during COVID-19 and V1A.3 prior to COVID-19 were mainly driven by eastern China, multiple provinces acted as major sources for 3a2 (Fig. 5).

Figure 5.

Transmission patterns of influenza B/Victoria 3a1, 3a2, and V1A.3-China lineages in China during 2019–2021. The number of transmission routes with definitive support (BF > 100, solid arrows) and sufficient support (BF > 3, dashed arrows) are indicated on the left, and all migration routes within and between regions are shown by directional arrows in the circular flow diagrams on the right.

Altered age distribution of influenza B virus sequenced cases in China

We found that children under 10 years of age accounted for the overwhelming majority of reported cases regardless of subtype, while the age distribution of B/Victoria showed bimodal peaks. Furthermore, an increased number of cases were detected among individuals 20–30 years of age, and fewer cases were detected among 10- to 20-year-olds compared to the age distribution of all sequenced human seasonal influenza cases in China before COVID-19 (2011–2019) (Fig. 6A, Supplementary Fig. S7). During 2020–2021, the proportion of infections in individuals over 10 years of age increased from less than 35 per cent to almost 45 per cent (Fig. 6B). When compared by year of birth, we also found fewer B/Victoria infections among the early 2000s birth cohort (Fig. 6A).

Figure 6.

Age distribution of sequenced influenza B virus cases in China. (A) Distribution of cases during 2011–2019 by age at time of infection (in years) and year of birth. The fraction of cases was calculated relative to all sequenced influenza B cases. (B) B/Victoria age distribution across four seasonal waves of B/Victoria epidemics in China illustrated in Supplementary Fig. S1.

Discussion

In non-pandemic years, seasonal influenza epidemics in mainland China are characterized by the co-circulation of four influenza virus lineages causing short and intense annual winter epidemics in northern provinces and semi-annual or year-round epidemics in central and southern provinces (Yu et al. 2013; Diamond et al. 2022). Since the emergence of SARS-CoV-2 in late 2019, nearly all influenza cases in China have been attributable to B/Victoria lineage viruses. Although A(H3N2), A(H1N1), and B/Yamagata viruses were detected during late 2019 and early 2020, A(H1N1) and A(H3N2) have since been only sporadically detected, and nearby A(H3N2) epidemics in Southeast Asia have not seeded local outbreaks as international borders have remained closed to non-resident travel. Multiple B/Victoria clades survived despite intermittent COVID-19 control measures through persistent regional circulation. However, since late 2020, seasonal influenza epidemics in China have been caused by two B/Victoria clades, 3a1 and 3a2. The resurgence of B/Victoria lineages during 2020–2021 in China can likely be attributed to the recent acquisition of HA amino acid mutations in 3a1 and 3a2 that render them fitter than previous variants through partial immune escape.

The spatiotemporal patterns of seasonal influenza spread have been intensively studied (Cheng et al. 2013; Pollett et al. 2015; Htwe et al. 2019), yet estimation of phylodynamic patterns in individual countries or continents is often hampered by repeat introductions and co-circulation of multiple transmission lineages. As such, the closure to international travel and the isolated circulation of B/Victoria in mainland China during 2020–2021 offered a unique opportunity to record patterns of B/Victoria transmission in the absence of competition from B/Yamagata, A(H1N1), or A(H3N2) subtype viruses. We identified regions, especially the east, central, northwest, and southwest, capable of sustaining local B/Victoria transmission lineages over several months. Jiangxi province in the southeast acted as a major amplifier and source for seasonal influenza B virus transmission across China, while the northwest and southwest regions became secondary epicentres for sustained 3a2 and 3a1 epidemics, respectively. Prolonged control measures in northern provinces with local SARS-CoV-2 transmission (Zhang et al. 2021), including Beijing, Heilongjiang, Xinjiang, Hebei, and Jilin, may have contributed to the absence of influenza in these regions during COVID-19. Temporal phylodynamics suggest increases in N since July 2020, yet outbreaks remained largely undetected by surveillance systems until December 2020. Taken alongside evidence of sustained low-level transmission in multiple provinces, as shown for east, northwest, and southwest China, this indicates that influenza transmission lineages can survive and cryptically circulate under COVID-19-associated NPIs. Moreover, we also found that 3a2 viruses outside China are reassortants containing five segments from V1A viruses that circulated in 2019. Clades V1A, 3a2, and 3a (ancestral to 3a2) were co-circulating in China during 2019 (Fig. 1D) and may explain the origin of the reassortant 3a2 clade, which potentially emerged from China prior to COVID-19.

Contrasting patterns of age distribution for B/Victoria and B/Yamagata have been shown in many countries outside China in the 2010s, whereby a bimodal age distribution was observed in B/Yamagata cases, with greater rates of infection among adults (Dhanasekaran et al. 2015; Caini et al. 2019). Interestingly, based on sequenced samples in China from 2011 to 2019, a bimodal age distribution of B/Victoria cases was observed, analogous to what has been reported for B/Yamagata globally (Caini et al. 2019). The higher disease burden of B/Victoria viruses among adults may have facilitated its persistence and spread. We also found that the cohort born in the early 2000s in China likely has stronger protection against influenza B/Victoria viruses. A modelling study of influenza age distribution from New Zealand showed additional protection against B/Yamagata in individuals immunologically imprinted with B/Yamagata as their first influenza B virus infection (Vieira et al. 2021). Notably, virus genetic records suggest that the B/Victoria lineage emerged in China in the early 1980s (Chen et al. 2007) and, following global circulation for a decade (Rota et al. 1990), remained endemic in China in the 1990s with only rare detections outside Asia until a global resurgence in 2001 (Shaw et al. 2002). We speculate that B/Victoria immune imprinting may have similarly occurred in China, although this requires further analysis and demonstration based on detailed epidemiological records and immunological studies.

Further, our longitudinal comparison of age profiles among influenza viruses in China found a shift towards individuals above age 10 during 2020–2021, which could be attributable to changes in social contact patterns, improvements in surveillance coverage, changes in healthcare-seeking behaviour, and/or loss of lineage competition. Immune profiles and effective contact patterns of a susceptible host population play a key role in viral transmission dynamics, shown previously as social mixing patterns for global movement of influenza lineages (Bedford et al. 2015) and early transmission dynamics of SARS-CoV-2 in China (Kraemer et al. 2020).

Analysis of influenza B viruses collected during the past two decades (Dhanasekaran et al. 2015; Dudas et al. 2015; Langat et al. 2017; Virk et al. 2020) showed that B/Victoria and B/Yamagata maintain distinct HA and polymerase genes (PB2, PB1, and PA) (Dudas et al. 2015) and show distinct patterns in their evolution and epidemiology. We identified a lineage defining HA1 mutation, H122Q, among B/Victoria clade 3a2 viruses in China, which was common among B/Yamagata viruses. Recent studies (Huang et al. 2021; WHO 2021) showed that the majority of B/Victoria viruses circulating in China during 2020–2021 are not well recognized by ferret antisera raised against the 2020–2021 vaccine strain B/Washington/02/2019 (clade V1A.3), which does not possess the H122Q mutation. Likewise, the current clade 3a2 vaccine strain, B/Austria/1359417/2021, lacks the H122Q mutation seen uniquely among 3a2 viruses in China. The 3a2 lineage circulating in China is genetically and antigenically distinct from the B/Victoria viruses circulating in the rest of the world. This could elicit new epidemics when borders reopen, as people outside China will not have immunity to this genotype. Therefore, we recommend the spread of H122Q B/Victoria viruses be carefully monitored and impacts on age profiles measured upon introduction to susceptible populations, particularly amid the potential for reduced vaccine efficacy. Finally, these findings provide cautionary evidence that long-term regional isolation can promote the emergence of new viral variants. Such unique strains could lead to patterns of viral competition and co-circulation or reassort to form new, more transmissible or pathogenic strains with potentially dire consequences.

Our study has the following limitations. (1) Our phylodynamic inference was based on the concatenated whole genome, as there was insufficient temporal signal to analyse gene segments individually. To avoid potential bias introduced due to reassortment, we compared the phylogenetic relationships of each of the segments to confirm that the branching patterns and clades with high support were consistent across the segments. (2) Studies have shown that phylogeographic inference of viral movement can be affected by sampling bias (De Maio et al. 2015; Kalkauskas et al. 2021). We therefore repeated our analysis on a systematically down-sampled dataset to confirm our results. (3) We were unable to investigate the epidemiological dynamics during early to mid-2020 due to a lack of intensive sequence data. (4) All age distribution analyses were based on sequenced samples, which could have introduced sampling bias. (5) Not all reported cases were sequenced, and lineage information was not available for all reported cases.

Overall, the genomic data analysed here provide important insights into the epidemiology of influenza viruses circulating in humans in China during a unique period of prolonged border closure during 2020–2021. In particular, we find that distinct B/Victoria viruses that circulated during the 2019–2020 epidemic maintained low-level circulation in multiple regions within China, despite COVID-19-associated NPIs. Two of these B/Victoria lineages re-emerged in late 2020 to cause outbreaks of increasing intensity throughout 2021, while local circulation of A(H1N1), A(H3N2), and B/Yamagata virus lineages perished. This demonstrates that while COVID-19-associated NPIs could limit the circulation of influenza virus lineages, some low-level circulation may persist and resurge upon relaxation of control measures. We also find that multiple provinces can maintain influenza transmission lineages for several months. Jiangxi province was identified as a major epicentre during 2020–2021, and we show that lineage expansion and interprovincial dispersion coincided with major holidays, which could inform regional and national control policies. While the observed shift in age-susceptibility profile requires further study, clade 3a2 B/Victoria viruses with the HA1 H122Q mutation should be monitored to prevent major outbreaks outside China as travel measures are relaxed.

Click here for additional data file.