Detection of baloxavir resistant influenza A viruses using next generation sequencing and pyrosequencing methods.

Baloxavir, a new antiviral drug targeting cap-dependent endonuclease activity of polymerase acidic (PA) protein of influenza viruses, is now approved in multiple countries. Several substitutions at isoleucine 38 in PA protein (e.g., PA-I38T) have been associated with decreased baloxavir susceptibility in vitro and in vivo. In recent years, next generation sequencing (NGS) analysis and pyrosequencing have been used by CDC and U.S. Public Health Laboratories to monitor drug susceptibility of influenza viruses. Here we described an improved pyrosequencing assay for detecting influenza A viruses carrying substitutions at PA-38. Cyclic and customized orders of nucleotide dispensation were evaluated, and pyrosequencing results were compared to those generated using NGS. Our data showed that the customized nucleotide dispensation has improved the pyrosequencing assay performance in identification of double mixtures (e.g., PA-38I/T); however, identification of PA-38 variants in triple mixtures remains a challenge. While NGS analysis indicated the presence of PA-I38K in one clinical specimen and isolate, our attempts to detect this mutation by pyrosequencing or recover the virus carrying PA-I38K in cell culture were unsuccessful, raising a possibility of a rarely occurring sequencing error. Overall, pyrosequencing provides a convenient means to detect baloxavir resistant influenza viruses when NGS is unavailable or a faster turnaround time is required. Published by Elsevier B.V.

Introduction

