Influenza A virus infection in turkeys induces respiratory and enteric bacterial dysbiosis correlating with cytokine gene expression

Turkeys and ethics statement

Specific pathogen free (SPF) turkeys were obtained from the Food Animal Health Research Program (OARDC, OSU) flock (Ngunjiri et al., 2019a) and maintained free from known turkey pathogens including reovirus, avian encephalomyelitis virus, infectious bursal disease virus, infectious bronchitis virus, Newcastle disease virus, avian influenza virus, hemorrhagic enteritis virus, Mycoplasma gallisepticum, M. synoviae, M. meleagridis, Bordetella avium and Salmonella pullorum. The turkeys were housed and acclimatized in a BSL-2 facility with forced air ventilation and adequate air exchanges to prevent ammonia build up. Air entering or leaving the facility is HEPA filtered. At 2 weeks of age, birds were randomly picked from the brooders, assigned to one of the experimental groups, tagged, transferred to isolators (Model 934–1, 900 sq. inch, Federal Designs Inc., Comer, GA) inside the BSL-2 facility, and experiments were conducted in accordance with protocol #2009AG0002-R2 approved by The Ohio State University Institutional Animal Care and Use Committee (IACUC). Room and isolator temperatures were maintained at 25 ± 3 °C. The animals were provided with ad libitum access to feed and water and their welfare monitored twice daily. They were humanely euthanized at 5 (1st sampling time point) or 14 days post-infection (2nd sampling time-point and experiment endpoint), or when they displayed symptoms such as ruffled feathers and reluctance to move, not moving when prodded, respiratory distress, or injuries that were not related to experimental treatment. Administration of analgesia drugs was unnecessary since the birds were euthanized soon after the emergence of clinical symptoms. Humane euthanasia was actualized by exposure to carbon dioxide (CO₂) as previously described (Elaish et al., 2017; Elaish et al., 2019; Ghorbani et al., 2019; Jang et al., 2018). Three poults were placed in the euthanasia chamber connected to a CO₂ source, with CO₂ flow set at 10–30% displacement of chamber volume/minute. Birds were observed for respiratory arrest and lack of heartbeat and the CO₂ flow was maintained for at least one minute after the respiratory arrest and lack of heartbeat were observed. Upon observation of death, an additional secondary physical euthanasia (cervical dislocation or removal of a vital organ) was performed before collection of tissues and carcass disposal (Elaish et al., 2017; Elaish et al., 2019; Ghorbani et al., 2019; Jang et al., 2018).

Experimental design

A total of 96 poults were placed in six isolators at 10 days of age. Each bird was randomly picked from the brooder and assigned to one of isolators. The first bird was assigned to Isolator #1, second bird to isolator #2, and so on, and the sequence was repeated until 16 birds were placed in each isolator. After an 8-day acclimatization period, the birds were assigned into three treatment groups (2 isolators/group; 32 birds/group) and inoculated with PBS-only (“Mock”) or H5N2 low-pathogenicity avian influenza viruses: CK/PA/13609/93 (CKPA) and TK/MN/10734-2/95 (TKMN). The treatment group was the experimental unit. Sixteen birds per group were sufficient for statistical comparisons of viral shedding, gene expression, and microbiota data based on our published studies (Ghorbani et al., 2020; Jang et al., 2018; Johnson et al., 2018; Ngunjiri et al., 2019b; Taylor et al., 2020). Median egg infectious dose (EID₅₀) for each virus stock was calculated prior to normalization of virus concentration on the day of challenge. Each bird was intrachoanally inoculated with 200 µL of PBS or virus (10⁶ EID₅₀). At 5 days post infection (dpi), 8 birds from each isolator were euthanized for necropsy. At 14 dpi, the remaining birds were necropsied. Five birds in the CKPA-treated group developed severe clinical signs and were euthanized at 4 dpi and between 6 and 7 dpi. These birds were not sampled for microbiota because our target sampling time-points were at 5 and 14 dpi. Animal handling procedures were done starting from the Mock group, then the infected groups to reduce confounders. All personnel involved in this study were aware of group allocation of the animals. The research question, experimental design, and data analysis protocols were not registered before the study.

Quantification of hemagglutination inhibition (HI) antibodies

