Exploring retropharyngeal lymph node metabarcoding as a tool for characterizing microbial signatures in white-tailed deer from urban and non-urban areas

Introduction

Wildlife living in urban environments are exposed to a variety of novel stressors that can affect the establishment and composition of species communities (Faeth et al., 2011). For example, urban environments may expose wildlife to increased levels of chemical pollutants, higher temperatures, and pathogens (Bradley and Altizer, 2007, Sepp et al., 2019; Swank et al., 2021). Human land-use modification has also previously been associated with changes to disease prevalence, including chronic wasting disease (CWD) in mule deer (Farnsworth et al., 2005). In gut microbiomes, resident microbial communities can influence pathogen colonization through competition for resources, inhibitions, and interactions with host immune responses (Pickard et al., 2017). In white ibises, urban-associated shifts in gut microbiome composition have also coincided with lower microbial diversity, and individuals with lower diversity were more likely to shed Salmonella (Murray et al., 2020). Evaluating the differences in microbial community composition across urban and non-urban areas can therefore aid in identifying shifts associated with urbanization and potentially pathogenic taxa of concern for disease surveillance.

White-tailed deer (Odocoileus virginianus) are one species commonly found in both urban and non-urban areas of Florida. White-tailed deer in urban environments have demonstrated changes in behavior compared to their non-urban counterparts, with higher likelihood of temporal and spatial overlap with human activity (Péter et al., 2025). Urban deer can have overlapping home ranges, potentially increasing the rate of disease transmission (Walter et al., 2011). Additionally, simulations have suggested that fragmented landscapes can concentrate white-tailed deer within remaining habitat fragments and accelerate the transmission of CWD (Kjaer and Schauber, 2022), which has also been positively correlated with human land-use modification in mule deer (Farnsworth et al., 2005). Although previous research in Odocoileus has found limited differences in body condition and pathogen prevalence between urban and non-urban individuals (Mathieu et al., 2018), few studies have leveraged microbial community characterization across urban gradients to establish baseline community profiles and monitor shifts in microbial communities, including the detection of potentially pathogenic taxa that may require more focused disease surveillance. This monitoring can support more responsive pathogen-specific screening and also contribute to One Health initiatives focused on mitigating spillover events (Flint et al., 2015, Eby et al., 2023, Hayman et al., 2023).

Evaluation of bacterial community composition and diversity is often done through methods such as metabarcoding, which enables low-cost and rapid screening of diverse sample types with varying quality (Taberlet et al., 2012, Cabodevilla et al., 2024). This flexibility can be particularly useful for analyzing opportunistically collected tissues such as lymph nodes, which have previously been used to assess microbial community composition in harvested wildlife and commercial livestock, and to monitor food sources for pathogens prior to human consumption (Wittekindt et al., 2010, Mann et al., 2015, Queirós et al., 2019). Lymph nodes are also the primary sites of infection by various pathogens, including several Salmonella species (St John and Abraham, 2009). Because lymph nodes process lymph draining from tissues, microbial signatures may reflect exposure to local environmental factors (Hirosue and Dubrot, 2015). This makes them particularly valuable for wildlife and disease monitoring programs, such as CWD surveillance, where they are regularly collected and stored.

Here, we assessed microbial community structure in lymph nodes obtained from white-tailed deer in Florida. We performed 16S metabarcoding of individuals collected across urban and non-urban areas and assessed differences in community diversity across urban and non-urban groups. We found geographically correlated patterns of diversity between locations and identified taxa that may be of interest for disease management and monitoring. We also provide suggestions for improving the methodology of future studies that leverage lymph node sample collections for research in microbial community assessment or surveillance.

Methods and Materials

Sample collection & DNA extraction

We obtained retropharyngeal lymph nodes from 77 O. virginianus sampled by the Florida Fish and Wildlife Conservation Commission (FWC) as part of CWD screening in the state. These samples included individuals that were primarily identified as roadkill and hunter harvested, as well as those that died of other causes between the years 2019 and 2021, with lymph nodes later harvested and frozen. To assign individuals to urban or nonurban areas, we used the United States Census Bureau’s 2020 Urban Area shapefile, which delineates urban and metropolitan areas with high population densities and urban land use (U.S. Census Bureau, 2022).

We generated a buffer with a 1 km radius, equivalent to a home range area of approximately 3.14 km², around each occurrence using the Buffer geoprocessing tool in QGIS v3.40.14 (QGIS Development Team, 2026). We chose this buffer as its area falls within the range of published white-tailed deer home range estimates in Florida, which include a mean adult female range of 1.94 km² and mean annual adult male home range of 4.54-7.0 km² (Land et al., 1993, Sargent and Labisky, 1995). We calculated the proportion of urban land within each buffer based on its intersection with the U.S. Census Bureau urban area polygon (Figure 1).