A new anti-influenza drug, baloxavir marboxil (baloxavir) was approved for treatment of influenza A and B infections in Japan and the U.S. in 2018, and since then it has entered markets in numerous countries. Baloxavir binds to polymerase acidic (PA) subunit of viral RNA polymerase and hinders the cap-dependent endonuclease activity, which is crucial for generation of capped RNA primers for viral transcription (Noshi et al., 2018). Influenza A viruses have been shown to develop resistance to baloxavir by acquiring substitutions at one of the highly conserved residues in the PA catalytic site - isoleucine 38 (I38X). Influenza A viruses carrying I38T (threonine) or I38S (serine) in PA can be readily selected in cell culture in presence of baloxavir (Chesnokov et al., 2020; Noshi et al., 2018). Notably, natural polymorphisms at PA residue 38 (PA-38) are rare among influenza A viruses (Gubareva et al., 2019; Omoto et al., 2018; Stevaert et al., 2013). Substitutions at PA-38 produce a varying effect on baloxavir susceptibility, depending on the substituted amino acid and virus subtype. Valine (V) does not seem to affect baloxavir susceptibility, while leucine (L) conferred 10-fold reduced susceptibility in A(H1N1)pdm09 and ∼3-fold in A(H3N2) viruses (Gubareva et al., 2019; Noshi et al., 2018; Omoto et al., 2018, CDC unpublished data). The amino acids threonine (T), methionine (M) or phenylalanine (F) confer 10- >50-fold reduced susceptibility to baloxavir (Omoto et al., 2018) and detected in viruses collected from ∼2 to 23% of baloxavir recipients infected with influenza A viruses (Hayden et al., 2018; Hirotsu et al., 2019). Virological surveillance conducted in Japan during 2018–19 season detected viruses carrying PA-38 substitutions with frequency of 2.3% in A(H1N1)pdm09 and 8.0% in A(H3N2) viruses, while no I38-substituted influenza B virus was reported (https://www.niid.go.jp/niid/images/flu/resistance/20191227/dr18-19e20191227-1.pdf). In ferret transmissibility experiments, influenza A and B viruses carrying I38T and I38M were shown to efficiently transmit to naïve ferrets by direct contact and airborne routes (Imai et al., 2020; Jones et al., 2020). Moreover, limited human-to-human transmission of I38T-carrying influenza A viruses have been documented, raising public health concerns (Imai et al., 2020; Takashita et al., 2019a). Substitutions at other PA residues (E23G/K, A37T and E199G) were also shown to decrease baloxavir susceptibility, although to a lesser degree; and these changes are less likely to emerge following drug treatment compared to PA-38 substitutions (I38F/M/N/T) (Ince et al., 2020; Omoto et al., 2018).

One of the goals of virological surveillance is to timely detect emergence and spread of antiviral resistant viruses. Testing for drug resistance has become an integral part of influenza virological surveillance conducted worldwide by laboratories participating in the World Health Organization Global Influenza Surveillance and Response System (WHO-GISRS). Over the years, various genotypic methods have been used to detect markers associated with resistance to M2 blockers and neuraminidase inhibitors by WHO Collaborating Centers and other laboratories of GISRS (Deyde et al., 2010; Nakauchi et al., 2011, https://www.who.int/influenza/gisrs_laboratory/antiviral_susceptibility/en/). One such assay, pyrosequencing (Qiagen platform), was utilized by the Centers for Disease Control and Prevention (CDC) to conduct influenza antiviral resistance testing of viruses collected globally (Bright et al., 2005, 2006; Deyde et al., 2010). Many U.S. Public Health Laboratories (PHLs) have also used pyrosequencing to conduct surveillance at the state level (Storms et al., 2012). Pyrosequencing is a convenient method to detect predefined drug resistance markers in virus genomes. Short sequences in the form of pyrograms can be generated in high throughput manner (96-well format) and rapidly analyzed using built-in PyroMark software. Another salient feature of this technology is the high resolution of single nucleotide polymorphism (SNP ∼5%) (Tsiatis et al., 2010). For the last decade, CDC has been using pyrosequencing to diagnose resistance of influenza viruses to neuraminidase inhibitors as a public health service.

In recent years, CDC has implemented the Sequence First Initiative to advance influenza virological surveillance. The next generation sequencing (NGS) analysis of codon-complete genomes has become the primary method of influenza virus characterization (Jester et al., 2018). Typically, PHLs submit upper respiratory specimens that are influenza-positive and of a known subtype or lineage to CDC. These specimens are subjected to NGS prior to their culturing and phenotypic testing. In addition, large subset of phenotypically characterized virus isolates is also sequenced because virus can mutate during culturing. Since 2014–2015, Illumina MiSeq platform has been utilized for sequencing of influenza viruses submitted to CDC. Genomic sequence analysis has been performed using the CDC-developed IRMA (iterative refinement meta-assembler) approach and other bioinformatic tools (Shepard et al., 2016). The term “NGS analysis” is used throughout the text to describe results obtained using this approach. In collaboration with the U.S. Association of Public Health Laboratories (APHL), CDC implemented NGS at three National Influenza Reference Centers (NIRCs) that have been actively participating in the U.S. national influenza surveillance (Jester et al., 2018). Combined CDC and NIRCs data are used to regularly update reports posted on FluView, including surveillance data on drug susceptibility of influenza viruses circulating in the U.S.

NGS analysis enables detection of all known markers of antiviral resistance in the viral genome, however, this method is not yet widely used by PHLs for surveillance and clinical care testing. Therefore, we investigated the applicability and limitations of pyrosequencing for detecting various amino acid substitutions at PA-38 that have been previously identified in baloxavir-resistant seasonal influenza A viruses.

Material and methods

Viruses

Influenza A viruses used in this study were submitted to the WHO Collaborating Centre for Surveillance, Epidemiology and Control of Influenza at the CDC by U.S. PHLs and other laboratories participating in GISRS (Table S1). A reverse genetically engineered influenza A(H1N1)pdm09 virus containing I38F, was kindly provided by WHO Collaborating Centre for Reference and Research on Influenza at Melbourne (Koszalka et al., 2019). Additionally, A(H1N1)pdm09 viruses carrying either I38S or I38T were selected by culturing A/Illinois/08/2018 (I38) in Madin-Darby canine kidney (MDCK) cells for five passages in presence of baloxavir (Chesnokov et al., 2020). Codon complete genome sequences of both derivative viruses were analyzed by NGS. Besides I38S, no additional substitutions were detected in the I38S carrying virus; while the genome of I38T virus had one additional substitution, HA-K119N (Chesnokov et al., 2020). A(H1N1)pdm09 viruses were propagated in MDCK cells, whereas A(H3N2) viruses were propagated in MDCK-SIAT1 cells (kindly provided by M. Matrosovich). Prior to use in this study, all viruses were propagated by infecting cells at a low multiplicity of infection (MOI) and incubating for short period of time, except two viruses containing I38V, A/California/153/2016 and A/Hawaii/89/2016, were propagated using procedure for routine surveillance to produce high hemagglutination titers (https://www.who.int/influenza/gisrs_laboratory/manual_diagnosis_surveillance_influenza/en/) (Table S1).

Table S1

Influenza A viruses used in this study.

Subtype	Virus name	Amino acid at PA-38	Codona	GISAID # EPI_ISL_	Source	Reference
H1N1pdm09	A/Illinois/08/2018	I	ATT	315855	Surveillance	Gubareva et al., 2019
	RG-A/Perth/261/2009, I38F	F	TTT	N/A	Reverse Genetics	Koszalka et al., 2019
	A/Illinois/37/2018	L	CTT	315856	Surveillance	Gubareva et al., 2019
	A/Illinois/08/2018	S	AGT	365522	In vitro selection	Chesnokov et al., 2020
	A/Illinois/08/2018	T	ACT	348120	In vitro selection	Chesnokov et al., 2020
	A/California/153/2016	V	GTT	241691	Surveillance	Gubareva et al., 2019
H3N2	A/Louisiana/50/2017	I	ATA	315857	Surveillance	Gubareva et al., 2019
	A/Louisiana/49/2017	M	ATG	315858	Surveillance	Gubareva et al., 2019
	A/Bangladesh/3007/2017	T	ACA	286069	Surveillance	Chesnokov et al., 2020
	A/Hawaii/89/2016	V	GTA	239181	Surveillance	Gubareva et al., 2019

N/A: Not available.

Letters in boldface indicate the nucleotide change compared to wildtype codon for isoleucine (I).

To evaluate accuracy of detection and quantitation of viral subpopulations by pyrosequencing, we prepared subtype-specific artificial mixtures of two or three viruses carrying different amino acids at PA-38. The mixtures were prepared using virus stock preparations, and ratios were determined using NGS analysis.

Next generation sequencing

Viral RNA was extracted using QIAamp Viral RNA Mini Kit (Qiagen) and codon-complete influenza genome was amplified using Uni/Inf primer set and Super-Script III One-Step RT-PCR with Platinum Taq High Fidelity enzyme (Invitrogen). Indexed paired end libraries were generated using Nextera XT Sample Preparation Kit (Illumina) following the manufacturer's protocol. Illumina MiSeq was used to generate sequences that were analyzed by the IRMA approach (Shepard et al., 2016). The current cut-off for reportable SNP for depositing NGS data to public databases (GISAID and GenBank) and for routine surveillance purposes is set at 20–25%, although sequence data are stored and can be re-analyzed to provide a more detailed information on viral quasispecies composition. To analyze artificially prepared virus mixtures, RNA extraction, PCR amplification and library preparation were carried out three times independently to evaluate any possible variation that the sequencing process may have on PA-38 variant detection.

Pyrosequencing

Viral RNA was extracted from 100 μl of sample using MagNA Pure Compact Nucleic Acid Isolation Kit I (Roche Diagnostics). A SuperScript III One-Step RT-PCR system with Platinum Taq High Fidelity enzyme (Invitrogen) was used for cDNA synthesis and amplification. Type A-specific primers, forward InfA-F58 and reverse InfA-R280biot, were used for RT-PCR (Table S2). Subtype-specific reverse primers (InfH1-R1141biot and InfH3-R1126biot) were used to generate long amplicons (Table S2). Pyrosequencing reactions were performed on PyroMark Q96 ID (Qiagen) instrument as described by the manufacturer's protocol. Both sequence analysis (SQA) and allele quantitation (AQ) modes of the instrument were utilized for analysis of substitutions at PA-38. Two different nucleotide dispensation orders, cyclic (GCAT)6G and customized GCAAGCTTC(GCAT)4, were used to perform pyrosequencing by SQA mode. SNP analysis by AQ mode was conducted to determine the proportion of each PA-38 variant in a mixture. The dispensation order for AQ analysis was generated by PyroMark software depending on the targeted SNP sequence. A 10% cut-off was used to determine the detection of an SNP (Deyde et al., 2009).

Table S2

Primers used for RT-PCR, pyrosequencing and Sanger sequencing for analysis of PA-38 codon.

Name of primer	Sequence (5’ to 3’)
RT-PCR
InfA -F58	GCAATGAAAGARTATGGGG
InfA-R280biot	biot-TACTGTTYACYACTGTCCAGGCCA
InfH1-R1141biot	biot-AGTCCACTTTTTCTGGTGCCATAT
InfH3-R1126biot	biot-GTGCCATGTTTTCACCAAGAG
Pyrosequencing
InfH1-F91	GAAACTAATAAGTTTGCTGC
InfH3-F95	CCAACAAATTTGCAGC
Sanger sequencing
InfH3-F70	GAAAAAGCAATGAAAGAGT

InfA: influenza type A specific, InfH1: influenza A(H1N1)pdm09 specific, InfH3: influenza A(H3N2) specific, biot: biotinylated at 5’ end.

RNA extraction followed by RT-PCR and pyrosequencing, was carried out three independent times to calculate proportions (%) of each PA-38 variant in artificially prepared virus mixtures. All pyrograms were visually inspected by the operator. Pyrograms were interpreted as: “pass” if variant sequences at PA-38 and downstream sequence were correctly determined; “fail” if PA-38 variant sequences and downstream sequence could not be correctly interpreted; “indeterminant” if PA-38 variant sequences could not be conclusively determined, but the downstream sequence was correctly identified.

Results

Optimization of pyrosequencing to detect PA-38 variants in the double mixtures

A panel of ten seasonal influenza A viruses with and without amino acid substitutions at PA-38 was assembled for this study (Table S1). This panel comprised of viruses detected by surveillance, selected in vitro, or generated using reverse genetics. Viruses were tested using, 1) type A-specific primers for RT-PCR amplification, 2) subtype-specific sequencing primers, and 3) the cyclic order of nucleotide dispensation (GCAT)6G, as described previously (Methods and Table S2) (Chesnokov et al., 2020; Gubareva et al., 2019). With this approach, amino acids I, F, M, L, S, T, or V at PA-38 were readily identified (data not shown).