To determine homologous serum HI antibody titers, blood was drawn from live birds via the brachial (wing) vein for serum collection prior to euthanasia at 14 dpi. The blood was kept overnight at room temperature for serum separation to occur. The separated sera were decanted into two mL tubes, inactivated at 60 °C for 30 min, and stored at −20 °C until used for HI test. The HI test was conducted as described previously (Jang et al., 2018). In brief: 50 µl of heat-inactivated serum was serially diluted in 96-well plates and mixed with an equal volume of homologous virus preparation containing eight hemagglutinating units of the virus, the serum-virus mixture was incubated for 30 min at room temperature, and 50 µl of 1% turkey erythrocyte suspension was added in each well and mixed with the virus-antibody complexes using a benchtop shaker. The HI titer was determined as the reciprocal of end-point dilution showing complete inhibition of hemagglutination.

Transcriptional analysis

Transcriptional analysis of several innate immune genes was performed with trachea samples collected during necropsy at 5 dpi, as described in our previous studies (Ghorbani et al., 2020; Jang et al., 2018; Jang, Ngunjiri & Lee, 2016; KC et al., 2020). The uppermost 0.5 cm of tracheal tissue was excised using sterile scissors and forceps, placed into a tube with one mL TRIzol, snap-frozen in liquid nitrogen, and stored at −70 C until processing. Later, 0.6 mL dry volume mixed 0.5 mm and 0.1 mm ZR BashingBead lysis matrix (Zymo Research, Tustin, CA, USA) was decanted into tubes containing the tracheal samples. The tissues were then homogenized using a TissueLyser II (Qiagen Sciences Inc., Germantown, MD) at a frequency of 30 Hz for 5–7 min. Once the tissue was completely homogenized, chloroform was added to the TRIzol at a ratio of 1-part chloroform to 5-parts TRIzol. The mixture was vortexed vigorously for 15 s and centrifuged (16,000× g, 4 °C, 20 min). Crude RNA (the top, transparent layer) was transferred to a new tube, taking care not to disturb the white DNA layer below. An equal volume of EtOH was added to the tube and the mixture was vortexed before being transferred in entirety to a Zymo-Spin IIC Column (Zymo Research, Tustin, CA, USA) and centrifuged (10,000× g, 30 s). 400 µL RNA Wash Buffer were added to the column, which was centrifuged (16,000× g, 30 s), transferred to a new collection tube, and centrifuged again (16,000× g, 2 min) to ensure that the filter was dry. Further steps in the Direct-zol RNA MiniPrep Plus Quick Protocol (Zymo Research, Tustin, CA, USA) were followed, including the recommended DNase I treatment step. Centrifugation for Steps 3–4 was conducted at 16,000× g for 30 s. And additional centrifugation (to ensure filter dryness) was added after Step 4 (16,000× g, 2 min). Final elution volume was changed to 50 µL, incubation time was increased to 1 min (room temperature), and final centrifugation was conducted at 10,000× g for 1 min. The RNA was stored at −80 °C until use.

Complementary DNA (cDNA) was synthesized using GoScript reverse transcriptase (Promega, Madison, WI, USA) with oligo(dT)15 primer. RNA samples were thawed and quantified using NanoDrop™ 2000 (Thermo Fisher Scientific, Waltham, MA, USA) and diluted to a standard concentration of 250 ng per 10 µL aliquot before adding the 2 µL of Oligo(dT)₁₅ primer (pre-diluted to 250 ng/ µL). The RNA-primer mixture was incubated at 72 °C for 5 min and immediately placed on ice. Next, 8 µl of master mix containing 4 µl of GoScript Reaction Buffer (5X), 2 µl of MgCl2 (25 mM), 1 µl of deoxynucleoside triphosphate (dNTP) mix (containing 10 mM of each dNTP), and 1 µl of reverse transcriptase enzyme (160 U/ µL) was added to each tube and synthesis of the first cDNA strand was allowed to take place at 42 °C for 1 h. Then, quantitative PCR (qPCR) was performed using PerfeCTa SYBR green FastMix Low ROX kit (QuantaBio, Beverly, MA, USA) in a 7500 real-time PCR system (Applied Biosystems, Waltham, MA, USA). Briefly, cDNA was diluted 1.5625 times and 5 µl of the diluted cDNA was mixed with 15 µl of master mix containing 10 µl of SYBR Green FastMix (2X) and 2.5 µl of each primer (20 µM) in a 0.2-ml tube. The cycle conditions of quantitative PCR were 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. Afterward, melting curve analysis was performed from 60 °C to 90 °C. Primers for amplification of IFN-α, IFN-β, IFN-γ, interleukin-6 (IL-6), and lipopolysaccharide-induced tumor necrosis factor-α factor (LITAF) genes were published in an earlier study (Sharafeldin et al., 2015), while those for amplification of 2′–5′-oligoadenylate synthase (OAS) and myxovirus (influenza virus) resistance 1 (Mx), IFN-λ, and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) genes were designed in this study. All primer sequences, accession numbers for the target genes, and sizes of the amplified products are provided in Table 1. Differential gene expressions were normalized using GAPDH gene as the internal control. Relative mRNA expression levels were calculated as fold changes over the mock group using the 2−ΔΔCT method (Ghorbani et al., 2020; Jang et al., 2018).