Figure 1. Sampling locations of white-tailed deer in central and southern Florida. Pie charts show the proportion of urban (purple) and non-urban (green) land cover within a 1 km radius buffer around each sampling location. The locations of overlapping pies have been slightly adjusted for visualization.

Library preparation and sequencing

We extracted DNA via spin-column using a Quick-DNA Miniprep Plus Kit (Zymo Inc., Irvine, California, USA) per the manufacturer’s instructions, using an overnight digestion of a lymph node fragment with Proteinase K at 55°C. We targeted the V3-V4 hypervariable regions of the 16S rRNA following the Illumina 16S Metagenomic Sequencing Library Preparation workflow. We amplified genomic DNA in 25 μL PCR reactions containing 12.5μL KAPA HiFi HotStart ReadyMix (Kapa Biosystems Inc., Wilmington, Massachusetts, USA), 5 μL of forward and reverse primers with Illumina overhang adapters (forward: 5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGA
CAGCCTACGGGNGGCWGCAG-3’; reverse: 5’-GTCTCGTGGGCTCGGAGATGT
GTATAAGAGACAGGACTACHVGGGTATCTAATCC-3’) and 2.5 μL of DNA extract. We performed initial PCR under the following conditions: 90 °C for 3 min, 25 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and final extension at 72 °C for 5 min. PCR was followed by magnetic bead cleanups to remove primers and short fragments.

Sequencing libraries were created by attaching Illumina sequencing adapters and unique dual indices through indexing PCR. Reactions contained 25 μL KAPA HiFi ReadyMix, 10 μL water, 10 μL IDT for Illumina unique dual indices (Illumina Inc., San Diego, California, USA), and 5 μL of purified PCR product. Indexing PCR was run under the following conditions: 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and final extension at 72 °C for 5 min. We also purified indexed libraries with a second magnetic bead cleanup to remove residual adapters and primers.

We quantified final libraries using qPCR with the NEBNext Library Quant Kit for Illumina. We prepared serial dilutions of each library and quantified them using manufacturer-provided DNA standards using Illumina adapter-specific primers. We performed qPCR under the following conditions: initial denaturation at 95 °C for 1 min, 35 cycles of 95 °C for 15 s, and 63 °C for 45 s. We performed equimolar pooling of libraries prior to sequencing. Sequencing was performed on an Illumina MiSeq platform (v3 chemistry, 2 x 300 bp paired-end reads) at the University of Florida Interdisciplinary Center for Biotechnology Research (ICBR) NextGen DNA Sequencing Core.

Bioinformatic analyses

We processed raw FASTQ files using QIIME2 v2025.10 (Bolyen et al., 2019) and assessed sequencing quality and read depth across samples. Forward and reverse primer sequences and adapters were removed with qiime cutadapt trim-paired plugin (Martin, 2011). Amplicon sequence variants (ASVs) were inferred using the DADA2 algorithm in qiime dada2 denoise-paired (Callahan et al., 2016). We truncated forward and reverse reads to 234 bp and 229 bp, respectively, based on quality score distributions while maintaining a median quality score of 30 across the read length (Supplementary Material). We then filtered and merged paired reads using qiime dada2 denoise-paired to generate a ASV feature table and representative ASV sequences. We identified and removed contaminants using qiime quality-control decontam-identify, using negative controls defined in sample metadata. Contaminant filtering was performed using a decontam score threshold of 0.10, and negative controls were removed from the decontaminated feature table.

We assigned taxonomy to ASVs using the SILVA ribosomal RNA reference database (Chuvochina et al., 2025) with qiime feature-classifier classify-sklearn. We retained taxonomic assignments at the highest confidence supported by the classifier and filtered the feature table to remove sequences annotated as mitochondria or chloroplasts, restricting analyses to bacterial taxa. We also reconstructed phylogenetic relationships in ASVs prior to diversity analyses by aligning representative sequencing using MAFFT (Katoh and Standley, 2013) and constructing a phylogenetic tree using FastTree (Price et al., 2010). We generated rarefaction curves using qiime diversity alpha-rarefaction and selected a rarefaction cutoff of 1000 sequences per sample to balance sample retention and stabilization of diversity estimates across samples (Supplementary Material).

We evaluated alpha microbial diversity using the Shannon diversity index and observed features (ASV richness) and beta microbial diversity using Bray-Curtis dissimilarity and Jaccard distance across all rarefied samples. Significance of associations between alpha diversity and the proportion of urban land within the 1 km radius buffers was tested using Spearman correlations. For beta diversity, we evaluated associations based on the proportion of urban land using permutation multivariate analysis of variance (PERMANOVA) (Anderson, 2017) with 999 permutations. We also tested whether alpha diversity and beta diversity differed by sex and age class using Kruskal-Wallis tests. To assess spatial structure and isolation-by-distance, we generated pairwise geographic distance matrices from sample coordinates and performed Mantel tests to test correlations between geographic distance and community dissimilarity for Bray-Curtis and Jaccard distance matrices.

