Effects of high dietary zinc supplementation timing on the biological responses of gestating sows and their piglets
Data files
Nov 01, 2023 version files 64.89 KB
-
Dataset_S1.csv
-
Dataset_S2.csv
-
Dataset_S3.csv
-
Dataset_S4.csv
-
README.md
Abstract
High dietary zinc (Zn) fed to gestating sows may have utility as a fetal imprinting strategy to decrease pre-weaning mortality of piglets. However, the biological action that Zn may work through is unknown. High Zn may modulate the microbiome of the sow and the microbial seeding of the offspring’s gut microbiome. Additonally, high dietary Zn and piglet birth weight may alter gene expression in piglet whole blood. Sows (n = 267) were fed 1 of 3 dietary treatments: 1) Control: a corn-soybean meal-based diet containing 125 ppm total supplemental zinc, 2) Breed-to-Farrow: as Control + 141 ppm supplemental Zn as ZnSO4 fed from 5 days post-breeding to farrowing; and 3) Day 110-to-Farrow: as Control + 2,715 ppm supplemental Zn as ZnSO4 starting on day 110 of gestation until farrowing. A subset of third parity sows (n = 30) were selected to assess the microbiome of colostrum, milk, and rectal and vaginal surfaces of sows. At farrowing, 4 pigs per litter (n = 120) were selected based on birthweight (BiW), as 2 average BiW pigs and 2 pigs with BiW below the litter average were selected for assessing the piglet gut microbiome on the day of birth (day 0) and day 5 of age. 16S rRNA sequencing were implemented for milk and colostrum samples while all other sample types were sequenced using shotgun metagenomics to determine taxonomic and functional profiles. On a different subset of pigs, whole blood was collected from 9 LBW pigs per treatment and 8 ABW Control pigs for RNA-sequencing to evaluate differentially expressed genes (DEGs) and pathways. Only 2 to 3 genes were differentially expressed between Control LBW and LBW pigs born to sows fed high Zn. However, 262 DEGs were identified when comparing LBW and ABW pigs, mostly reflecting pathways associated with translation, ribosome biogenesis, and amino acid and protein synthesis. Measures of alpha diversity (richness and Shannon’s H Index) and beta diversity (Bray-Curtis, PERMANOVA) were conducted along with indicator species analyses. Species with an indicator value of > 0.50 were confirmed with a generalized linear mixed model as each P-value was corrected for false discovery rate (FDR) to generate a Q-value. For piglet samples, the MaAsLin2 R package was used to determine multivariate associations between dietary treatment and piglet BiW. High dietary concentrations of Zn fed to gestating sows did not affect the colostrum, milk, or vaginal microbial diversity or populations of sows. Pathogenic bacteria such as Shigella flexneri and Salmonella enterica were less abundant in fecal samples from Breed-to-Farrow sows compared to Control sows. For piglets born to Breed-to-Farrow sows, their gut microbiome favored fiber fermenting, short chain fatty acid generating microbial species compared to Control pigs. Day 110-to-Farrow piglets demonstrated a lower abundance of SCFA producing bacteria compared to Control piglets. Gene families and pathways playing roles in central metabolic functions (starch, pyruvate, sucrose, amino acid metabolism) were more abundant in Breed-to-Farrow piglets compared to pigs born to Control sows. In conclusion, high Zn fed to gestating sows may influence SCFA-producing species and may reduce the abundance of potential pathogenic bacteria in the sow and piglet. Piglet birth weight may have greater effects on gene expression of neonatal pigs.
README: Effects of high dietary zinc supplementation timing on the biological responses of gestating sows and their piglets
https://doi.org/10.5061/dryad.q573n5tq5
The original work was published by Hammers et al (2023). In addition to this work, supplementary datasets, figures, and tables are included here to describe the effects of high zinc supplementation to gestating sows on gene expression in piglet whole blood and on the modulation of the sow and piglet microbiome. These datasets are comprised of the summary statistics regarding 16S rRNA and shotgun metagenomic sequencing and data processing for colostrum, milk, vaginal, and fecal microbiomes of sows fed high dietary zinc and the intestinal microbiomes of their piglets. Additionally, supplemental figures and tables are included from data analysis. For gene expression data, RNAseq sequencing data, all differentially expressed genes based on piglet birth weight, and a supplementary pathway analysis figure is included.
Description of the data and file structure
Supplemental Datasets includes 16S rRNA sequencing data for colostrum and milk samples collected from sows feed high dietary concentrations of zinc. Supplemental Datasets also includes Shotgun metagenomic sequencing data for vaginal and fecal samples collected from sows and fecal samples collected from piglets on day 0 and day 5 of age. RNA sequencing data are also included in a Supplemental Dataset along with differentially expressed genes based on piglet birth weight. Supplementary figures and tables were generated during statistical analysis to determine differential gene expression and pathway analysis and statistical analysis of microbiome data.
Datasets include:
- Effects of high zinc supplementation to gestating sows on piglet whole blood gene expression:
- Dataset S1: Summary of RNA sequencing data
- Dataset S2: List of differentially expressed genes comparing piglet birth weight of low and average birth weight
- Dataset includes gene Feature ID, gene names, raw read counts for low and average birth weight (kg) pigs, the fold change, Bonferroni correction and FDR p-value corrections generated from the EdgeR statistical analysis.
- Figure S1: Piglet birth weight pathway analysis of differentially expressed genes. Located in Zenodo.
- Effects of high zinc supplementation to gestating sows on the microbiome of sows and piglets
- Dataset S3: Summary of 16S rRNA sequencing data for milk and colostrum samples
- Treatment column corresponds to the dietary treatment assigned to each sow: A) Control diet containing 125 ppm Zn, B) Breed-to-Farrow containing 266 ppm Zn fed 5 days after breeding until farrowing, and C) Day 110-to-Farrow treatment which was the Control diet fed until day 110 of gestation and then 2,840 ppm supplemental Zn was fed until farrowing.
- Samples with "NA" in the treatment column are negative control samples and do not have an assigned treatment
- Dataset S4: Summary of shotgun metagenomic data for sow vaginal, fecal, and piglet fecal samples
- Figure S2: Microbial diversity of sow vaginal samples
- S2a-c: Shannon's H index, microbial richness, and beta diversity
- S2d: Discriminant species of sow vaginal samples
- Figure S3: Diversity of sow vaginal gene families
- S3a-c: Shannon's H index, gene family richness, and beta diversity
- S3d: Discriminant gene family of sow vaginal samples
- Figure S4: Microbial diversity of sow fecal samples
- S4a-c: Shannon's H index, microbial richness, and beta diversity
- Figure S5: Diversity measures of sow fecal gene families
- S5a-b: Shannon's H index and gene family richness
- S5c: Community composition of gene families in fecal samples of sows
- Figure S6: Microbial diversity measures of colostrum samples
- S6a-c: Shannon's H index, microbial richness, and beta diversity
- S6d: Discriminant family identified in colostrum samples
- Figure S7: Microbial diversity measures of milk samples
- S7a-b: Shannon's H index and microbial richness
- S7c: Community composition of milk samples
- Figure S8: Discriminant species identified for low birth weight (LBW) compared to average birth weight piglets (ABW)
- Figure S9: Diversity measures of pig fecal sample gene families at day 0 of age
- S9a-b: Shannon's H index and gene family richness
- S9c: Community composition of gene families of pig day 0 fecal samples
- Figure S10: Diversity measures of pathways in piglet fecal samples at day 0 of age
- S10a-b: Shannon's H index and richness of pathways
- S10c: Pathway beta diversity of piglet samples at day 0 of age
- Figure S11: Discriminant pathways identified for Breed-to-Farrow piglets at day 0 of age compared to Control pigs
- Figure S12: Diversity measures of piglet gene families at day 5 of age
- S12a-b: Shannon's H index and richness
- S12c: Community composition of gene families at day 5 of age for pigs
- Figure S13: Diversity measures of pathways in piglet fecal samples at day 5 of age
- S13a-b: Shannon's H index and pathway richness
- S13c: Community composition of pathways at day 5 of age for pigs
- Figure S14: Discriminant pathways identified for Day 110-to-Farrow piglets at day 5 compared to Control pigs
- Table S1: Discriminant pathways in sow vaginal samples
- Table S2: Discriminant pathway in sow fecal samples
- Table S3: Discriminant pathway in piglet fecal samples at day 0
- Table S4: Discriminant pathway in Breed-to-Farrow piglet fecal samples at day 5
- Dataset S3: Summary of 16S rRNA sequencing data for milk and colostrum samples
- All listed figures and tables are located in Zenodo.
Sharing/Access information
Data is originally published and supplemental data corresponds to:
- Hammers, K. L. 2023. Feeding high zinc to gestating sows and implications on sow and piglet performance, their biological responses, and the environment. Ph.D. Dissertation. University of Minnesota, St. Paul, MN.
Code/Software
RNA reads were checked for quality and trimmed using Trimmomatics (Bolger et al., 2014) and FastQC. High quality reads were aligned to the Sus scrofa genome using HISAT2 (Zhang et al., 2021). Over-represented globulin transcripts (Hbb and Hbz) were manually removed prior to data normalization. Mapped RNAseq data were analyzed using the edgeR (Robinson et al., 2010) statistical software within the R interface to determine differentially expressed genes using the criteria of an absolute fold change of > 2 and a False Discovery Rate (FDR) of < 0.05. Comparisons between Control LBW pigs and the LBW pigs born to sows fed high Zn (Breed-to-Farrow and Day 110-to-Farrow) and Control LBW pigs and Control ABW pigs were evaluated, which resulted in the generation of three datasets highlighting differentially expressed genes based on sow dietary Zn treatment and pig birth weight. After analysis of differential gene expression, ShinyGo (Ge et al., 2020) and gProfiler (Raudvere et al., 2019) were used to conduct a pathway enrichment analysis based on differentially expressed genes identified for the pairwise comparison between birth weights categories.
Sequence data were generated targeting the V4 variable region of the 16S rRNA gene on the MiSeq sequencing platform (MiSeq 2×300bp sequencing lane) using the primers 515F (59-GTGCCAGCMGCCGCGGTAA-39) and 806R (59-GGACTACHVGGGTWTCTAAT-39). Copy numbers of the 16S rRNA bacterial gene were quantified by qPCR using Kapa HiFi polymerase (Kapa Biosystems, Woburn, MA) through 25 cycles. Perl scripts to process raw sequence data created by Law et al. (2021) were used to remove primer sequences and filter low quality reads. Processed sequences were then run through the QIIME2 pipeline (Bolyen et al., 2019) and assigned amplicon sequence variants (ASVs) using the DADA2 plug-in (Callahan et al., 2016) and the Greengenes database, version 13_8 (McDonald et al., 2012).
Shotgun metagenomic sequencing was applied to vaginal and fecal samples. The Kneaddata pipeline (version 0.12.0; McIver et al., 2018) and FastQC (Babraham Bioinformatics) were used for sequence data pre-processing. Briefly, raw reads were then trimmed using Trimmomatic (Bolger et al., 2014) based on a sliding window approach where reads were cut if the average base Phred quality score within a four-base sliding window dropped below 20. Reads were also discarded when the length of the read was shorter than 90 base pairs. Remaining high quality reads were mapped against the Ensembl (Sscrofa11.1) Sus scrofa reference genome (Warr et al., 2020) to identify and remove swine contaminated reads using Bowtie2 (version 2.5.1; Langmead and Salzberg, 2012). Bowtie2 was used with the default parameters and decontaminate-pairs alignment. The remaining high quality reads were used for taxonomic classifications with the kraken2-bracken (Lu et al., 2017; Wood et al., 2019) pipeline and functional profiles using HUMAnN3 (Beghini et al., 2021). Taxonomic and functional profiles were generated in relative abundances from phylum to species level. HUMAnN3 output data were transformed to relative abundance and used the Unifrac90 reference protein database.
Methods
Experimental design
This experiment was conducted on a commercial sow farm (2,500 sows) located in the Upper Midwest. Three consecutive weekly farrowing groups that included 267 females (parity 0 to 6; PIC Camborough, Hendersonville, TN) were assigned randomly to 1 of 3 dietary treatments within parity about 5 days post-breeding. Microbial samples were collected from a subset of third parity sows (n = 30) and were selected based on parity to prevent the confounding effects of parity on performance response criteria and sample analyses. After farrowing, 4 piglets per litter were selected based on average litter birth weight for fecal sampling at days 0 and 5 of age. For each subsampled sow, 2 low birth weight (LBW) pigs (below the average litter birth weight of the respective litter) and 2 average birth weight (ABW) pigs were selected.
Dietary treatments consisted of: 1) Control – sows fed a corn-soybean meal-based diet containing 125 ppm total supplemental zinc supplied by zinc hydrochloride (Intellibond Z, Micronutrients, Indianapolis, IN); 2) Breed-to-Farrow – as Control + 141 ppm supplemental Zn as ZnSO4 fed from 5 days post-breeding to farrowing; and 3) Day 110-to-Farrow – as Control + 2,715 ppm supplemental Zn as ZnSO4 starting on day 110 of gestation until farrowing. Final supplemental Zn concentrations of the 3 dietary treatments were: 1) Control – 125 ppm; 2) Breed-to-Farrow – 266 ppm; and 3) Day 110-to-Farrow – 2,840 ppm. The treatments were imposed by feeding 45 g of the Breed-to-Farrow topdress or 80.5 g of the Day 110-to-Farrow topdress in addition to the gestation and lactation diets. These topdresses were formulated to contain the same total amount of Zn either fed for the duration of gestation or 5 days before farrowing. Control sows did not receive any topdresses. After farrowing, all sows were fed a standard lactation diet.
During gestation, sows were housed in individual stalls on partially slatted concrete floors. Treatments were assigned to a block of gestation stalls to avoid cross-contamination of treatment diets among adjacent sows; thus a “buffer” sow was placed at the end of each block to receive the same dietary treatment but was not included in the statistical analysis. At approximately day 110 of gestation, sows were moved to a separate farrowing stall until litters were weaned. Before sows entered the farrowing room, rooms were power washed, disinfected, and allowed adequate time to dry before sows were loaded into farrowing stalls within each room. Each farrowing room contained 39 stalls and was equipped with 1 stainless steel feeder and 1 nipple waterer on a fully slatted floor over a deep manure collection pit. Within 12 hours of birth, piglets (n = 3,990) were weighed, and ear tagged (LeeO, PrairiE Systems, Spencer, IA). Piglets were weighed again the day before weaning and pigs were weaned about 19 ± 6 days of age.
Piglet Whole Blood Collection for RNAseq and RNA Extraction
Whole blood was collected from 9 LBW pigs and 8 ABW pigs born to Control sows. Blood was also collected from 7 LBW piglets born to Breed-to-Farrow sows and 8 LBW piglets from Day 110-to-Farrow sows. Whole blood was collected from piglets to assess the transcriptome profile of offspring from sows fed high dietary Zn. Additionally, the effects of birth weight on differential gene expression (DGE) were also measured. Pigs were restrained and blood (2.5 mL) from each piglet was collected into PAXgene Blood RNA tubes (PreAnalytiX, Hombrechtikon, Switzerland) from the anterior vena cava. Blood was collected with PAXgene tubes because RNA degradation is minimized and the RNA expression profiles of samples are preserved to allow for accurate analysis of gene expression. Each tube was then inverted 8 to 10 times to facilitate adequate mixing of blood with tube contents (proprietary agent for intracellular RNA stabilization). According to instructions by the manufacturer, blood was incubated for 2 hours at room temperature to allow complete lysis of blood cells. Blood samples were then frozen at -20ºC for 24 hours according to the manufacturer and then placed on dry ice and stored at -80ºC until RNA extraction. Tubes were allowed to thaw for 2 hours and RNA was extracted manually using the PAXgene Blood RNA Kit (IVD) following manufacturer’s instructions. After extraction, RNA quality and quantity was assessed using a NanoDrop 2000 spectrophotometer (NanoDrop, ThermoFisher Scientific, Waltham, MA). Integrity of extracted RNA was evaluated by the University of Minnesota Genomics Center (UMGC; St. Paul, MN) and RNA integrity (RIN) was determined for each sample, which yielded an average RIN score of 9.0. Samples were submitted to the UMGC where 32 unique dual-indexed TruSeq Stranded mRNA libraries were created and sequenced using the NovaSeq sequencing platform. A total of 1,867,364,569 sequence reads were generated and samples averaged 29,177,571 reads ± 3,346,237 with an average sequence quality (Q) score of 36.
Statistical Analysis
RNA reads were checked for quality and trimmed using Trimmomatics (Bolger et al., 2014) and FastQC. High-quality reads were aligned to the Sus scrofa genome using HISAT2 (Zhang et al., 2021). Over-represented globulin transcripts (Hbb and Hbz) were manually removed prior to data normalization. Mapped RNAseq data were analyzed using the edgeR (Robinson et al., 2010) statistical software within the R interface to determine differentially expressed genes using the criteria of an absolute fold change of > 2 and a False Discovery Rate (FDR) of < 0.05. Comparisons between Control LBW pigs and the LBW pigs born to sows fed high Zn (Breed-to-Farrow and Day 110-to-Farrow) and Control LBW pigs and Control ABW pigs were evaluated, which resulted in the generation of three datasets highlighting differentially expressed genes based on sow dietary Zn treatment and pig birth weight. After analysis of differential gene expression, ShinyGo (Ge et al., 2020) and gProfiler (Raudvere et al., 2019) were used to conduct a pathway enrichment analysis based on differentially expressed genes identified for the pairwise comparison between birth weights categories.
Microbial sample collection from sows and piglets
Rectal and vaginal surfaces were sampled the day before expected farrowing date for each sow. Briefly, a sterile swab was used to scrape the inner rectum wall of each sow and then placed in a sterile 5ml tube. For vaginal samples, the outer vaginal area was cleaned with sterile gauze and a PBS (phosphate buffered saline) solution to remove any fecal material before sample collection. Subsequently, a sterile swab dipped into PBS was inserted approximately 2 inches into the vaginal opening avoiding contact with other skin surfaces. To collect colostrum and milk, each sow’s udder and sample collector’s gloves were disinfected with a rubbing alcohol wipe to prevent contamination of samples with environmental bacteria. Sow’s teats were hand-stripped to collect colostrum within 24 hours of the onset of farrowing and milk on day 2 of lactation. Samples were collected directly into a sterile 15 mL tube. An intramuscular injection of oxytocin (10 IU) was administered to sows for the collection of milk samples. Piglet fecal samples were collected from 4 piglets on day 0 before cross-fostering and 5 postpartum with a sterile swab as described for sows. All samples were placed immediately on dry ice after collection until they could be frozen at -80ºC prior to DNA extraction.
DNA extraction and sequencing
Microbial DNA was extracted from sow and piglet samples using the Qiagen PowerSoil DNA extraction kits (QIAGEN, Venlo, Netherlands) for all samples. Prior to DNA extraction, colostrum and milk samples were spun at 3,000 × g for 20 min to form a pellet. The subsequent pellet formed after centrifugation was collected on a sterile cotton swab and used for extractions. For milk and colostrum samples, 16S rRNA sequencing was implemented to concentrate microbial DNA as milk is a relatively low microbial biomass matrix and rich in host DNA (Perez et al., 2007; Cheema et al., 2021). Sequence data were generated targeting the V4 variable region of the 16S rRNA gene on the MiSeq sequencing platform (MiSeq 2×300bp sequencing lane) using the primers 515F (59-GTGCCAGCMGCCGCGGTAA-39) and 806R (59-GGACTACHVGGGTWTCTAAT-39) and dual-indexing library preparation (Gohl et al., 2016). Briefly, copy numbers of the 16S rRNA bacterial gene were quantified by qPCR using Kapa HiFi polymerase (Kapa Biosystems, Woburn, MA) through 25 cycles. Perl scripts to process raw sequence data created by Law et al. (2021) were used to remove primer sequences and filter low-quality reads. Raw sequence data contained an average of 81,485 ± 35,763 forward/reverse reads per sample (range: 1,095 to 136,182 reads/sample), which was reduced to an average of 60,909 ± 31,639 reads per sample (range: 330 to 114,806 reads/sample) after processing and quality control procedures. Processed sequences were then run through the QIIME2 pipeline (Bolyen et al., 2019) and assigned amplicon sequence variants (ASVs) using the DADA2 plug-in (Callahan et al., 2016) and the Greengenes database, version 13_8 (McDonald et al., 2012).
Shotgun metagenomic sequencing was applied to vaginal and fecal samples. The Kneaddata pipeline (version 0.12.0; McIver et al., 2018) and FastQC (Babraham Bioinformatics) were used for sequence data pre-processing. Briefly, raw reads were then trimmed using Trimmomatic (Bolger et al., 2014) based on a sliding window approach where reads were cut if the average base Phred quality score within a four-base sliding window dropped below 20. Reads were also discarded when the length of the read was shorter than 90 base pairs. Remaining high-quality reads were mapped against the Ensembl (Sscrofa11.1) Sus scrofa reference genome (Warr et al., 2020) to identify and remove swine contaminated reads using Bowtie2 (version 2.5.1; Langmead and Salzberg, 2012). Bowtie2 was used with the default parameters and decontaminate-pairs alignment. The remaining high-quality reads were used for taxonomic classifications with the kraken2-bracken (Lu et al., 2017; Wood et al., 2019) pipeline and functional profiles using HUMAnN3 (Beghini et al., 2021). Taxonomic and functional profiles were generated in relative abundances from phylum to species level. HUMAnN3 output data were transformed to relative abundance and used the Unifrac90 reference protein database. Raw sequence data for sow fecal, sow vaginal, and piglet fecal samples are presented in Table 5.1, and sequence data after data pre-processing steps utilizing the Kneaddata pipeline are presented in Table 5.2 of the manuscript.
DNA extraction kits possess a distinct collection of microbes and have been defined as the “kitome” (Salter et al., 2014). Negative controls (n = 17) consisting of a sterile swab and each kit’s reagents were added during DNA extraction for the identification and subsequent removal of potential reagent or environmental contamination. Negative control samples (n = 13) extracted to be compared against metagenomic data did not produce any genomic material identified during library prep and thus did not produce a sufficient library to be sequenced.
Statistical analyses
All statistical analyses of microbiome data were performed using the R statistical interface (R Core Team, 2014). Negative controls created from DNA extraction kit reagents and sterile swabs were used to screen ASV-level sequence data for potential contamination using the R decontam package (Davis et al., 2018) for the milk/colostrum 16S rRNA data. This method is based on both the presence or absence of bacterial taxa in samples compared to the corresponding control samples and for the frequency of appearance. The ASVs identified as contaminants were removed after this procedure from milk and colostrum samples. Negative controls were also created in the metagenomic dataset; however, a library was not created for sequencing due to low microbial content and negative control samples could not be compared. 16S rRNA data were then filtered using the R labdsv package (Robeson et al., 2021) to remove ASVs that were likely sequencing artifacts as these ASVs were present at very low frequencies (n < 5) or in only 3 or fewer samples.
The R vegan package (Oksanen et al., 2022) was used to perform alpha diversity analyses, beta diversity Bray-Curtis distances, and PERMANOVA calculations. The R ape package (Paradis et al., 2023) was utilized to conduct Principal-coordinate analyses based on Bray-Curtis distances among dietary treatments. Discriminant taxa for each treatment was identified using the labdsv package (Roberts, 2023) in R using a threshold of indicator values of > 0.5. Indicator values range 0 to 1 and represent taxon mean abundances and frequencies of occurrence as an indicator value of 1 indicates a genus is present in all samples of a treatment and occurs in high mean abundances (Legendre and Legendre, 2012). Heatmaps showcasing discriminant indicator species or gene families were created using the R package pheatmap. Piglets selected were chosen based on their dam’s dietary treatment and birth weight, thus a generalized linear mixed model (GLM) was fitted to alpha diversity matrices, the first axis scores of each generated Principal Coordinate Analysis plot, and each taxon, gene family, and pathway with an indicator value of 0.50. The GLM model used for sow samples considered dietary treatment as a fixed effect and sow farrowing group as a random effect. The GLM model used for piglet fecal samples considered dam dietary treatment and birth weight as fixed effects and sow farrowing group as a random effect. Additionally, the MaAsLin2 R package (Mallick et al., 2021) was utilized to determine multivariate associations between microbial features and grouping variables that occurred for piglets and sows. The MaAsLin2 model for sows (colostrum, milk, vaginal, and fecal samples) included the fixed effect of dietary treatment with the random effect of farrowing group. Within sow samples, no significant associations were identified regardless of sample type. The MaAsLin2 model utilized for piglet fecal samples explored four different models including no random effects, the random effect of sow farrowing group, the random effect of sow, and the combination of group and sow as random effects. No significant associations were reported using the piglet MaAsLin2 model with random effects including sow or sow and group combined. Individual box plots were created using base R functions or using MaAsLin2. Testing statistical significance was performed using Kruskal-Wallis tests, Wilcoxon tests, MaAsLin2 significant associations, or PERMANOVAs for all nonparametic microbiome data. Discriminant taxa, gene families, and pathways identified using the MaAsLin2 package with an FDR adjusted P-value (q-value) < 0.10 were considered. All data were considered statistically significant at P ≤ 0.05 and marginally significant at 0.05 < P ≤ 0.10. Statistical significance in all figures is denoted with three asterisks if the P-value < 0.001, two asterisks if the P-value < 0.01, one asterisk to denote a P-value < 0.05, and a cross when the P-value is < 0.10.