It is common for respiratory specimens collected from drug-treated patients to contain a mixture of sensitive and resistant viruses. Indeed, viruses with and without amino acid substitution at PA-38, e.g., I38 T/I, were detected in baloxavir-treated patients (Takashita et al., 2019a, 2019b; Uehara et al., 2019). Although rare, virus culturing can also result in emergence of PA polymorphic variants, including at PA-38 (Chesnokov et al., 2020; Gubareva et al., 2019).

We wanted to investigate whether PA-38 variants in the mixtures could be reliably identified using pyrosequencing assay. To this end, we prepared artificial mixtures containing two viruses, one of which had I38 (wildtype) (Table 1 ). Using cyclic dispensation (SQA mode), we tested five different double mixtures of A(H1N1)pdm09 subtype. Notably, only one, containing I/F, was identified correctly after visual inspection of pyrogram and passed the test (Table 1). Pyrograms for other mixtures - I/L, I/S, I/T, or I/V - failed quality assessment due to misalignment of downstream sequence peaks (see a representative pyrogram for I/T in Fig. 1 A) and could not be interpreted by an operator. Therefore, we wanted to know whether the assay performance can be improved by using a customized order of nucleotide dispensation. Based on analysis of PA gene sequence alignment for both the subtypes (not shown) and an approach previously described (Levine et al., 2011), we customized the dispensation by including nine non-cyclic nucleotides prior to four rounds of cyclic dispensations, which should result in synchronized extension of downstream sequence reads for the variants. Addition of A and T at positions four and eight of customized dispensation, respectively, was carried out to capture any unincorporated nucleotides due to homopolymers present in certain variants (e.g., F→TTT, K→AAA) (Zaheer et al., 2012). Thus, the following dispensation order was derived: GCAAGCTTC(GCAT)4. When the single variants were tested using this customized dispensation, all different amino acids at PA-38 were easily identified in both subtypes (representative pyrograms in Fig. 1B and C). For all five double mixtures of A(H1N1)pdm09 tested, visual inspection of pyrograms was required to pass the test and to confirm the presence of both the variants (representative pyrogram in Fig. 1D, Table 1).