Results

The 77 sampled deer included 48 males and 29 females, which were primarily adults (N = 45), with vehicle collisions (N = 41) and hunter harvests (N = 18) as the primary causes of death (Table 1). We generated 7,250,147 paired-end read pairs across 77 samples after adapter trimming, averaging 94,158 read pairs per sample. Following DADA2 quality filtering, paired-end merging, and removal of contaminant reads and negative controls, 5,353,462 reads remained across 72 samples, representing 3,805 ASVs with lengths ranging from 234-451 nucleotides. Resulting sample sequencing depth had a median frequency of 8,177 (IQR 1348-33,666). Of the 3,805 ASVs, 2,709 remained after removal of mitochondrial and chloroplast sequences. At the phylum level, the taxa with the highest prevalence across the 72 samples were Firmicutes (N = 70), Proteobacteria (N = 70), Actinbacteriota (N = 53), Bacteroidota (N = 49), and Patescibacteria (N = 21). At the genus level, the taxa with the highest prevalence included Pseudomonas (N = 61), Psychobacter (N = 50), Escherichia-Shigella (N = 50), Streptococcus (N = 47), Acinetobacter (N = 42), Cutibacterium (N = 40), Vagococcus (N = 33), Clostridium sensu stricto (N = 30), Staphylococcus (N = 28), and Carnobacterium (N = 26). Other notable taxa that were present included Erysipelothrix (N = 23) and Bartonella (N = 20).

Table 1. Sex, cause of death, and age classes of samples white-tailed deer grouped by urban classification, with urban samples described as those having at least 50% urban land within the 1 km radius buffer around collection locations.

After rarefaction to 1000 sequences per sample, 50 (69.4%) of the 72 samples were retained, including 23 from urban areas and 27 from non-urban areas. The mean number of observed features (ASV richness) was 41.5 (σ = 24.8) for urban samples and 53.6 (σ = 46.1) for non-urban samples, and there was no evidence of differentiation between groups (H = 0.27, p = 0.61). Mean Shannon diversity was 3.47 (σ = 1.23) for urban and 3.19 (σ = 1.54) for non-urban, with no differences between groups (H = 0.72, p = 0.40). No patterns or clustering associated with the proportion of urban land or other metadata were apparent in Bray-Curtis or Jaccard ordinations (Supplementary Material). PERMANOVA showed no significant association between the proportion of urban land and community composition for Bray-Curtis (F = 0.73, R2 = 0.01, p = 0.89) and Jaccard (F = 1.12, R2 = 0.02, p = 0.11). Additionally, we did not detect significant differences in alpha or beta diversity by sex (Bray-Curtis: F = 1.01, R2 = 0.02, p = 0.42; Jaccard: F = 0.90, R2 = 0.2, p = 0.87) or age class (Bray-Curtis: F = 0.68, R2 = 0.02, p = 0.99; Jaccard: F = 0.96, R2 = 0.03, p = 0.74). Results of our Mantel test revealed significant associations between geographic distance and community dissimilarity based on both Bray-Curtis (r = 0.26, p < 0.01) and Jaccard (r = 0.25, p < 0.01) beta diversity, and there was some geographical clustering based on principal component ordination, particularly near the Tampa metropolitan urban area (Figure 2).

Figure 2. Geographic distribution of variation in bacterial community composition in white-tailed deer samples (following rarefaction) from central and southern Florida. Sampling locations are colored by scores along the first principle coordinate axis (PCoA1) from Jaccard dissimilarities (A) and Bray-Curtis dissimilarities (B). Samples with similar colors have more similar bacterial community composition, whereas contrasting colors indicate greater differences. The locations of overlapping points have been slightly adjusted for visualization.

Discussion

