and probability values. Multiple t-tests were conducted post-hoc in the emmeans package (version 1.4.8; Lenth, Singmann, Love,  Buerkner & Herve, 2018) in 'R' (v4.0.2) to test for pairwise hybrid-cross differences, and corrections for multiple tests were made using Bejanmini & Hochberg adjusted p-values (Benjamini & Hochberg, 1995).

To explore differences in the microbiome community composition, the rarefied ASV table was used to measure pairwise distances using an abundance-based metric (Bray-Curtis), and a presence/absence-based metric (Jaccard). Patterns of differences among crosses were visualized by plotting 2D principal coordinate analyses (PCoA) using the two most variance- explaining eigenvectors for Bray-Curtis and Jaccard distance matrices. To simplify visualization, principal coordinate scores were extracted using PAST (v3.25; Hammer, Harper & Ryan 2001), averaged for each cross and plotted as a centroid, and then, 95% confidence interval were computed to represent dispersion error bars. To visualize ASV presence/absence patterns among crosses, we listed ASVs present in the rarefied table for each cross and then created an intersection plot in the UpSetR package (v1.4.0; Lex, Gehenborh, Strobelt, Vuillemot & Pfister, 2014) in R (v4.2.0; R Core Team, 2016). The intersection plot shows the ASVs found exclusively in certain stocks, and ASVs found commonly among multiple stocks.

To measure mean differences in community composition, permutational analysis of variance models (PERMANOVAs; Anderson, 2001) tests were conducted using an 'R' vegan package wrapper in QIIME2. First, global nested-PERMANOVAs were run to test for differences in pairwise distance matrices (Bray-Curtis & Jaccard) among crosses, among sires (nested in crosses) and between replicate pens (nested in sires, with sires nested in crosses). To further explore specific patterns in community composition differences in the microbiome among crosses, post-hoc pairwise PERMANOVAs were conducted for pairs of crosses to determine hybrid cross-specific differences. Pairwise PERMANOVA tests were corrected for multiple comparisons among pairs of crosses using false rate discovery (Benjamini & Hochberg, 1995). Models were also run for each cross to determine the statistical significance of microbiome compositional differences due to sire and pen effects (nested within sire) for each cross separately. To determine whether significant statistical differences in crosses were attributed to differences in locational (mean composition) or dispersion (spread) differences, ad-hoc permutational homogeneity of dispersion (PERMDISP; Anderson, 2004) tests were conducted globally and pairwise among crosses. Each non-parametric model was run with 9,999 permutations, and pairwise tests were corrected for multiple comparisons among pairs of crosses using false rate discovery (Benjamini & Hochberg, 1995).

To determine the presence and the nature of taxa driving community composition differences among crosses, linear discriminant analysis effect size (LEfSE; Segata et al., 2011) tests were run in MicrobiomeAnalyst (Dhariwal et al., 2017) on relative abundances of ASVs. To prevent rare ASV effects, we filtered ASVs occurring in 5 samples or less, and then removed ASVs with 200 reads or less across all remaining samples. The remaining ASVs corresponded to 96.4% of all reads in our filtered ASV table.