Viral titration by quantitative real-time RT-PCR

Virus quantification was performed using supernatants from gut sample (cecum and ileum tissues and digesta contained therein) homogenates and respiratory tract (nasal and tracheal) washes. Procedures for respiratory washes and gut sample are detailed in our recent reports (Abundo et al., 2020; Ngunjiri et al., 2019b; Taylor et al., 2020). Gut sample supernatants were generated by adding 500 µL of sterile PBS to tubes containing 0.3 g of homogenates followed by vigorous vortexing for 2 min and centrifugation at 10,000× g at 4 °C for 5 min. Respiratory sample supernatants were generated by centrifugation of tracheal and sinus washes at 10,000× g at 4 °C for 10 min. Viral RNA (vRNA) was extracted from 200 µL of supernatant using RNeasy Mini Kit (Qiagen Sciences Inc., Germantown, MD, USA). The amount of viral genomes was determined using quantitative reverse transcription PCR (qRT-PCR) with type A influenza virus matrix gene specific primers and probes, and EID₅₀ viral titers were extrapolated from standard curves generated with the stocks used for infection as described in our previous reports (Ghorbani et al., 2019; Jang, Ngunjiri & Lee, 2016). To generate a standard curve, viral RNA was extracted from serial 10-fold dilutions of the same virus stock (with known EID₅₀ titer) used to inoculate the turkeys and subjected to qRT-PCR to generate cycle threshold (Ct) values. A scatter plot of Ct values (x-axis) against log₁₀-transformed virus dilution factors (y-axis) was made and fitted with a linear trendline to generate the equation y =bx + a, where y is the virus concentration (log₁₀ EID₅₀/ml), b is the slope of the trendline, x is the Ct value, and a is the y-intercept. This equation was then used to extrapolate gut and respiratory tract viral titers by substituting the Ct values obtained from sample-derived viral RNA for x.

Bursal atrophy

Bursal atrophy was determined based on bursa-body weight ratio (Cazaban et al., 2015). Body weight of individual birds were taken immediately after euthanasia before any tissues were excised. Full-size bursas were carefully excised during necropsy and weighed individually. The bursa-body weight ratio was calculated as follows: bursa-body weight ratio = [bursa weight (g)/body weight (g)] ×1000 (Stoute et al., 2013).

Sample collection and processing for 16S rRNA gene sequencing

Samples were collected from euthanized birds at 5 and 14 dpi. Collection, storage, and processing of respiratory and gut samples, including sinus and tracheal washes and gut (ileum and cecum) sample homogenization, were performed exactly as detailed in our recent reports (Abundo et al., 2020; Ngunjiri et al., 2019b; Taylor et al., 2020). For lower respiratory tract lavage collection, we assembled an apparatus by fitting a sterile non-filtered 200 µL pipet tip to a sterile 25 mL Luer Lock syringe. The bottom 3–5 mm of the tip was then trimmed with sterile scissors to broaden the channel. The syringe was filled with of sterile PBS before inserting the pipet tip into the lower trachea (about 0.5 cm above the syrinx) and securing it airtight with forceps before pumping the PBS into the lower respiratory tract. This process was repeated until about 40 mL of PBS filled the lower respiratory tract. After the last volume of PBS was added, the PBS was flushed in and out of the lower respiratory tract several (∼15–20) times until it became cloudy. As much of the PBS was retrieved as possible (typically yielding ∼30 mL) and transferred to a sterile tube. The tubes were kept on ice until processing. Later, the tubes were centrifuged at 7,830 ×g for 20 min. The supernatant was decanted, and pellets transferred to two mL for storage at −70 °C until DNA extraction.

