Assessment of two DNA extraction kits for profiling poultry respiratory microbiota from multiple sample types

Characterization of poultry microbiota is becoming increasingly important due to the growing need for microbiome-based interventions to improve poultry health and production performance. However, the lack of standardized protocols for sampling, sample processing, DNA extraction, sequencing, and bioinformatic analysis can hinder data comparison between studies. Here, we investigated how the DNA extraction process affects microbial community compositions and diversity metrics in different chicken respiratory sample types including choanal and tracheal swabs, nasal cavity and tracheal washes, and lower respiratory lavage. We did a side-by-side comparison of the performances of Qiagen DNeasy blood and tissue (BT) and ZymoBIOMICS DNA Miniprep (ZB) kits. In general, samples extracted with the BT kit yielded higher concentrations of total DNA while those extracted with the ZB kit contained higher numbers of bacterial 16S rRNA gene copies per unit volume. Therefore, the samples were normalized to equal amounts of 16S rRNA gene copies prior to sequencing. For each sample type, all predominant bacterial taxa detected in samples extracted with one kit were present in replicate samples extracted with the other kit and did not show significant differences at the class level. However, a few differentially abundant shared taxa were observed at family and genus levels. Furthermore, between-kit differences in alpha and beta diversity metrics at the amplicon sequence variant level were statistically indistinguishable. Therefore, both kits perform similarly in terms of 16S rRNA gene-based poultry microbiome analysis for the sample types analyzed in this study.


Introduction
The importance of gut and respiratory microbiota in poultry health and production performance has resulted in an increase of studies aimed towards developing microbiome-based intervention  [17]. The chickens were provided with ad libitum access to feed and water and their welfare was monitored twice daily. All birds were clinically healthy and were all included in the study. Prior to the collection of invasive samples, the birds were humanely euthanized through carbon dioxide exposure as previously described [17].

Experimental design and sample collection
A total of 72 four-week-old specific-pathogen-free chickens were hatched in the same incubator, then co-housed together for four weeks prior to sampling. At the day of sampling, the birds were randomly assigned to three groups for sample collection based on the target sample types and sampling techniques (Fig 1). The sample size was decided based on our recent reports showing that four or more respiratory samples were sufficient for 16S rRNA genebased taxonomic and diversity analysis [17]. An assortment of invasive and non-invasive respiratory samples including, nasal cavity wash, upper and lower tracheal wash, lower respiratory lavage, as well as choanal and tracheal swabs, were obtained from either live or euthanized birds. For the collection of tracheal and choanal swab samples, the swabs were gently inserted into the trachea or choanal cleft, moved 5 times along the cleft or the trachea, and then subsequently immersed in 1 mL of 1X PBS. To collect tracheal washes, the upper and lower tracheal tissues were aseptically excised, stored in a sterile tube, then immediately snap-frozen in liquid nitrogen. Each frozen trachea was thawed, and the luminal mucosa was washed using 1X PBS in a biosafety cabinet. For the collection of nasal cavity wash, the outer nares were thoroughly wiped with a cotton ball containing 100% ethanol to remove any visible organic material. The nasal cavity was then flushed through the nostrils with 1X PBS, until a total of approximately 1 mL of nasal wash was collected from each bird. For the collection of lower respiratory lavage, a 25 mL syringe containing 1X PBS with a fitted pipet tip was used to wash the lower respiratory  The samples were pooled as follows: 24  samples per pool for tracheal swabs from live birds, choanal swabs from live birds, tracheal  swabs from euthanized birds, choanal swabs from euthanized birds, upper tracheal wash pellets, lower tracheal wash pellets; and 48 samples per pool for nasal cavity wash pellets and lower respiratory lavage pellets. Pooling of the lower respiratory lavage and nasal cavity wash samples depended on whether the trachea was to be excised, or not. Birds in group 3 could not be used for lower respiratory lavage for practical reasons, since the entire trachea was excised leaving no place to attach the lavage collection apparatus. On the other hand, in groups 1 and 2, the trachea was cut at about 5 mm above the syrinx to support the attachment of the apparatus. Likewise, the nasal cavity wash was collected only from swabbed birds after swabbing, to prevent cross contamination between the nasal cavity and the trachea, since the procedure of washing the nasal cavity can have some leakage to the trachea. The pooled samples were homogenized and divided into 8 replicates for tracheal and choanal swabs and tracheal wash pellets, and into 16 replicates for the nasal cavity wash pellets and lower respiratory lavage pellets. The sample pools were evenly divided to be extracted with either the BT or ZB kit. This study did not involve treatment groups and, therefore, no control animals were required.

