## 2 Materials and Methods
### Sequencing and Phylogenetic Analyses  
Sequencing of extracted DNA was performed at MR DNA (www.mrdnalab.com, Shallowater, TX, USA on 18 February 2019 and 8 July 2019) using MiSeq sequencing following the manufacturer's guidelines. Sequence data were processed using MR DNA's analysis pipeline (MR DNA, Shallowater, TX, USA). The primers 515F (5'-GTGYCAGCMGCCG TAA-3') and 806R (5'-GGACTACNVGGGTWTCTAAT-3') were used to target the V4 region of the 16S rRNA gene of bacteria. These primers were used in a single-step, 30-cycle PCR using the HotStarIaq(tm) Plus Master Mix Kit (Qiagen, Germantown, MD, USA) as follows: 94 \({}^{\circ}\)C for three minutes, followed by 30-35 cycles of 94 \({}^{\circ}\)C for 30 s, 53 \({}^{\circ}\)C for 40 s, 72\({}^{\circ}\)C for one minute, and a final elongation step at 72 \({}^{\circ}\)C for five minutes. After amplification, PCR products were checked in a 2% agarose gel to determine the success of amplification and the relative intensity of bands. Pooled samples were purified using calibrated Ampure XP beads (BECKMAN COULTER Life Sciences, Indianapolis, IN, USA). The pooled and purified PCR product was used to prepare the Illumina DNA library.  
A Quantitative Insights into Microbial Ecology 2, version 2022.11 (QIIME2) [33], pipeline was implemented using the SILVA release 138 database ([[https://www.arb-silva.de/documentation/release-138/](https://www.arb-silva.de/documentation/release-138/)]([https://www.arb-silva.de/documentation/release-138/](https://www.arb-silva.de/documentation/release-138/))) [34], as well as DADA2, version 1.10.0 [35], to prepare the data for analysis. QIIME2 was used to demultiplex samples and assign barcodes based on sample identification. DADA2 was used to separate demultiplexed samples into amplicon sequence variants (ASVs), and the primers were removed. ASVs were selected over OTUs for reasons outlined in [PERSON], et al. [36]. Visualizations and artifacts were then produced by QIIME2 in order to describe the data. Sampling depths were set to 40,000 features per sample for all bacterial analyses. Subsequent analyses were performed in rStudio, version 2023.09.1+494, using vegan (version 2.6-6.1), fossil (version 0.4.0), ggplot2 (version 3.5.1), reshape2 (version 1.4.4), phyloSeq (version 1.48.0), metagenomeSeq (version 1.46.0), and DESeq2 (version 1.44.0) packages [37; 38; 39; 40; 41; 42; 43; 44; 45], as well as Microsoft Excel (version 16.86).  
Statistical significance in alpha diversity was determined using a non-parametric, pairwise, Wilcoxon test by means of the function wilcox.test under the core R stats package. The alpha was set to 0.05, and we rejected the null hypothesis (i.e., no statistical difference between the two datasets existed) in favor of the alternate hypothesis (i.e., there existed a statistical significance between the two datasets), if the comparison yielded a \(p\)-value greater-than or equal-to our predetermined alpha.  
Statistical testing in our beta diversity was performed using a PERMANOVA, as outlined by [PERSON] [46; 47]. Multivariate homogeneity of group distances was calculated using the R package, \"vegan\", and the function \"betadisper\" as well as the functions \"anova\" and \"TukeyHSD,\" to test for homogeneity of dispersions, as the design was unbalanced (unequal number of samples in each group). Results exhibiting homogeneity of group distances were then tested with a PERMANOVA test, using the adonis2 function in vegan. Results were considered significant when \(p\)-values were at or below our alpha (0.05).