DNA extraction for 16S rRNA gene analysis

DNA was extracted from respiratory wash pellets using the DNeasy Blood & Tissue kit (Qiagen Sciences Inc., Germantown, MD, USA) as described previously (Ngunjiri et al., 2019b) with some specific modifications, which were optimized based on sample type to improve DNA yield and minimize PCR inhibitors. In particular, enzymatic lysis buffer was used instead of Buffer ATL. Lysis buffer and reconstituted lysozyme were prepared beforehand and stored separately. Lysozyme was added to lysis buffer immediately before use at concentrations of 0.02 M for nasal wash samples, and 0.031M for tracheal and lower respiratory lavage samples. The disparity in reagent concentrations is due to the volumes of PBS added to samples to resuspend pellets.

Lysozyme was reconstituted by slowly adding approximately 1.75 mL sodium acetate (C₂H₃NaO₂) to 1.5 g lysozyme desiccate in a 15mL round-bottom tube and vortexing until thoroughly dissolved. Once the mixture was suspended, three mL of sterile glycerol was added and vortexed. The final volume was six mL, with a concentration of 250 mg/mL, was aliquoted to two mL screwcap tubes (250 µL/tube) and stored at −20 °C until use.

For nasal wash pellets, lysis buffer was prepared by combining Tris reagent (pH 8), EDTA, and Triton X to final concentrations of 0.02 M, 0.002 M, and 1.2%, respectively, in nuclease-free water. Lysis buffer for tracheal wash and lower respiratory lavage pellets was prepared in a similar manner except that the final concentrations were 0.031 M, 0.0031 M, and 1.867%, respectively. Lysis buffer was stored at room temperature until use.

Nasal wash pellets were resuspended in 180 uL of enzymatic lysis buffer (lysis buffer plus lysozyme at 20 mg/ml) by adding 2 sterile glass beads and vortexing to break up the pellet. Samples were incubated at 37 °C for 2 h. After adding 25 µL of proteinase K and 200uL Buffer AL, samples were incubated again at 56 °C for 30 min. A 200 µL volume of 100% ethanol was added to each sample, and the total sample volume was transferred to a filter spin column and centrifuged (10,000×g for 1 min). The column was washed by applying 730 µL of Buffers AW1 and AW2, with centrifuging between each wash (10,000×g, 1 min). To elute DNA from the column filters, 50 µL of Buffer AE was applied directly to the center of the column, then the samples were incubated at room temperature for 5 min and centrifuged (10,000×g, 1 min). The eluate was reapplied to the filter, incubated, centrifuged a second time.

The procedure for tracheal wash and lower respiratory lavage pellets differed from the above protocol for nasal wash samples in the following ways. First, 500 µL of sterile 1X PBS was used to resuspend the pellet, then 900 uL enzymatic lysis buffer was added. Second, one mL of Buffer AL was added after the first incubation step (37 °C for 45 min), and one mL 100% ethanol was added after the second incubation step (56 °C for 30 min). Running the total volume of solution through the spin column required 5-6 rounds of centrifugation. Finally, 500 µL of Buffers AW1 and AW2 were used to wash the filter.

The protocol for DNA extraction from ileum and cecum homogenate samples was modified from the DNeasy PowerSoil Kit and the extraction was performed exactly as detailed in our recent report (Ngunjiri et al., 2019b).

16S rRNA gene amplification, quantitation, and sequencing

Our optimized methods for sequencing and sequence processing are detailed in our recent reports (Abundo et al., 2021; Taylor et al., 2020). For each sample, DNA was quantified using NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), adjusted to a concentration of 100 ng/µL, and shipped to the University of Minnesota Genomics Center (Minneapolis, MN, USA) for sequencing. The number of 16S rRNA gene copies in each sample was determined through quantitative PCR (qPCR) of the V4 hypervariable region using primer sets 515F (5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGGACTACHVGGGTWTCTAAT-3′) Nextera primers (Caporaso et al., 2012). The conditions for qPCR were as follows: one cycle of 95 °C for 5 min; followed by 35 cycles of 98 °C for 20 s, 55 °C for 15 s, and 72 °C for 1 min; and one cycle of 72 °C for 5 min using QuantStudio Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA). The Cycle Threshold (CT) values were then converted to 16S rRNA gene copy numbers using a standard curve. The number of copies per ng of DNA was estimated by the ratio of total 16S rRNA copies to the DNA concentration/uL multiplied by the volume of DNA used in the copy quantification step (5 µL) and used as a proxy for bacteria content. Group comparisons for the bacterial content were performed using the Wilcoxon rank sum test with Holm correction (‘stats::pairwise.wilcox.test’).

