An 81-year-old male had a slight cough and chest distress on 1/12/2016 at his home in Xihu District, Hangzhou City, Zhejiang Province, China. On the morning of 1/15/2016, the symptoms worsened with a nasty cough, chest distress and fever. That afternoon, the man went to the local community hospital, where the temperature was 38.5 C, and then he was sent to the Hangzhou First People’s Hospital for medical treatment. The examination showed that the white blood cell count (WBC) was 13.4 × 10⁹/L, the percentage of neutrophils (N%) was 92.2%, and the C-reactive protein (CRP) was 110 mg/L. The chest radiographs showed an infection of the right lower lung. The mezlocillin sodium and sulbactam sodium were given for intravenous injection as an anti-infective therapy. The patient was sent to the respiratory department for hospital treatment and pneumonia was confirmed on 1/16/2016.

The next day, a PCR result from a throat swab was positive on influenza A virus and the patient was sent to infection ward for further treatment which included oseltamivir (75mg/bid) and meropenem drugs. On 1/19/2016, patient’s symptoms worsened further and chest radiographs confirmed that infections spread on both lungs. The patient received endotracheal intubation and then admitted to intensive care unit (ICU)The treatment was continued, the dose of oseltamivir was doubled (150mg/bid) and meropenem was changed to imipenem. On the same day, the patient’s RT-PCR test taken from the throat was positive on the M gene, H7 gene, N9 gene and H1 gene of influenza A virus. Next day, the patient was transferred to Hangzhou Xixi Hospital for treatment in isolation where the anti-viral oseltamivir (150mg/bid) continued to be given until 2/20/2016. Although, the symptomatic treatments such as diuresis, analgesia, vasodilation and nutritional support has been given in the Xixi hospital, the patient did not show signs of improvement, and passed away on 2/28/2016.

In the patient’s anamnesis it is stated that the patient went to local live poultry market and bought a live duck at a merchant’s site about a week before the first symptoms appeared on 1/12/2016. The live duck was slaughtered, depilated, and bellied by the merchant at his site. After returning home, the patient salted the duck. The patient had a history of hypertension and denied the history of diabetes, viral hepatitis, tuberculosis, and other diseases. There is no history of trauma, surgery, or blood transfusion. Denying any history of drug or food allergies.

Ten serial samples were collected from this patient from 1/19/2016 to 2/19/2016, including seven throat swabs and three sputa. The swabs and sputa were placed into 1 ml viral transport medium, transported to the laboratory within 24 h at 4°C, and then frozen at −80°C. Viral RNAs were extracted from samples using a RNeasy Mini Kit (QIAGEN, Germany). Identification of influenza A virus was achieved by RT-PCR using specific primers targeting the M, H7, N9, and H1 gene according to the protocol provided by WHO manual (Organization, W.H, 2002). This study was approved by the Institutional Review Board of BGI (NO.BGI-IRB 16008).

Top eight samples were taken to perform single-molecule real-time sequencing (SMRT). The cDNAs were synthesized from viral RNAs by reverse transcription using Uni12 and Uni13 primers (Bi et al., 2016). The PCR was performed using a Phusion High-Fidelity PCR Kit (New England Biolabs) utilizing the barcoded influenza A virus general primers (Supplementary Table S1) (Mei et al., 2016). The concentration of PCR product was quantified by the Agilent Technologies 2,100 bioanalyzer. The two corresponding volumes of PCR products (containing equal mass of dsDNA) were mixed into one sample and quantified in the bioanalyzer again. About 2–3 μg mixed sample was used to SMRTbell library construction following the 2 kb template preparation protocol (Roberts et al., 2013b). The sequencing was performed on a PacBio RS II instrument (Pacific Biosciences, USA) with one SMRT Cell used for each library, using P6/C4 chemistry with a 4 h movie (Bull et al., 2016). SMRTbell adapter sequences were removed and circular consensus sequence (CCS) reads were achieved with SMRT Analysis v2.3 (Roberts et al., 2013a).