Table 1

Pyrosequencing and NGS analysis to detect and quantitate proportions of PA-38 variants in the artificially prepared influenza A virus double mixtures.

Subtype	Virus name and amino acid at PA-38	Codona	PA-38 mixture	Pyrosequencing, SQA resultsb		Percentages of detected variants (Mean ± SD)c
						NGS	Pyrosequencing, AQ resultsd
				Cyclic	Custom		Short amplicon	Long amplicon
H1N1pdm09	A/Illinois/08/2018-I38RG-A/Perth/261/2009-I38F	ATTTTT	I/F	pass	pass	I: 76.7 ± 3.7F: 23.3 ± 3.7	I: 76.3 ± 1.2F: 23.7 ± 1.2	–
	A/Illinois/08/2018-I38A/Illinois/37/2018-I38L	ATTCTT	I/L	fail	pass	I: 40.6 ± 0.5L: 59.4 ± 0.5	I: 45.4 ± 1.2L: 54.6 ± 1.2	–
	A/Illinois/08/2018-I38A/Illinois/08/2018-I38S	ATTAGT	I/S	fail	pass	I: 62.3 ± 2.0S: 37.7 ± 2.0	I: 54.8 ± 2.5S: 45.2 ± 2.5	–
	A/Illinois/08/2018-I38A/Illinois/08/2018-I38T	ATTACT	I/T	fail	pass	I: 44.2 ± 0.1T: 55.8 ± 0.1	I: 47.3 ± 1.3T: 52.7 ± 1.3	–
	A/Illinois/08/2018-I38A/California/153/2016-I38V	ATTGTT	I/V	fail	pass	I: 78.2 ± 1.1V: 21.8 ± 1.1	I: 54.3 ± 5.8V: 45.7 ± 5.8	I: 71.0 ± 0.8V: 29.0 ± 0.8