Before sequencing, the samples were normalized to 1.67 × 10⁵ 16S rRNA gene copies/ µL then subjected to two rounds of PCR. In the first round, the 515F and 806R Nextera primers were used for the amplification of V4 region from 3 μL of normalized DNA using the following cycling parameters: one cycle of 95 °C for 5 min, followed by 25 cycles of 98 °C for 20 s, 55 °C for 15 s, and 72 °C for 1 min. The PCR product was then diluted to 1:100, and 5 µL was used in the second round of PCR using forward (5′-AATGATACGGCGACCACCGAGATCTACAC(i5)TCGTCGGCAGCGTC–3′) and reverse (5′-CAAGCAGAAGACGGCATACGAGAT(i7)GTCTCGTGGGCTCGG–3′) indexing primers (Integrated DNA Technologies, Coralville, IA). In the second PCR, the following cycling parameters were used: 1 cycle at 95 °C for 5 min, followed by 10 cycles of 98 °C for 20 s, 55 °C for 15 s, and 72 °C for 1 min. Pooled size-selected samples were denatured using NaOH, diluted to 8 pM in Illumina’s HT1 buffer, spiked with 20% PhiX, and heat denatured at 96 °C for 2 min immediately before loading. MiSeq 600 (2 × 300 bp) cycle v3 kit (Illumina, San Diego, CA, USA) was used to sequence the samples. This protocol is an optimization of the Earth Microbiome Project protocol (Gilbert & Meyer, 2012; Gohl et al., 2016).

The barcoded Illumina sequencing data files were demultiplexed to generate fastq files for each sample followed by merging of R1 and R2 reads for each sample, removal of the amplification and Kit adapter primer sequences, and removal of very short (<245 bp) and very long (>260 bp) sequences.

The merged sequences were imported into QIIME2 (qiime2-2019.1) (Bolyen et al., 2019) for denoising through removal of PCR chimeras, clustering of amplicon sequence variants (ASVs) with 100% identities, removal of non-16S rRNA gene ASVs, generation of a feature table with ASV occurrence frequencies, generation of a mid-rooted phylogenetic tree, and taxonomic classification using SILVA 16S rRNA gene database (release 132). The elaborate protocol including all software used is described in our recent report (Taylor et al., 2020).

16S rRNA gene sequence analysis

Using QIIME2 (qiime2-2019.1) (Bolyen et al., 2019), the feature table was rarefied to 5,000 reads per sample by calling the ‘core-metrics-phylogenetic’ command in the ‘diversity’ plugin. All microbiome analyses reflect the numbers of samples/group given in Table 2. All other analyses were performed in RStudio (R Core Team, 2019) using the indicated packages. Alpha diversity was quantified as the number of unique ASVs observed per sample (‘vegan::specnumber’) (Oksanen et al., 2013). Group comparison was performed using pairwise Wilcoxon rank sum tests with Holm correction (‘stats::pairwise.wilcox.test’). Beta diversity was quantified as the total abundance-weighted phylogenetic distance between the constituent ASVs of one sample and another (‘GUniFrac::GUniFrac’) (Chen, 2018). Group comparisons were performed using pairwise PERMANOVA with Holm correction (‘pairwiseAdonis::pairwise.adonis’) (Martinez Arbizu, 2017).

Predominant taxa were identified for each body site as those ASVs that were present in all birds of at least one treatment group at a mean relative abundance of >1%. In order to draw more general conclusions from these taxa, the definition was expanded to encompass all ASVs belonging to the same lowest taxonomic classification (assigned by SILVA) as the identified predominant taxa. The site, infection group, and number of ASVs for each identified taxon is described in Table S1.

Differential abundance analysis was performed on this taxonomic set by taking the ratio of mean relative abundance (percent) of a genus for each virus-treated group to the mean relative abundance of the Mock group. The second log (Log₂) of this value describes the fold change of mean abundance in the virus-infected groups from the Mock. Statistical comparisons were performed using the Wilcoxon rank-sum test on untransformed values. P-values were adjusted post-hoc using the Holm method (‘stats::p.adjust’).