The raw CCS reads were filtered by removing low quality reads (length<800bp, passes<5 or estimated accuracy<99.9%). The 800 bp length near the lower limit of influenza A virus genomic RNAs was used to exclude non-full-length genomic sequences. The other two criteria ensured reads with at least 99.9% estimated accuracy and necessary passes (Korlach, 2015). The sequencing error bases (frequency<0.3%) were additionally corrected to improve sequence reliability as follows: First, the remaining sequences after filtration were split into corresponding samples by 100% base match with barcodes. Then, the sequences of one sample were grouped by subtypes and genomic segments according to the sequence annotation result against an influenza virus genomes database downloaded from NCBI (https://ftp.ncbi.nih.gov/genomes/INFLUENZA) using BLASR (v5.1 with options: -bestn 1) (Chaisson and Tesler, 2012). All full-length genomic sequences of one group were aligned end to end using MUSCLE (v3.8.31 with default options) (Edgar, 2004). Following, the number and percentage of base A/C/T/G in the same nucleotide locus of genomic sequences were stated. Finally, the very low frequency base which percentage was less than 0.3% in the number of four type bases of the same nucleotide locus, was replaced with the dominant type of base with the largest proportion in this nucleotide locus.

The diversity index (Shannon entropy) of one group sequences is calculated by the formula (Crooks and Brenner, 2004): S = − 100 ∗ ∑ i = 1 n P i ∗ log 2 P i In which S is the Shannon entropy and P i is the ration of the number of one type of sequence to the number of total types of sequence in one group.

In order to screen out nucleotide loci with high-diversity base composition, we stated the number and percentage of base A/C/T/G in one nucleotide locus and screened out the loci in which the percentages of at least two types of bases were more than 10%.

The clinical symptoms and therapeutic schedule of this patient were recorded in Table 1. Ten samples were collected from this patient on different days, including seven throat swabs (S1-4, S7-10) and three sputa (S5-7). The collected date, sample types and Ct value of RT-PCR for H1N1 and H7N9 of each sample were listed in Table 2. The top eight samples (S1-8) were performed using SMRT with four SMRT cells (S9 and S10 were RT-PCR negative for influenza A virus). A total of 142,496 CCS reads (221.72 Mb) were generated from four SMRT cells, of which 82,471 high quality reads (≥99.9% estimated accuracy) were selected for further analysis according to strict filtering criteria. The IAV mutation rate was about 0.018–0.025%, which means that each replicated influenza genome (∼13 kb) contained an average of 2–3 mutations (Pauly et al., 2017). However, the estimated base sequencing error rate of high-quality CCS reads (∼0.1%) in this study, was notable higher than the normal replication mutation rate (0.018–0.025%) of IAV genomes. Therefore, it was necessary to correct sequencing error bases in the prevention of taking them for real mutations. Finally, 69 group sequences were clustered from four SMRT cells according to samples, subtypes, and genomic segments (Figure 1A). Among them, 58 group sequences were with a satisfactory abundance (the sequences number ≥20). In case that the base type (A/C/G/T) in one nucleotide locus of one group sequences was less than 0.3%, it was considered as sequencing error base. The base sequencing error rates of 58 group sequences ranged from 0.13 to 0.21%, approximate to the estimated sequencing error rate of 0.1%, which are composed of 78.37% mismatches, 13.78% deletions and 7.16% insertions (Figures 1B,C and Supplementary Table S2). Thus, we set 0.3% as the cutoff value to screen out sequencing error base types in nucleotide loci of each group sequences and correct them with dominant base types in these loci.

All the influenza virus genomic sequences produced by SMRT were clustered into 69 groups by eight samples, two subtypes, and eight genomic segments (Figure 1A). To monitor the composition and dynamic of H1N1 and H7N9 genomic sequences, we calculated the number of sequence reads, the number of sequence types, and the diversity index of sequence types in each group (Figure 2). The 16 groups (eight from H1N1 and eight from H7N9) in the first sample S1 confirmed that this patient was coinfected with H1N1 and H7N9 IAV. Interestingly, in sample S1, although the number of H1N1 and H7N9 sequence reads were almost equal, the diversity index of H7N9 sequence types was obviously higher than that of H1N1 (Figure 2). The reasons for the diversity difference of sequence types between H1N1 and H7N9 were unclear. One possible explanation was that in the upper respiratory tract (URT) environment with the dominant α-2,6-SA receptors preferentially recognized by H1N1 virus; the H7N9 virus had to generate more various genomic sequences to adapt to the relative hostile environment (Chen et al., 2013).

Furtherly, there was a sharp decrease of H7N9 viral load in the URT samples from S1 to S4. But the H7N9 viral load was still relative abundant in the subsequent LRT samples from S5 to S7 (Table 2 and Figure 2B). This might indicate that the H7N9 virus had transferred to LRT environment from URT (Gao et al., 2013). Meanwhile, there was an obvious increase of H1N1 viral load in the URT samples from S1 to S2 (Figure 2A). The reason for this increase might be due to the transfer of H7N9 from URT to LRT that contributes freeing up more cellular resources for H1N1 growth in URT environment.

Fifteen and eleven high diversity nucleotide loci were detected in H1N1 and H7N9 genomes, respectively, in which the A-G transversion was the most abundant substitution type (67% in H1N1 and 64% in H7N9) (Table 3). In H1N1, another two substitutions were C-T and A-C transversion (20 and 13% respectively). IN H7N9, the C-T, A-C and G-T took the same 0.09% proportion respectively (Table 3). In terms of amino acid, the percentage of non-synonymous mutation was greater than synonymous mutation, especially the percentage of non-synonymous mutation was up to 91% (10/11) in H7N9, comparing with the non-synonymous mutation of 66% (10/15) in H1N1(Table 3). In H1N1, all five synonymous mutations are in the internal genes of IAV, including two in NS gene (R78R and L69L), one in PB2 gene (T25T), one in NP gene (E64E) and one in M gene (R134R). In H7N9, the only one synonymous mutation was R753R in internal PB2 gene (Table 3). The H7N9 with a high frequency of non-synonymous mutation might indicate that H7N9 as avian-origin influenza virus was subjected to higher selection pressures in the human host (Wei et al., 2012; Xu et al., 2019). It is noteworthy to mention that among of these high-diversity loci, four mutations in H1N1 were involved in the evolution or viral drug-resistance, such as the R222K in the HA protein and the H275Y, A204T and K207R in NA protein (Table 3). In H7N9, three mutations were related to host adaption or viral drug-resistance, including the G685R and E627K in PB2 protein, and K432E in NA protein (Table 3).

We achieved 25 and 22 haplotypes for H1N1 and H7N9 genomic segments in all eight samples based on combining high diversity nucleotide loci in the same genomic sequences, respectively. For instance, we obtained seven haplotypes of NA in H1N1 based on four high diversity nucleotide loci (Table 4). The bases on high diversity nucleotide loci of each haplotype and its abundance in each sample were listed in Supplementary Table S4. Then the composition and dynamics of viral quasispecies in this patient co-infected with H1N1 and H7N9 IAV were displayed along the clinical treatment process (Figure 3). In NA gene, the haplotype Hap_2 replaced Hap_1 as the dominant haplotype with more than half proportion (53.11%,1,119/2,107) in sample S2, which contain H275Y drug-resistant mutation. However, Hap_2 in NA failed to be continuously dominant and Hap_1 return the dominant haplotype in following samples (S3 to S6) (Figure 3A). Similarly, we found Hap_3 of HA with antigenic drift mutation R222K transiently become dominant in sample S6 (56.19%, 127/226) (Figure 3A). Compared with the genomic segment integrity of H1N1 viral quasispecies among almost all samples, The H7N9 quasispecies had a poor integrity except for the first sample S1 (Figure 3B). In conclusion, the haplotype provides a more concise overview of viral quasispecies composition and dynamics than whole genomic sequences (Figure 2 and Figure 3).

The data that support the findings of this study have been deposited into CNGB Sequence Archive (Guo et al., 2020)of CNGBdb (Chen et al., 2020)with accession number CNP0001131.