H3N2	A/Louisiana/50/2017-I38A/Louisiana/49/2017-I38M	ATAATG	I/M	pass	pass	I: 64.3 ± 0.9M: 35.7 ± 0.9	I: 57.7 ± 0.5M: 42.3 ± 0.5	–
	A/Louisiana/50/2017-I38A/Bangladesh/3007/2017-I38T	ATAACA	I/T	pass	pass	I: 68.9 ± 5.5T: 31.1 ± 5.5	I: 76.5 ± 2.2T: 23.5 ± 2.2	–
	A/Louisiana/50/2017-I38A/Hawaii/89/2016-I38V	ATAGTA	I/V	pass	pass	I: 56.4 ± 0.2V: 43.6 ± 0.2	I: 39.2 ± 1.0V: 60.8 ± 1.0	I: 60.0 ± 2.4V: 40.0 ± 2.4

AQ: allele quantitation, RG: reverse genetically engineered virus, SD: standard deviation, -: not tested.

Letters in boldface indicate the nucleotide change compared to wildtype codon for isoleucine (I).

Pyrosequencing was performed in SQA mode using cyclic (GCAT)6G and customized GCAAGCTTC(GCAT)4 nucleotide dispensation orders to detect PA-38 variants in the mixtures. All the pyrograms were visually inspected by the operator: pass = both variants were correctly identified; fail = variant sequences could not be determined.

Results are derived from three independently carried out RNA extraction, RT-PCR and pyrosequencing; and three independent RNA extraction, PCR amplification and library preparation for NGS.

Dispensation order for AQ analysis was automatically generated by PyroMark software depending on the targeted sequences and a position of SNP.

Fig. 1

Representative pyrograms and software readouts for influenza A(H1N1)pdm09 viruses containing one or two PA-38 variants. Pyrosequencing was carried out in SQA mode using either (A) cyclic (GCAT)6G or (B–D) customized GCAAGCTTC(GCAT)4 nucleotide dispensation orders. (D) Arrow points to a nucleotide “C” identified at the second position of the triplet encoding the amino acid at PA-38; indicating the presence of 38 T. (E) Pyrogram of 38I/T mixture by AQ mode with nucleotide percentages provided by the PyroMark ID software (blue box). Expected sequences were based on NGS analysis. Underlined nucleotides indicate codon for PA-38. Lower case nucleotides indicate the single nucleotide polymorphism (C and T) at the second position of codon for PA-38. Highlighted readout is determined by the software: Blue = pass; Yellow = check (visual inspection by the operator is required to interpret the sequence and determine pass or fail); Red = fail.

Next, we tested double mixtures of A(H3N2) viruses using both cyclic and customized nucleotide dispensations. Regardless of dispensation order used, both PA variants in the mixtures were correctly identified (Table 1). Because customized dispensation produced good results for both virus subtypes, it was used in all subsequent experiments.

Comparison of PA-38 variant proportions in the double mixtures using pyrosequencing and NGS analysis

The proportions of PA-38 variants were determined in eight artificial mixtures by both, NGS analysis and pyrosequencing (Table 1). For the latter assay, testing was done using AQ mode, which requires knowledge of targeted sequences and a position of SNP. Of note, NGS and pyrosequencing were performed independently, including RNA extraction and downstream steps. Considering this experimental set-up, the proportions determined by both assays correlated well with differences ranging 0.4–8%, except for two mixtures containing I38V, which showed difference of 17–24% (Table 1). This discrepant outcome required further investigation. Among the viruses used to prepare mixtures, only two viruses carrying I38V were grown differently (see Methods). It is known that infection of cells at high MOI and prolonged incubation can lead to accumulation of defective interfering (DI) particles. Such DI particles often contain defective gene segments encoding polymerase proteins (Dimmock and Easton, 2014). The defective gene segments mostly retain the gene termini but lack the middle region. Notably, NGS method utilizes universal primers to 5′ and 3’ termini of all gene segments, which results in amplification of DI particle sequences. However, such defective PA sequences would be recognized at the NGS assembly step and purged. Conversely, both the full-size and the defective PA gene segments, would be amplified and sequenced by pyrosequencing as DNA amplicon is short (223 bp), and the targeted nucleotide triplet (PA-38) resides close to the gene terminus. Therefore, we wanted to test whether use of a long amplicon would improve the assessment of virus subpopulations by pyrosequencing. To this end, we used new subtype-specific reverse primers, InfH1-R1141biot and InfH3-R1126biot, to generate long amplicons (1084 and 1069bp, respectively) (Table S2). Indeed, sequencing of the long amplicons yielded a decreased proportion of I38V in both mixtures (Table 1), bringing the results of pyrosequencing closer to those of NGS. Of note, use of a short DNA amplicon for pyrosequencing is preferred as it is likely to improve testing results for samples with low virus load.