Sample processing and DNA extraction
All the respiratory samples were pelleted by centrifugation at 10,000 × g at 4˚C for 10 mins.
After the supernatant was removed, the pellet was resuspended and divided accordingly into replicates (Fig 1). The general steps for BT kit extraction were as follows: all the respiratory samples were pretreated for gram-positive bacteria using an enzymatic lysis buffer containing 20 mM Tris�Cl, pH 8.0, 2 mM sodium EDTA, 1.2% Triton X-100, and 20 mg/mL of lysozyme. After incubation time (3 hours for nasal cavity wash and swab samples, 45 min for the trachea wash and lower respiratory lavage samples) at 37˚C, 100% ethanol was added. The mixture was then added to a spin column for DNA isolation, and each column was washed with buffers AW1 and AW2, respectively. Finally, 50 μL of elution buffer AE was used to elute DNA from the spin column membrane.
Since our samples were collected from multiple respiratory sites, subtle modifications to the BT kit extraction protocol were implemented. For extraction of nasal cavity wash, tracheal swab, and choanal swab samples, the manufacturer's protocol was modified to include the lysozyme-based enzymatic lysis buffer option for disruption of gram-positive bacteria and increasing the lysis incubation time from 30 min to 3 hours (MBT-A protocol) [29]. The longer incubation time was implemented due to the high levels of organic material in these samples [17]. Upper tracheal wash, lower tracheal wash, and lower respiratory lavage samples were extracted using another modified protocol where the enzymatic lysis incubation time was increased from 30 min to 45 min (MBT-B protocol) as previously described [17]. In the MBT-A protocol, proteinase K, Buffer AL, and 100% ethanol volumes were as recommended by the manufacturer's protocol, while AW1 and AW2 Wash Buffer volumes were increased from 500 μL to 730 μL. The volumes of these reagents were increased 5 times in the MBT-B protocol to improve DNA quantities and purity as previously described [17]. In both protocols, DNA was eluted using 50 μL of elution Buffer AE, reapplied in the column, incubated for 5 min, and eluted once more.
Extraction of DNA from respiratory samples using the ZB kit was performed as per the manufacturer's direction [30]. Identical ZB kit extraction protocol was implemented regardless of the sample type. The steps for the ZB kit extraction were as follows: initially, 730 μL of lysis solution and beads (size = 0.1 & 0.5 mm) were added to the respiratory pellets. The tubes were fastened to a Disruptor Genie vortex (Scientific Industries, Bohemia, NY) using a vortex adapter. The samples were vortexed at maximum speed for 20 min. The samples were then centrifuged at 8,000 × g for 1 min, and 400 μL of supernatant was added to III-F filter and centrifuged once more at 8,000 × g for 1 min. Afterwards 1,200 μL of DNA binding buffer was added to the filtrate and the mixture was transferred to IICR Column and centrifuged at 10,000 × g for 1 min. Once the flow-through was discarded, the column was washed with 400 μL of Wash Buffer 1 and 700 μL of Wash Buffer 2, centrifuged at 10,000 × g, and the flowthrough was discarded. Next, an additional 200 μL of Wash Buffer 2 was added to the column and centrifuged for 10,000 × g for 1 min. For the elution process, 50 μL of DNase/RNase free water was directly added to the column, incubated for 1 min, and centrifuged at 16,000 × g for 3 min. As a final step, the DNA eluate was passed through the III-HRC prepped filter and centrifuged at 16,000 × g for 3 min.
Both a mock community control (ZymoBIOMICS Microbial community standard) and negative extraction controls (1X PBS) were extracted alongside the samples in both kits. At least two negative extraction controls were collected with each sample type for each kit. The negative extraction controls were subjected to each step, all the way from sample collection to 16S rRNA gene sequencing.