The distribution of predominant taxa among body sites and treatments was profiled according to presence at a treatment-level consensus relative abundance of >0.5%. The distribution of the taxa thus defined was necessary to ensure that correlations were performed with sufficient non-zero data. Relative abundance data were used for Spearman correlation (‘stats::cor’) with local viral titers (in nasal cavity, trachea, cecum, and ileum) and with local gene expression (trachea only). Significance of the correlation coefficient was assigned from a pairwise t-test comparison built into ‘psych::corr.test’ (Revelle, 2019), adjusted with the Holm method.

Sequence data and software availability

Raw data files and metadata are publicly available in the NCBI BioProject database under accession number PRJNA644054. Raw data files are also available through the Sequence Read Archive (SRA) under accession numbers SAMN15446589– SAMN15447100. Pipe scripts for QIIME2 sequence processing and generation of the ASV table and R scripts for analysis are publicly available at https://github.com/kjmtaylor22/cwl-cats.

Differential viral replication parallels the levels of survival, immune-related gene expression, bursal atrophy, and serum hemagglutination inhibition antibodies

To enable differential assessment of quantitative and qualitative differences in microbiota between virus infected and mock infected poults, we chose two H5N2 viruses (CKPA, isolated in chicken, and TKMN, isolated in turkey) that replicate at different levels in turkeys (Pillai et al., 2010a). In accordance with our previous observation (Pillai et al., 2010a), five poults in the CKPA group developed severe clinical signs such as gasping and swollen heads and were euthanized at 4, 6, and 7 days-post-infection (dpi). These birds were not sampled because they were inadequate for statistical tests and the experimental design specifically called for quantitative analysis of microbiota at 5 and 14 dpi and tracheal gene expression at 5 dpi. Birds in the Mock and TKMN groups did not exhibit clinical signs throughout the duration of the experiment (14 days). At 5 dpi, poults in the CKPA group shed significantly higher viral titers in the upper respiratory tract (URT) (nasal cavity and trachea) than the TKMN group (Fig. 1A). Viral titers in the lower intestinal tract (LIT) (ileum and cecum) were significantly lower than URT and statistically indistinguishable between CKPA and TKMN groups. Viral clearance in the infected birds was confirmed with nasal wash samples collected from the CKPA group at 14 dpi (Fig. 1A). Accordingly, TKMN virus was assumed to be cleared since its titers tended to be equal or lower than CKPA virus. No virus was detected in gut or respiratory samples from the Mock group.

LPAI virus infection results in local and systemic immune-related gene expression in avian species (Bergervoet et al., 2019; Xing et al., 2008). Here, we focused on local gene responses in trachea due to its tendency to shed higher titers of each virus than the other sites (Fig. 1A). The immune-related genes encompassed three categories: antiviral cytokines (IFN-α, IFN-β, IFN-γ, IFN-λ) (Schneider, Chevillotte & Rice, 2014), interferon inducible genes (OAS and Mx), and proinflammatory cytokines (IL-6, IFN-γ, LITAF) (Zhang & An, 2007). Except for lipopolysaccharide-induced tumor necrosis factor-α factor (LITAF), the genes examined here may be induced by double-stranded RNA produced as a replication intermediate by RNA viruses, including the influenza virus. Relative to the expression in mock-infected birds at 5 dpi, there was significant upregulation of IFN-α, IFN-β, IFN-λ, OAS, Mx, and IL-6 by both viruses, and IFN-γ by CKPA virus (Fig. 1B). It is worth noting that the levels of IFN-γ and IFN-λ mRNA induced by CKPA virus were about eightfold higher than those induced by the TKMN virus (p < 0.01). Only minor non-significant changes in the expression of LITAF were observed in virus-infected poults (Fig. 1B). Data from infected (CKPA and TKMN) and non-infected (Mock) birds were analyzed to demonstrate the relationship between virus replication and host gene expression through a series of regression plots (Fig. S1), with viral titer (Log₂ EID₅₀/ml) and gene expression level (Log₂ fold change over Mock) as variables. For the CKPA group, the strength of correlation decreased in the following order: OAS (r² = 0.94) >IFN-λ (r² = 0.84) >IFN-γ (r² = 0.71) >IFN-β (r² = 0.53) >IL-6 (r² = 0.51) >Mx (r² = 0.27) >IFN-α (r² = 0.25) >LITAF (r² = 0.17). The order of the strength of correlation was different for TKMN virus: OAS (r² = 0.87) >IFN-β (r² = 0.43) = IFN-α (r² = 0.43) >IFN-λ (r² = 0.41) >IL-6 (r² = 0.37) >Mx (r² = 0.27) >IFN-γ (r² = 0.22) >LITAF (r² < 0.01). Noteworthy, the CKPA virus showed significantly higher serum HI antibody titers at 14 dpi compared to TKMN virus (Fig. 2A), while the TKMN group displayed significant bursa atrophy at 5 dpi compared to the Mock group, although a complete recovery was observed at 14 dpi (Fig. 2B). Additionally, LPAI virus infection resulted in significant suppression of body weight at 5 dpi, with TKMN and CKPA being statistically indistinguishable (Fig. 2C). However, no significant body weight suppression was observed at 14 dpi (Fig. 2C).

