PICRUSt2 metagenome prediction

16S rRNA gene sequence data processed using the software package mothur (v.1.40.2) [120,121] was used for metagenome prediction analysis using the PICRUSt2 pipeline [59]. 16S rRNA gene sequences were aligned using HMMER (v.3.3, hmmer.org), placed in the default 16S rRNA gene reference tree, which is comprised of 20,000 16S rRNA gene sequences in the Integrated Microbial Genomes database [127], using EPA-NG [128], and then the complete tree was constructed with GAPPA [129]. Following tree construction, unknown lineages were inferred, and KEGG pathway copy number was predicted using castor [130]. Finally, metabolic pathway abundances were inferred from MetaCyc using MinPath [131]. For analysis, metabolic pathway abundances were rounded to the nearest whole number and normalized across each sample to determine the relative abundance of each gene family and inferred metabolic pathway. Principle coordinate analysis was performed in R (v.3.6.3) using the “ggplot2” and “ggfortify” packages to visualize differences in inferred metabolic pathway relative abundances between experimental groups. AMOVA P values [124] were calculated using the “vegan” package and used to determine statistically significant differences between experimental groups [47].