Virome Sequencing Identifies H5N1 Avian Influenza in Wastewater from Nine Cities

Sample collection and processing

All collection, processing, and sequencing methods have been described in an earlier publication¹⁷ and were followed precisely. To restate, between 100–500 mL of raw wastewater influent was collected into 500 ml leak-proof prelabeled sample bottles at 23 wastewater treatment plants in 10 Texas cities, the identities of these cities are masked on the request of local public health officials and water utilities. From South Texas City H, sewage from three treatment plants is combined into a single sample before shipment to the sequencing laboratory; all other sites had samples processed individually. Treatment plants were coded upon the request of public health officials. The surface of the sample bottles was decontaminated with 10% bleach and moved to a “clean” zone, where the samples were sealed into biohazard bags in shipping boxes with absorbent pads and ice packs for overnight shipping to the Alkek Center for Metagenomics and Microbiome Research at Baylor College of Medicine, Houston, TX.

Wastewater samples were barcoded upon arrival and stored at 4 °C until processing. First, 50 mL of wastewater was decanted and centrifuged at 3374 x g for 10 min, separating the solid and liquid fractions. The supernatant was then vacuum filtered using an ion-based cellulose filter paper and the virus-containing cellulose filter was placed into a bead-beating tube with lysis buffer. The tube was run on a homogenizer for 1 min at 5 m/s, rested for 1 min, then run on the homogenizer for 1 more minute. Following bead beating the samples were centrifuged at 14– 17×1000 RPM for 2 min. DNA and RNA were extracted using the Qiagen QIAamp VIRAL RNA Mini Kit.

Library preparation, probe-based virome capture and sequencing

RNA extracts were converted to cDNA using Protoscript II First Strand cDNA Synthesis Kit (New England Biolabs Inc.), NEBNext Ultra II Non-Directional RNA Second Strand Module (New England Biolabs Inc.), and Random Primer 6 (New England Biolabs Inc.). A total of 25 ng of the cDNA and DNA mix was used for library construction using Twist Library Preparation EF 2.0 Kit and Twist Universal Adaptor System (Twist Biosciences). The libraries were pooled, a maximum of 18 samples per pool, at equal mass to a total 1,500 ng per pool. The Twist Comprehensive Viral Research Panel (Twist Biosciences) was used to hybridize the probes at 70 °C for 16 h. The post-capture pool was further PCR amplified for 12 cycles and final libraries were sequenced on Illumina NovaSeq 6000 SP flow cell, to generate 2×150 bp paired-end reads. Following sequencing, raw data files in binary base call (BCL) format were converted into FASTQs and demultiplexed based on the dual-index barcodes using the Illumina ‘bcl2fastq’ software.

Influenza A read detection and serotype assignment

Reads were processed and analyzed using the publicly available iav_serotype (v0.1.1) software package (https://github.com/mtisza1/influenza_a_serotype). For context, iav_serotype utilizes a database of 207,488 publicly available, complete Influenza A segments from all known serotypes (Influenza_A_segment_sequences database v1.1, available on Zenodo). In detail, paired-end reads from wastewater virome sequencing are quality filtered using fastp (v0.23.4) ¹⁹ with default settings. Reads are aligned to Influenza_A_segment_sequences database v1.1 using minimap2²⁰ (v2.28-r1209) and flags “-cx sr --secondary=yes” to allow secondary alignments. Average nucleotide identity (ANI) and alignment fraction (AF) of each alignment are calculated. Then read alignments are assigned to a particular serotype if the best alignment is exclusive to one serotype and (ANI*AF >= 0.9). Else, reads are assigned as “ambiguous”.

Manual Inspection of Putative H5N1 Sequences

After automated (iav_serotype) detection of H5N1 in wastewater samples was noted by MJT, sequencing reads were shared for independent validation by co-authors BMH and JRC. After sequence inspection confirmed the general findings, parameters for final validation were agreed upon. Using web BLAST^21,22 against all of NT (final inspection on May 4^th and 5^th), forward and reverse reads from each pair were inspected separately. PCR duplicates were assumed to share origin. Read pairs were marked as “H5N1 specific”, otherwise marked as “ambiguous”. “H5N1 specific” label was used only if neither forward or reverse read alignments had a higher bit score to non-H5N1 genomes and either of the following conditions were met: 1) forward and/or reverse reads aligned with 100% identity and alignment fraction to H5N1genomes and alignments to other serotype were worse by 1 or more bases. 2) forward and/or reverse reads aligned with less than 100% identity to H5N1 and non-H5N1 alignments were worse by 2 or more bases. Of 49 sequencing libraries with automatically assigned reads (per iav_serotype), 46 libraries were manually validated as containing H5N1 specific reads.

Variant Analysis of Pooled Reads

The variant analysis was conducted on pooled reads, independent of the serotype assignment. Here, raw reads underwent quality trimming to a score of 20 and reads with fewer than 30 base pairs were discarded using BBDuk (v38.84)²³. Duplicate reads were removed using Dedupe (v38.84). The quality-controlled reads were then pooled and aligned to the GISAID EPI_ISL_19027114 lineage using the Geneious Assembler in Geneious Prime 2024.0.5, with medium sensitivity settings to detect variations of any size.

We excluded variant calls with a Strand Bias P-value less than 0.05 and those resulting from ambiguities in the reference sequence. Strand Bias P-value assesses whether observed strand bias is likely due to chance or systematic effects. Variants with significant strand bias were filtered out to maintain data quality.

Influenza Hospitalization Data

Weekly statewide and regional hospitalization data was downloaded directly from the Texas Respiratory Illness Interactive Dashboard (https://texas-respiratory-illness-dashboard-txdshsea.hub.arcgis.com/) on May 4^th, 2024.