16S rRNA gene sequencing
Our optimized methods for pre-sequencing and sequence processing are detailed in our recent reports [15][16][17]. Briefly, DNA concentration was quantified using NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA) and adjusted to a concentration of 100 ng/μL. Next, 16S rRNA genes in the DNA samples were PCR-amplified and sequenced at the University of Minnesota Genomics Center. The 16S rRNA gene copies in the samples were determined by using quantitative PCR (qPCR) of the V4 hypervariable region using primer sets 515F (5 0 -TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGTGCCAGCMGCCGCGGT AA-3 0 ) and 806R (5 0 -GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGGACTACHV GGGTWTCTAAT-3 0 ) Nextera primers [31]. 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). The Cycle Threshold (CT) values were then converted to 16S rRNA gene copy numbers using a standard curve. The samples were then normalized to 1.67 × 10 5 16S rRNA gene copies/μL then subjected to two rounds of PCR. In the first round, 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 0 -AATGATACGGCGA CCACCGAGATCTACAC(i5)TCGTCGGCAGCGTC-3 0 ) and reverse (5 0 -CAAGCAGAAGACGG CATACGAGAT(i7)GTCTCGTGGGCTCGG-3 0 ) 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) was used to sequence the samples. This protocol is an optimization of the Earth Microbiome Project protocol [32,33].

Data processing and statistical analysis
Fastq files for each sample were generated after sequence demultiplexing. BBTools (v35), a suite of DNA/RNA analysis software provided by the Joint Genome Institute, was used to trim both proximal and distal primers and adapters from the sequences using the BBDuk module, and paired-end reads were merged using the BBMerge module [34]. Cutadapt (v1.4.2) was used to remove 16S rRNA gene amplicon primers [35]. The BBMap module of BBTools were used to filter sequence lengths of less than 245 and greater than 260 base pairs [34]. QIIME 2 (qiime2-2019.1) and software implemented therein, were used for additional data processing and statistical analysis [36]. Using a Casava 1.8 single-end demultiplexed format, sequences were imported into QIIME 2. DADA2, pooled chimera filtering method, was used to denoise the sequences [37]. VSEARCH [38] was used for chimeric sequence detection, non-16S rRNA gene identification, and open reference clustering of amplicon sequence variants (ASVs) with a similarity threshold of 100% using the SILVA (release 132) reference database [39]. Furthermore, ASVs with a total frequency of < 10 in the entire table and present in only 1 sample were removed. In addition, to generate a feature table using MAFFT, a de novo alignment of the representative set of sequences was used, which eliminated non-conserved and highly gapped alignment columns by applying the 'mast' option of the alignment plugin. The aligned filtered sequences were used to reconstruct a phylogenetic tree using Fast Tree, that was then rooted at the mid-point prior to beta diversity analysis.
A Naïve Bayes classifier was trained using the 16S rRNA gene sequences covered in the 515F (GTGYCAGCMGCCGCGGTAA) and 806R (GGACTACNVGGGTWTCTAAT) primer pair region for taxonomic classification [40]. For the respiratory samples, the classifier was trained on the reference sequences derived from the SILVA 16S rRNA gene database (release 132) based on matches with this primer pair and a taxonomic identification table (majority classification) based on 99% identity clustered reference sequences. For the mock community samples, the classifier was trained using the 16S rRNA sequences provided by ZymoBIOMICS with identical parameters as described above. Finally, the classifier was used in our experimental data to assign taxonomy to the observed ASVs. 5000 reads per sample was the rarefaction cut-off in the final feature table [15][16][17].
Weighted and unweighted UniFrac distance metrics were used to assess beta diversity variations in the bacterial community structure between extraction kits [41]. Statistical differences in the beta diversity for each sample type between the kits were quantified by pairwise permutational multivariate analysis of variance (PERMANOVA) with a Benjamini-Hochberg correction for false discovery [42]. The number of observed ASVs [43] was used for alpha diversity analysis to evaluate the richness of bacterial species for each sample, and Pielou's evenness [44] was used to evaluate the relative evenness of species richness. Statistical differences between DNA quality and quantity, 16S rRNA gene copies in total DNA per sample, number of 16S rRNA gene sequences per sample, avian mitochondria sequences per sample, and alpha diversity were examined using a nonparametric Kruskal-Wallis test with Benjamini-Hochberg correction for false discovery [45]. Differences were considered significant at p < 0.05.
Additionally, differentially abundant analysis (DESeq2) and contamination removal (decontam) were conducted in R [46]. DESeq2 was used to detect statistically significant differentially abundant taxa between kits for each respiratory sample type [47]. Adjusted p values < 0.05 (Wald statistic, with a Benjamini-Hochberg correction for false discovery) were considered to be significant [47]. Decontam frequency method (version 1.2.1) was used to identify taxa with significant inverse correlations with 16S rRNA gene concentration as measured via qPCR [48]. Moreover, respiratory samples were filtered based on the number of 16S rRNA gene copies found in the negative extraction control. Samples that contained lower 16S rRNA gene copies compared to the respective negative extraction controls were removed from analysis.