Analysis of A(H3N2) virus isolates containing mixture of PA-38 variants

Next, using pyrosequencing we tested three A(H3N2) virus isolates containing PA-38 polymorphic variants I/M, I/T, or I/K, as determined by NGS analysis conducted during routine surveillance (Chesnokov et al., 2020; Gubareva et al., 2019). For A/Massachusetts/04/2019 (passage S2), the proportion of I38M variant determined by NGS and pyrosequencing was similar, 27% vs. 25% (Table 2 ). For the second isolate, A/Bangladesh/3007/2017 (passage S2), NGS detected ∼10% lower proportion of I38T variant (39% vs 50%). It is worth noting that I38-substituted viruses identified in both isolates were successfully purified by conventional limiting dilution and tested phenotypically. Substitutions I38M and I38T conferred 11- and 116-fold reduced baloxavir susceptibility, respectively (Gubareva et al., 2019). Notably, neither assay detected I38-substituted variants in respective clinical specimens, suggesting that the PA-38 mutants were selected during virus culturing. For the third isolate, A/Hawaii/28/2017 (passage S2), there was an apparent discordance between the results of NGS and pyrosequencing. While NGS analysis detected I38K variant in both instances, clinical specimen (14%) and the virus isolate (27%); pyrosequencing did not detect it in either sample (Table 2). In addition, presence of I38K in the virus isolate could not be ascertained using Sanger sequencing (Table S2, data not shown). Lysine (K) can be encoded by nucleotide triplets AAG or AAA. For A/Hawaii/28/2017, the nucleotide substitution was detected by NGS analysis at the second position in the triplet, ATA → AAA. We attempted to recover the I38K-containing virus by limiting dilution. A total of 62 biological virus clones were generated from the S2 virus isolate and tested using pyrosequencing, but all clones showed only ATA at PA-38 (wildtype). It was shown that the substitution of a nonpolar isoleucine at 38 with a polar amino acid (T or S) conferred greater baloxavir resistance than with other nonpolar amino acids (F, M, L, or V) (Chesnokov et al., 2020; Omoto et al., 2018). Therefore, positively charged amino acid lysine at PA-38 is likely to adversely affect baloxavir binding to the catalytic site and confer drug resistance. Therefore, we cultured the S2 isolate containing 27% I38K (according to NGS) in presence of baloxavir (0.055–4.5 nM) in an attempt to increase the proportion of I38K in virus population. However, this approach did not yield such enrichment. Because we could not find biological evidence for the presence of I38K variant, we concluded its detection by NGS analysis is likely to be erroneous. It is known that data obtained by different sequencing-by-synthesis technologies may contain biases or errors depending on the experimental procedures (Schirmer et al., 2016). One of such common challenges is interpretation of sequences containing a long stretch of homonucleotides (e.g., AAAA). However, detection of I38K in both clinical specimen and virus isolate remains puzzling.

Table 2

PA-I38-substituted variant proportions in influenza A(H3N2) virus samples determined using NGS analysis and pyrosequencing.

Virus name	PA GISAID accession #	Sample	Percentage of PA-I38-substituted variant (Mean ± SD)
			NGS	Pyrosequencing (AQ)
A/Massachusetts/04/2019	EPI1362029	clinical	ND	ND
	EPI1397514	S2 isolate	M (27.0)	M (24.8 ± 0.7)
A/Bangladesh/3007/2017	EPI1086149	clinical	ND	ND
	EPI1107334	S2 isolate	T (39.2)	T (49.7 ± 1.2)
A/Hawaii/28/2017	EPI1016323	clinical	K (14)	ND
	EPI1039207	S2 isolate	K (27)	ND

ND: not detected (below limit of detection), S2: passage 2 in MDCK-SIAT1 cells.

ND: not detected (below limit of detection), S2: passage 2 in lLine">MDCK-SIAT1 cells.

Detection of PA-38 variants in triple mixtures using pyrosequencing

