Environmental metagenomics enhances detection of circulating viruses from live poultry markets in Cambodia

Abstract

Environmental surveillance has emerged as a pivotal strategy for early detection of pathogens that pose a threat to humans. In Asia, live-bird markets (LBMs) are key human-animal interfaces for zoonotic virus transmission. Traditional sampling strategies are time-consuming, expensive and carry significant biosafety risks. Here, we assess the performance of metagenomics on environmental samples (ES) versus traditional poultry swabs for detecting viral pathogens in two Cambodian LBMs between January 2022 and April 2023. ES, including air (n = 35), cage swabs (n = 17), carcass wash water (n = 17) and drinking water (n = 9) are collected alongside oropharyngeal and cloacal swabs from chickens (n = 30) and ducks (n = 29). ES is sensitive in detecting 40 viruses from pathogen families including Orthomyxoviridae and Coronaviridae. Air samples capture the greatest diversity of poultry viruses. Viral contigs from ES show high sequence identity to poultry swab contigs when aligned to the same gene. We show ES outperforms poultry samples in detecting the highly pathogenic influenza A/H5N1, including clades 2.3.4.4b and 2.3.2.1c, which are found in the environment but are missed by poultry swabs. Our findings show metagenomics on ES replicates traditional surveillance, offering broader coverage and improved pathogen detection. This approach could be pivotal for mitigating zoonotic spillover and enhancing pandemic preparedness.

Fig. 1. Air detects more virus species than poultry swabs.

A Heatmap showing 84 different virus species which were detected at least once in poultry swabs over the course of the study. The red colour shows when a virus was identified in at least one poultry swab and at the exact same timepoint was recaptured in at least one environmental sample. The pale-yellow colour shows when a virus species was identified in poultry swabs but failed to be recaptured through environmental sampling. The blue colour shows virus species which were found in the environment but failed to be detected in poultry swabs. Rows of the heatmap represented virus species and columns represent individual environmental samples (annotated at the top of the heatmap). B Comparison of the alpha-diversity measure observed species between all ten groups included in this study (P value shown = 0.02 and n = 140). Lastly the number of viruses detected in each environmental sample is compared relative to (C) chicken oropharyngeal (n = 15; P = 0.004 versus air slaughter area and P = 0.007 versus air holding area), D chicken cloacal (n = 15), E duck oropharyngeal (n = 14; P value shown =0.01) and (F) duck cloacal (n = 15; P value shown = 0.04) swabs. Comparative environmental samples were collected from air (holding area, n = 13; slaughter area, n = 14; outside, n = 10), wash water (n = 19) drinking water (n = 9), and cage swabs (n = 17). Each group is represented with data from independent biological samples. Boxplots display the median (centre line), interquartile range (box), and whiskers extending to 1.5× the interquartile range; points represent individual samples. Statistics were calculated using a Kruskal-Wallis with Dunns post-hoc test. All P values obtained were corrected for false-discovery rate (FDR) using the Benjamini-Hochberg method. P values are annotated as follows: P  <  0.05 *; P  <  0.01 **; P  <  0.001***.

Fig. 2. Performance of environmental samples (ES) is dependent on poultry host species and sample type.

A Principal coordinate analysis (PCoA) of beta-diversity (Jaccard) of poultry virus species. PERMANOVA was calculated adjusting for both location and time as a confounder. The eigen values are reported showing the variation reported by PCo1 and PCo2. B Based on the PCoA in A, we calculated the median centroid of each group by determining the median PCo1 and PCo2 coordinates. Relative to poultry swabs, we calculated the distance of all samples from (C) chicken oropharyngeal (n = 15), (D) chicken cloacal (n = 15), (E) duck oropharyngeal (n = 14), and (F) duck cloacal (n = 15) median centroids, as shown in panel B. Comparative environmental samples were collected from air (holding area, n = 13; slaughter area, n = 14; outside, n = 10), wash water (n = 19) drinking water (n = 9), and cage swabs (n = 17). Each group is ranked from the shortest to longest distance from its respective poultry median centroid. Spearman correlations between virus abundance in reads per million (RPM) and the distance from (G) chicken oropharyngeal, H chicken cloacal (P value shown = 0.0003), I duck oropharyngeal (P value shown = 0.0001), and (J) duck cloacal (P value shown = 0.0001) median centroids were calculated for viruses either recaptured (n = 45) or undetected (n = 38) by environmental sampling, as shown in Fig. 1A. Each point represents an individual virus species, and therefore n values correspond to the number of distinct viruses rather than biological replicates. Boxplots display the median (centre line), interquartile range (box), and whiskers extending to 1.5× the interquartile range; points represent individual samples. K Ridge plot showing the probability of each poultry swab being the source of pathogenic viruses for all environmental samples. Statistics were calculated using a Kruskal-Wallis with Dunns post-hoc test. All P values obtained were corrected for false-discovery rate (FDR) using the Benjamini-Hochberg method. P values are annotated as follows: P  <  0.05 *; P  <  0.01 **; P  <  0.001***.

Fig. 3. Environmental contigs cover the majority of virus genes in poultry with a high sequence similarity.

A Presents a heatmap of detection rates for over 4,000 poultry virus reference genes (identified using DIAMOND blastx) in poultry swabs versus environmental samples. Red indicates that genes were found in both groups at the same timepoint, pale yellow represents detections only in poultry swabs, and blue for genes detected only in the environment. Rows correspond to individual virus genes annotated by genus (shown in the lower-right corner) and major clusters are labelled on the right. B For each gene (x axis) the difference in percent identity (PID) between poultry and environmental contigs is shown. C shows the number of poultry environment contig pairs at specified PID thresholds. D Distribution of percent identity (PID) values for environmental contigs aligning to the same reference genes as those detected in poultry. Boxplot displays the median (centre line), interquartile range (box), and whiskers extending to 1.5× the interquartile range; points represent individual samples. E Differences in percent coverage between poultry and environmental contigs are shown for each environmental sample type. Poultry samples comprised chicken oropharyngeal (n = 15), chicken cloacal (n = 15), duck oropharyngeal (n = 14), and duck cloacal (n = 15) swabs. Environmental samples were collected from air (holding area, n = 13; slaughter area, n = 14; outside, n = 10), wash water (n = 19), drinking water (n = 9), and cage swabs (n = 17). Each group represents independent biological samples. F For contig pairs aligned to the same gene(s), this plot displays the relative proportions of the major virus genera within each sample group. G For contigs that aligned uniquely to a single gene (unpaired), this plot shows the relative proportion contributed by each gene. Statistics were calculated using a Kruskal-Wallis with Dunns post-hoc test. All P values obtained were corrected for false-discovery rate (FDR) using the Benjamini-Hochberg method.

Fig. 4. Environmental sampling detected highly pathogenic influenza A contigs more often than poultry swabs.

A, B Maximum likelihood phylogenies of the hemagglutinin (HA) subtypes A) H5 and B) H9. C Heatmap showing the detection rate of different HA and neuraminidase (NA) subtypes across all twelve visits to the respective markets. Ward.d2 hierarchal clustering reveals HA and NA subtype combinations which co-circulated. Columns represent individual samples, annotated by location and sample type; rows represent HA and NA genes which were detected at LBMs. Red colour indicates that the influenza gene segment was detected, whereas white represents samples in which the segment failed to be detected. D–F Maximum likelihood phylogenies NA Influenza subtypes (D) N1, (E) N2, (F) N6. Each tip is coloured based on clade designation shown in the legend accompanying each phylogenetic tree.