BT and ZB kits differ greatly in DNA yield and 16S rRNA gene copy numbers
For both kits, DNA was eluted in a volume of 50 μL to allow between-kits comparison of DNA concentrations and purity. DNA quality was spectrophotometrically measured and expressed as ratios of absorbances at 260 nm and 280 nm wavelengths (A260/A280) and 260 nm and 230 nm wavelengths (A260/A230). The mean A260/A280 values for most of the samples were around 1.8 and were not statistically different between the kits. However, the mean value for the lower respiratory lavage samples extracted with the ZB kit was higher than 2.0 which was significantly higher compared to the BT kit (Fig 2A). On the other hand, higher mean A260/A230 ratios were observed for most of the samples extracted using the BT kit, with significantly higher ratios in the live choanal swab, euthanized choanal swab, euthanized tracheal swab, nasal cavity wash, and lower respiratory lavage (Fig 2B). DNA yield was spectrophotometrically measured and expressed as mass per unit volume (ng/μL). In general, the BT kit tended to have more DNA yield than the ZB kit, with DNA quantities being significantly higher for most of the sample types except for the lower trachea wash and lower respiratory lavage (Fig 2C).
The DNA was subsequently subjected to qPCR to determine the number of 16S rRNA gene copies. Negative extraction controls were included to facilitate filtering out of samples with low amounts of 16S rRNA gene copies. The highest amounts of 16S rRNA gene copies detected in the negative controls were set as elimination thresholds for the corresponding sample types. The elimination cut-offs for the BT kit were 1.6 × 10 3 , 1.4 × 10 4 , 2.0 × 10 3 , 3.1 × 10 3, and 1.1 × 10 4 16S rRNA gene copies per nanogram of total DNA for choanal swab, lower respiratory lavage, nasal cavity wash, tracheal swab, and lower and upper tracheal washes, respectively. The elimination cut-offs for the ZB kit were 20, 50, 20, 60, and 15 16S rRNA gene copies per nanogram of DNA for choanal swab, lower respiratory lavage, nasal cavity wash, tracheal swab, and tracheal wash, respectively. Using these criteria, 3 out of the 40 samples extracted using the BT kit were dropped from further analysis (1 live tracheal swab and 2 lower respiratory lavage samples). In contrast, all 40 samples extracted using the ZB kit were retained.
Among the retained samples, the nasal cavity wash, euthanized tracheal swab, and euthanized choanal swab samples extracted with the ZB kit had significantly higher numbers of 16S rRNA gene copies than those extracted using the BT kit ( Fig 3A). Subsequently, all samples were normalized to 1.67 × 10 5 16S rRNA gene copies/μL. This was done to achieve 5 × 10 5 16S rRNA gene copies in 3μL of normalized sample used for pre-sequencing PCR, which approximately translates to a target sequencing coverage of at least 10X. The PCR products were subsequently used for sequencing library preparation. No significant differences in the 16S rRNA gene sequence counts from the normalized samples were observed between the kits, except for the lower respiratory lavage samples (Fig 3B). Additionally, a substantial number of non-16S rRNA gene sequences, including host avian mitochondria gene sequences, were detected in samples extracted with either kit (Fig 4). Avian mitochondria gene sequences were 175, 176, 178, and 194 bp in length, while the 16S RNA gene sequence lengths were between 250-264 bp (Fig 4A). For all sample types, no significant differences were detected in the number of avian mitochondria gene sequences between the kits (Fig 4B).