Reportedly, virus samples collected from baloxavir-treated patients occasionally contain three PA-38 variants (Takashita et al., 2019a, 2019b; Uehara et al., 2019). Therefore, we prepared artificial virus mixtures containing three different PA variants and tested them using NGS and pyrosequencing (Table S3). In three instances, presence of all PA-38 variants was correctly identified upon visual inspection of pyrograms, which were deemed “pass”. For two other mixtures, not all PA-38 variant sequences could be conclusively determined and were deemed “indeterminant” (Table S3). Due to consecutive changes at more than two nucleotides of the same codon, built-in algorithm of the PyroMark software did not support SNP analysis of triple mixtures tested here using AQ mode. Overall, it was more challenging to determine PA-38 variants in triple mixtures of A(H1N1)pdm09 subtype. While all variants in the A(H1N1)pdm09 triple mixture containing I/S/T were identified, the other two mixtures (I/L/T and I/F/T) could not be resolved using pyrosequencing method described here. Visual inspection of pyrogram of the triple mixture I/L/T suggested the presence of ATT and ACT, while presence of CTT in the mixture was not evident. Similarly, it was difficult to ascertain the presence of all three PA-38 variants in the mixture I/F/T (data not shown).

Table S3

Assessment of PA-38 variants in artificially prepared influenza A virus mixtures.

Subtype	Virus name and amino acid at PA-38	PA-38 mixture	Codona	Percentages of variants by NGS (Mean ± SD)b	Pyrosequencing, SQA resultsc
H1N1pdm09	A/Illinois/08/2018-I38A/Illinois/37/2018-I38LA/Illinois/08/2018-I38T	I/L/T	ATTCTTACT	I: 47.9 ± 1.6L: 28.4 ± 0.8T: 23.7 ± 1.0	indeterminant
	A/Illinois/08/2018-I38RG-A/Perth/261/2009-I38FA/Illinois/08/2018-I38T	I/F/T	ATTTTTACT	I: 26.8 ± 0.7F: 8.0 ± 1.1T: 65.2 ± 1.7	indeterminant
	A/Illinois/08/2018-I38A/Illinois/08/2018-I38SA/Illinois/08/2018-I38T	I/S/T	ATTAGTACT	I: 57.7 ± 0.1S: 20.2 ± 0.7T: 22.1 ± 0.8	Pass (I/S/T)
H3N2	A/Louisiana/50/2017-I38A/Louisiana/49/2017-I38MA/Bangladesh/3007/2017-I38T	I/M/T	ATAAGTACA	I: 50.9 ± 1.2M: 22.9 ± 0.9T: 26.1 ± 2.0	Pass (I/M/T)
	A/Hawaii/89/2016-I38VA/Louisiana/49/2017-I38MA/Bangladesh/3007/2017-I38T	V/M/T	GTAATGACA	V: 21.4 ± 1.0M: 40.2 ± 2.7T: 38.3 ± 2.0	Pass (V/M/T)

RG: reverse genetically engineered virus.

Letters in boldface indicate the nucleotide change compared to wildtype codon for isoleucine (I).

Results are derived from three independent RNA extraction, PCR amplification and library preparation.

Pyrosequencing was performed in SQA mode using customized nucleotide dispensation order GCAAGCTTC(GCAT)4 to detect variants at PA-38 in the mixtures. All the pyrograms were visually inspected by the operator: pass = all variants were correctly identified; indeterminant = all variants could not be conclusively identified.

Discussion

Over the past few years, NGS analysis has become the primary method of influenza virological surveillance in the U.S., including drug resistance monitoring (Jester et al., 2018). Availability of viral genome sequences is especially desirable as it facilitates the establishment of molecular markers associated with decreased susceptibility for new drugs. Conversely, pyrosequencing assay is useful for detecting previously established markers of drug resistance by generating short, targeted sequence readouts. Once optimized, pyrosequencing can provide a formidable method to test hundreds of samples in a timely manner. Here we showed that pyrosequencing using both cyclic and customized nucleotide dispensation orders was suitable for identifying different amino acid substitutions at PA-38 of seasonal influenza A viruses. However, the cyclic dispensation was not optimal for detecting I38-substituted A(H1N1)pdm09 viruses present in a mixture with the wildtype (I38) virus. On the other hand, the customized dispensation improved the assay performance in detecting PA-38 variants in double mixtures. The improvement, however, was not sufficient to identify all variants in triple mixtures. Because double or triple mixtures are anticipated in specimens collected from baloxavir-treated patients, it is prudent not to rely on software sequence readout, but to verify pyrogram by visual inspection. In some instances of triple mixtures, even visual inspection of pyrograms could not ascertain the presence of all the PA-38 variants and therefore the result of these mixtures was deemed “indeterminant”. This limitation of the assay may not be critical for surveillance purposes, but it could complicate the interpretation of testing outcomes needed to make informed clinical care decisions. Moreover, pyrosequencing has an inherent limitation due to use of dATP analog deoxyadenosine α-thio triphosphate in the reaction as it produces a higher peak on the pyrogram than the other dNTPs (Harrington et al., 2013). This needs to be kept in mind, when visually inspecting pyrograms.