We used lymph nodes collected from CWD monitoring to evaluate differences in community composition in deer from urban and non-urban areas in Florida. Although we were unable to detect differentiation using alpha and beta diversity statistics, we found that patterns of beta diversity were geographically correlated using the Bray-Curtis and Jaccard indices. These results suggest that urbanization alone, when represented using broadly defined urban areas, may not be the only driver of lymph node-associated microbial community diversity across our sampled deer, and that the observed variation may be due to other processes, including more localized spatial or individual-level processes. The apparent geographical clustering near the Tampa metropolitan area, but lack of significant differentiation in diversity statistics, suggests that bacterial diversity may be driven more by a specific urban area’s population dynamics than by general urban versus non-urban patterns. Habitat fragmentation can concentrate deer within habitat patches and increase overlap among individuals (Water et al., 2011, Kjaer and Schauber, 2022), and we may not be able to detect variation in small urban habitat patches or at finer spatial scales where interactions might take place. Additionally, local differences in deer density, movement, and contact rates can also contribute to the observed clustering independently of, or in addition to, the proportion of urban land cover. We classified several taxa at the genus level, and the most common across samples were Pseudomonas, Psychrobacter, Escherichia-Shigella, Streptococcus, and Acinetobacter. Pseudomonas, Psychrobacter, and Acinetobacter are broad genera that include both pathogenic and non-pathogenic environmental and host-associated taxa (Wong et al., 2016, Welter et al., 2021, Girard et al., 2021). Escherichia-Shigella is a potential pathogen commonly found in gut microbiomes (Li et al., 2018, 2022), and its detection in lymph node tissue may reflect transfer during handling of deer carcasses rather than a lymph node-specific signal. This may be particularly relevant for deer collected after a vehicle collision, which may remain at the site of mortality for an unknown period prior to collection and have greater environmental exposure and leakage of bodily fluids. Streptococcus is common in vertebrate oral and upper respiratory microbiomes, but the genus also encompasses both commensal and pathogenic lineages (Abranches et al., 2018). Other taxa with potential relevance to disease monitoring included Erysipelothrix and Bartonella. Erysipelothrix includes the ungulate pathogen Erysipelothrix rhusiopathiae, which has been associated with widespread mortality in northern ungulates (Forde et al., 2016, Mavrot et al., 2020). Bartonella is a bloodborne vector-transmitted genus frequently detected in cervids and includes zoonotic species that can cause bacteremia in humans (Duodu et al., 2013, Razanske et al., 2018). Although these results provide a descriptive overview of retropharyngeal lymph node-associated bacterial taxa, interpretation of the function or impact of these taxa on communities of white-tailed deer should be approached cautiously without species- or strain-specific identification.

Disease monitoring by wildlife agencies can provide an opportunistic source of retropharyngeal lymph node samples and similar tissues, and we offer four recommendations for researchers interested in leveraging these for surveillance and monitoring. First, metabarcoding using multiple markers can broaden detection and improve taxonomic resolution. Standard 16S rRNA metabarcoding of the V3-V4 hypervariable region primarily targets bacterial DNA and is not intended to detect fungi, eukaryotic microorganisms, or viruses (Bukin et al., 2019). We recommend including additional metabarcoding targets such as fungal ITS (Beeck et al., 2014) and eukaryotic 18S rRNA (Hadziavdic et al., 2014) regions, as well as viral metabarcoding or metagenomic approaches (Gil et al., 2021). Second, QIIME databases such as SILVA are useful for broad taxonomic profiling, though we encourage researchers to develop custom or disease-specific databases to provide more relevant characterization for monitoring beyond a general community overview. These databases typically include curated sequences with taxonomic annotations and can be tailored to detect species or taxa of particular concern, for example through tools such as RESCRIPt, which support the development and creation of custom reference databases (Robeson et al., 2021). Third, sample collection and storage may be a major source of variation, especially for hunter-harvested samples and field collection with resource limitations. When possible, we recommend preservation by freezing or ethanol storage to mitigate post-mortem bacterial growth. We also encourage researchers to partner with agencies to engage hunters as collaborators by communicating the value of these samples for disease monitoring, which can promote more consistent handling, confidence in results, and help support healthy game populations through proactive disease management. Fourth, future work should focus on sampling multiple habitat patches within specific urban areas where managers or researchers have knowledge of local areas that support deer, to improve detection of patterns that may be missed with broad characterizations of urbanization. Habitat patches within urban and suburban areas can form mosaics that create complex dynamics relevant to potential disease transmission. Expanding the number of sampling sites in urban areas would improve characterization of microbial communities, particularly as deer can occur within major cities where suitable habitat patches remain (Ciach and Fröhlich, 2019, Potratz et al., 2019).

Conclusions

In this study, we characterized bacterial communities in retropharyngeal lymph nodes collected through CWD surveillance to evaluate patterns associated with urbanization and geography in Florida deer. We found no differences between urban and non-urban groups in alpha or beta diversity. However, geographic structuring in beta diversity suggests local spatial patterns. These results highlight the utility of wildlife disease surveillance programs, such as CWD screening, as an opportunistic source of tissues for research and monitoring beyond their original intent. Using these tissues for long-term monitoring may also provide a proactive means of detecting shifts in diversity associated with changing disease risk or pressure as part of an integrated wildlife monitoring approach (Barroso et al., 2023). We recommend that researchers use multi-marker approaches to broaden detection across taxa, develop disease monitoring-relevant databases for characterization, preserve samples as consistently as feasible to minimize post-mortem changes, and prioritize sampling in specific urban centers to better understand local patterns of diversity.