Commensal microbiota is disrupted during active virus replication

Although commensal microbiota may demonstrate inhibitory effects against virus pathogens including influenza virus (Yitbarek et al., 2018a; Yitbarek et al., 2018b), the microbiota itself can be disrupted by virus-induced innate immune responses and pathologic effects. Here, we examined changes in bacterial microbiota during active virus replication (5 dpi), focusing on changes to total bacterial content (measured by 16S rRNA gene copies per nanogram of DNA per sample), diversity, composition, and abundance. The virus infected groups displayed significant increase in total bacterial content of the nasal cavity and ileum when compared to the Mock group, but not to each other (Fig. 3A). The increase in total bacterial content in the nasal cavity and ileum tended to be inversely associated with bacterial species richness (ASVs observed in each sample) (Fig. 3B). In the nasal cavity, although species richness was significantly decreased by both viruses compared to mock-infected birds, the average number of ASVs observed in the CKPA group was significantly lower than the TKMN group. In the ileum, only the TKMN group displayed a significant reduction of species richness compared to either the Mock or CKPA groups.

Dissimilarities in bacterial composition amongst birds in different treatment groups (beta diversity) were measured using abundance-weighted UniFrac distances. Following the trends observed in total bacterial content and alpha diversity (Figs. 3A and 3B), the ileum, cecum, and nasal cavity displayed significant virus-induced changes in beta diversity. The compositions of nasal microbiota in all three treatment groups were significantly different from one another (Fig. 3C, Table S2). For ileal and cecal microbiota, both virus-infected groups were significantly different from the Mock group, but not from each other.

To determine the impacts of virus infection on taxonomic composition of respiratory and gut microbiota, we targeted the taxa that were present in all birds of at least one treatment group at a mean relative abundance of >1%. These taxa were selected for each group and body site independently, and then summarized to genus level. Those taxa not classified to genus level were collapsed to family or order level while those identified as contaminants were dropped from further analysis, resulting in 23 taxa (Table S1). Contaminants were conservatively identified as ASVs present in all negative extraction control samples using the same extraction kit.

Log₂-transformed values of relative abundances were used to calculate differential abundance of the predominant bacteria between the Mock group and each of the virus-infected groups. As revealed by the heatmap in Fig. 4A, the abundances of several taxa were disrupted (increased or decreased) during active virus replication (5 dpi). For the CKPA group, 16 taxa were significantly disrupted in the nasal cavity, five in trachea, 1 in lower respiratory tract, one in ileum, and one in cecum. For the TKMN group, 10 taxa were significantly disrupted in the nasal cavity. Remarkably, the relative abundances of the 23 taxa were generally not significantly different between the infection groups (Fig. S2). Interestingly, bacteria from class Gammaproteobacteria were significantly disrupted across gut and respiratory sites; they tended to increase during active virus infection at 5 dpi (Fig. 4 and Fig. S3). Overall, the relative abundances of predominant genera and sample compositional variance explained by site were more similar between the three respiratory sites in CKPA group (more severe) compared to TKMN (less severe) and Mock (uninfected) groups (Fig. S4).

Disruption of microbiota persists after virus clearance

After virus clearance at 14 dpi (Fig. 1A), the nasal bacterial content was restored to similar levels as the Mock group, but the ileal bacterial contents of the TKMN virus-infected birds remained significantly higher than the Mock group (Fig. 5A). Notably, in contrast to the earlier observation at 5 dpi (Fig. 3A), the virus infected groups showed a significant increase of bacterial content in the trachea (Fig. 5A).