Koszalka and colleagues recently published a study, where they explored utility of pyrosequencing for rapid detection of PA-I38X variants in influenza A and B viruses (Koszalka et al., 2020). There were several differences in the assay design compared to our study (different cyclic dispensation order, biotinylation of forward primer, generation of shorter amplicon (100bp), and use of DNA plasmids to prepare artificial PA-38 variant mixtures). The authors showed that pyrosequencing using AQ mode could estimate proportions of variants in known artificial mixtures containing I/T and I/M, but not I/F; and triple mixtures were not tested. Because of the intuitive nature of the pyrosequencing assay, it is possible to create hypothetical scenarios to determine the outcome of a sequence readout. Applying such an approach made us question whether testing unknown samples carrying two to three different variants would yield conclusive data. Taken together, we believe that a customized dispensation order of nucleotides provides an improvement in resolving samples carrying mixtures of PA-38 variants compared to cyclic dispensation orders (Gubareva et al., 2019; Koszalka et al., 2020). Our data indicate that DI particles containing defective PA gene segments may affect the proportion of a viral population detected by pyrosequencing compared to NGS analysis, at least in artificially prepared mixtures. However, this is unlikely to affect testing of clinical specimens. One limitation of our study is not propagating I38V containing viruses using a low MOI and incubating for short period of time as was done for other viruses.

Every sequencing technique has limitations and may introduce biases or errors (e.g., general increase of errors toward the end of reads for Illumina) (Schirmer et al., 2016). In our study, NGS analysis detected the presence of I38K variant in A/Hawaii/28/2017 (H3N2) virus; while pyrosequencing, Sanger sequencing and limiting dilution did not. To the best of our knowledge, viruses carrying I38K as a dominant variant have not been reported by global surveillance or detected after baloxavir exposure. However, virus specimen from a baloxavir-treated patient was reported to contain I38K as a mixture (I38 T/K/I) in Japan during 2018–19 season (Takashita et al., 2019a). Recently Jones and colleagues have shown that introduction of polar, positively charged lysine (K) at PA-38 reduced the polymerase activity by 50% compared to wildtype (I38) in minireplicon assay (Jones et al., 2020). However, they apparently did not investigate the effect of this amino acid substitution on virus properties. To improve the current knowledge, it would be desirable to generate a reverse genetically modified virus and determine effect of I38K on virus fitness and baloxavir susceptibility. This information would benefit the interpretation of PA-38 sequencing results for clinical care use. In future experiments, viruses with recently reported PA-38 substitutions (e.g., I38N) could be derived using reverse genetics and tested by our assay to further show its utility (Ince et al., 2020; Imai et al., 2020).

It is worth noting that the frequency of baloxavir resistance detection among circulating viruses remains low in the U.S. (https://www.cdc.gov/flu/weekly/weeklyarchives2019-2020/Week13.htm). However, PA-38 substitutions were detected in 2.3% of A(H1N1)pdm09 and 8.0% of A(H3N2) viruses collected during 2018–19 season of Japan, with the majority of mutant viruses being isolated from baloxavir-treated children (https://www.niid.go.jp/niid/images/flu/resistance/20191227/dr18-19e20191227-1.pdf,Takashita et al., 2019a). The widespread use of influenza antiviral can lead to emergence of resistant viruses and is a serious public health concern.

In the past, it was suggested that the widespread use of adamantanes, incited by the fear of ‘bird flu’ and the SARS coronavirus facilitated emergence of adamantane-resistant viruses in Asia in 2003–2004 (Bright et al., 2005) and their global spread in the following years (Bright et al., 2006; Deyde et al., 2007). Therefore, it is conceivable that severe influenza seasons and ongoing SARS-CoV-2 pandemic can also lead to rampant use of available antivirals, including baloxavir, which in turn may cause the emergence and spread of drug-resistant viruses. In view of this possibility, it would be prudent to strengthen laboratory capabilities to allow early detection of drug resistance emergence.

In conclusion, various genotypic assays have been used to detect drug resistance markers in influenza viruses, including baloxavir resistant viruses (Gubareva et al., 2019; Koszalka et al., 2020; Nakauchi et al., 2020). All methods have their advantages and disadvantages and decision which assay to use depends on the community needs and laboratory capabilities.

Declaration of competing interest

None.