Minimal contaminants were identified, and the kits displayed subtle differences in extraction efficiencies
Several bacterial taxa were found in the negative extraction controls, indicating the presence of contaminants originating from either the reagents, sample processing steps, or sequence library preparation. These taxa were analyzed using decontam R package to identify contaminants in the respiratory samples [48]. Of the 1487 ASVs detected among the samples, 14 were identified to be contaminants ( Table 1). The contaminants were detected in low abundances (< 0.15% relative abundance) independently of the kit used for extraction.
A mock community control was processed alongside the respiratory samples to assess the efficiency of extraction between the kits and the extent of contamination during the extraction process. Contaminants in the mock community samples were identified based on their absence in the compositional shotgun sequence data provided by ZymoBIOMICS in the certificate of analysis [49]. Consistent with our previous findings [17], the mock community had < 1% contaminant bacteria (Table 1). While all bacteria in the mock community were

PLOS ONE
detected in samples extracted with either kit, their relative abundances based on 16S rRNA gene sequencing varied compared to the theoretical 16S rRNA gene abundances (Table 2) [50]. Even though disparities in the relative abundance of 16S rRNA gene sequences were detected between the kits, only Listeria monocytogenes and Staphylococcus aureus, both grampositive bacteria, were found to be significantly higher in the ZB kit as detected using DESeq 2 [47]. However, this difference may be due to the low number of mock community replicates available for this study ( Table 2).

Method of extraction had minimal effects on taxonomic composition and diversity
After filtering the sequence data based on the 16S rRNA gene copy numbers and contaminants in negative extraction controls, the impacts of the extraction kits on basic microbiota parameters were assessed. Both of the extraction kits were able to detect the following predominant taxa at the class level: Bacilli, Clostridia, Gammaproteobacteria, Bacteroida, Erysipelotrichia, Actinobacteria, Mollicutes, and Alphaproteobacteria (Fig 5). Furthermore, several ASVs that were shared between the extraction methods within each sample type constituted > 92% of total abundance in the corresponding bacterial community (Table 3). Although several unique ASVs were detected by either kit, their relative abundances were low (S2 File). In addition,  (Table 4). While majority of these differentially abundant taxa had very low relative abundance, Vibrio prevalence was highly pronounced in the BT kit.
The characteristics of the microbial communities detected by each kit were evaluated further through alpha and beta diversity metrics. In majority of the samples, the number of observed ASVs in a given sample type were statistically indistinguishable between the kits, except in the live choanal swab where, on average, samples extracted using the BT kit had a higher number of ASVs (Fig 6A). Similarly, the extraction kit did not affect the species evenness in any of the sample types (Fig 6B). Beta diversity was assessed by comparing within sample type UniFrac distances between kits. No statistical differences were observed between the kits (Fig 7).