Virus-induced species richness depletion was still evident in nasal cavities at 14 dpi. However, there was a dramatic increase in the number of observed ASVs in ilea of TKMN virus-infected birds (Fig. 5B), despite diminishment of virus-induced disruption of bacterial beta diversity to non-significant levels in the guts of both infected groups. Bacterial community compositions of TKMN virus-infected birds remained disrupted in nasal cavities and new disruptions were observed in lower respiratory tracts (Fig. 5C). Even though the relative abundances of virtually all predominant genera in virus-infected birds remained disrupted at 14 dpi, only 3 taxa were significantly disrupted by the TKMN virus (Fig. 4B).

Microbiota disruption mainly correlates with local virus replication but not body weight

Taken together, the data presented above clearly show that influenza virus is responsible for the observed dysbiosis in all five body sites (Figs. 3, 4, 5, and Figs. S1–S3). Therefore, we performed a series of Spearman’s correlation tests to quantify the relationships between the abundance of site-associated predominant taxa and local viral titers (at 5 dpi) for the sites for which viral titers were measurable (excluding the lower respiratory tract; see Materials & Methods). Each infection group was compared with the Mock. The relative abundances of several taxa significantly correlated with viral titers depending on the body site.

In the nasal cavities, the abundances of unclassified Eggerthellaceae, Flavobacterium, unclassified Lachnospiraceae, Lactobacillus, Macrococcus, Romboutsia, Ruminococcus, Subdoligranulum, and unclassified Rhizobiaceae were negatively correlated with the titers of both TKMN and CKPA viruses, while those of Clostridium, Enterococcus, Escherichia-Shigella, and Pseudomonas were positively correlated (Fig. 6A). Further, Butyricicoccus and Staphylococcus abundances were negatively correlated with CKPA viral titers in the nasal cavities. Correlations observed in tracheas and ceca of TKMN virus-infected birds were generally weak and non-significant. However, in tracheas of CKPA infected birds, the abundances of Enterococcus, unclassified Erysipelotrichaceae, Escherichia-Shigella, Staphylococcus, and Pseudomonas were positively correlated with viral titers while Lactobacillus was negatively correlated. The abundances of Anaerostipes, unclassified Lachnospiraceae, Merdibacter, Sellimonas, and Subdoligranulum were negatively correlated with TKMN virus titers in the ileum. However, only Escherichia-Shigella (positive) and Lactobacillus (negative) were significantly correlated with ileal CKPA virus titers (Fig. 6A). Interestingly, despite the significant suppression of body weights in virus infected birds at 5 dpi (Fig. 2C), correlations between body weight and viral titers were generally weak and not statistically significant except for three correlations in nasal cavities, ceca, and ilea of TKMN virus-infected birds (Fig. 6B).

Tracheal microbiota disruption correlates with innate immune gene expression

Virus replication in trachea resulted in significant upregulation of IFN (α, β, γ, λ), IFN-stimulated (OAS and Mx), and IL-6 (proinflammatory) genes relative to mock-infected birds (Fig. 1B). Therefore, we paired the tracheal gene expression values (Log₂ fold-change) with relative abundance values for predorminant bacterial taxa in the same birds euthanized at 5 dpi and performed Spearman’s correlation tests as described above. None of correlations in the TKMN group was statistically significant (Fig. 7A). However, the abundances of 5 taxa in the CKPA group were statistically correlated with the expression levels of one or more immune-related genes (Fig. 7B). As expected, because the LITAF gene expression is not stimulated by influenza virus infection (Fig. 1B and Fig. S1), its correlations with microbiota abundances were very weak and not significant (Fig. 7B). Of the IFN-related genes examined in the CKPA group, IFN-α and Mx were poorly correlated with dysbiosis, again reflecting the low levels detected at 5 dpi. Upregulation of IFN-γ, IFN-λ and OAS genes was correlated with significant increase in the abundance of Clostridium, Pseudomonas, and Staphylococcus, and a significant decrease in Lactobacillus. Additionally, OAS gene upregulation was correlated with significant increase in the abundance of Escherichia-Shigella, while the expression of IL-6 proinflammatory gene correlated with significant increase of Staphylococcus and depletion of Lactobacillus (Fig. 7B).