Discussion
The DNA extraction process can be a major source of variation between microbiome studies [7]. This poses a potential problem with inter-study comparison among the growing number of poultry respiratory microbiome studies. Using the BT kit in our previous studies, we were able to characterize the microbial communities present in sinus cavities, tracheas and lower respiratory tracts of turkeys, chicken broilers, and chicken layers raised in commercial settings, and specific-pathogen free chickens raised for research [2,[15][16][17]. Although we optimized the BT kit for increased DNA yields and removal of PCR inhibitors from respiratory samples [2,15,17], the kit is not specifically certified for microbiome sample extraction [29]. Here, we compared the performance of the BT kit with that of the ZB kit, which is validated for extraction of low-bioburden microbiome samples [30]. Overall, while statistical differences were observed between kits in the pre-sequencing results, post-sequencing microbiota data obtained with normalized samples showed only minimal differences between the kits.
Both kits were very efficient in eliminating protein contaminants absorbing at 280 nm as indicated by A260/A280 ratios within the expected range of 1.8-2.0 [29,51,52]. However, the ZB kit seemed inefficient at removing non-protein contaminants absorbing at 230 nm from some sample types, especially the nasal cavity wash and lower respiratory lavage samples, as indicated by A260/A230 ratios that were below the expected range of 2.2-2.5 (Fig 2B) [53]. While the exact nature of these impurities remain undetermined, their impact on PCR amplification was negligible as samples extracted with the ZB kit tended to have similar or higher 16S rRNA gene copy numbers compared to those extracted with the BT kit ( Fig 3A). Nevertheless, the ZB kit may require further optimization since the A260/A230 ratios reported in other studies comparing the performances of several commercial kits for bacterial DNA extraction from different sample types were all within the expected range [22,54,55]. On the other hand, the extraction process with the BT kit tended to acquire higher levels of contaminating bacteria resulting in 16S RNA gene copy cut-off values that were 2 to 3 orders of magnitude higher compared to the ZB kit. We speculate that the minimal steps and low bioburden reagents in the ZB kit contributed to lower 16S copies found in the negative extraction controls. Our protocol for 16S rRNA sequencing on the Illumina platform calls for sample normalization to 5 × 10 5 16S rRNA copies/ 3μL before PCR amplification and library preparation to have a 10X target sequencing coverage. This normalization appeared to eliminate between-kit differences in the number of 16S rRNA gene sequences (Fig 3B). Studies evaluating kit effects on bacterial diversity typically focus on the differences observed in the concentration of 16S rRNA gene copies detected via quantitative PCR without paying much attention to the actual number of 16S rRNA gene sequences recovered [54,56]. The contribution of unnormalized DNA concentrations to the observed difference in microbiota profile between kits needs to be clarified. Furthermore, the number of non-16S rRNA gene sequences, most of which were derived from the host (Fig 4A), accounted for 18% and 24% of total sequences obtained through BT and ZB kit extraction processes, respectively, confirming that the 16S rRNA gene primer set used in this study was highly specific and PCR thermocycling conditions were optimal.
In addition to pre-sequencing data, bacterial compositions and diversity metrics are commonly used as proxies for performance in extraction kit-comparison studies [22,55,57,58]. Here, all the predominant ASVs detected in samples extracted with one kit were present in replicate samples extracted with the other kit and did not show significant differences at the genus level (Tables 3 and 4 and Fig 5). Probing further into shared taxa at higher levels of taxonomic resolution, a total of 11 genus level and 2 family level taxa were differentially abundant in one or more sample types (Table 4). Although majority of these differentially abundant taxa were low in relative abundance, Vibrio was more prevalent in samples extracted with the BT kit, which is intriguing since gram-negative bacteria are easy to lyse [59], and should therefore be represented at similar levels in samples extracted with either kit. However, between-kits   differences in alpha and beta diversity metrics were statistically indistinguishable (Figs 6 and  7). The high variability in the weighted UniFrac distances reflects differential abundances of several taxa between the kits (Table 4, S1 File).
Other studies have typically compared multiple extraction kits using a single sample type [22,55,57,58]. In this study, although post-sequencing analysis revealed comparable performances between the BT and ZB kits across multiple types of respiratory samples, the impact of sample type on other DNA extraction methods requires further investigation. For example, between-kit differences appear to be magnified in gut samples [54,60] and reduced in respiratory samples [56,57], possibly due to body site differences in microbial densities and other factors such as PCR inhibitors, complexity of matrix, high bacteria density and inhibitors present in the sample. Nevertheless, both BT and ZB kits performed similarly in terms of extracting specific bacteria in the mock community and showed differences in bacterial compositions and diversity metrics between different respiratory sites (Tables 2, 3), and have therefore validated our previous reports [2,[15][16][17]. While high numbers of unique ASVs were present in samples extracted with either kit, their total relative abundances were very low. The presence of unique ASVs in low biomass samples, such as respiratory samples, is not unexpected [23,24]. However, since there is no standard reference for avian respiratory microbiome, we cannot determine whether the unique ASVs are contaminants despite multiple steps of filtering, or rare taxa inhabiting the avian respiratory tract [61][62][63]. We suggest that these low abundance taxa be investigated further using a combination of culture dependent and independent approaches. With that said, the ZB kit is less likely to contain contaminant bacterial DNA, thereby reducing the likelihood of excluding samples with low DNA concentrations from analysis. Note that while both kits performed similarly in the detection of bacterial taxa from

PLOS ONE
different sample types, the number of ASVs from live choanal swab were exaggerated with the BT kit (Fig 6A), suggesting that this kit is not optimal for this sample type.
The design of this study resulted in uneven distribution of sample replicates between sample types. For each kit, there were 8 replicates/sample type for lavage and nasal wash samples and 4 replicates/sample type for tracheal wash and swab samples. Therefore, reasonable comparisons between sample types within each kit could not be achieved in this study.
In conclusion, both the ZB and BT methods of DNA extraction produced similar chicken respiratory microbiome data in terms of bacterial diversity and taxonomic profiles.
Supporting information S1 File. Data workbook. DESeq2 analysis for shared taxa within a sample type between kits. (XLSX) S2 File. Data workbook. Average relative abundance of shared and unshared taxa within a sample type between kits. (XLSX)