Total nr animals produced Table 2 ). Sequence tags were rarefied to ~20,000 per sample, and clustered in 32,032 OTUs 137 (97% similarity threshold). 12,054 OTUs present in at least 0.2% of the samples (with more 138 than two tags in at least two samples) and amounting to an average of 98.7% of sample reads, 139
were retained for further analysis. They were annotated to 41 phyla, 87 classes, 149 orders, 140 207 families, 360 genera and 150 species. The number of OTUs detected per sample averaged 141 1,575 (range: 54 to 3,112) (Suppl. Table 2 ). The first two Principal Coordinates (based on of 58 OTUs that were annotated to 21 taxa were identified in >95% of day 120 and 240 feces 155 and cecum content samples of both F6 and F7 generations, hence defined as core bacterial 156 taxa (Suppl. Fig. 1A ). -diversity (measured by Shannon's index) of fecal samples was lower 157 at day 25 than at days 120-240, reminiscent of the enrichment of the intestinal flora observed 158 between child-and adult-hood in humans (Yatsunenko et al., 2012; Radjabzadeh et al., 2020) . 159
It was also lower for ileal content than for cecal content and mucosa (and probably ileal 160 mucosa) (Fig. 2C ). Six percent of ileal content samples harbored less than 100 OTUs (Suppl. 161 Table 2 ).
-diversity (measured by pair-wise Bray-Curtis dissimilarities) tended to be 162 which has largest n). The combined p-value for the five traits combined and computed as 234 above was 0.013, hence providing a second line of evidence for an effect of genetics on 235 microbiome composition in this population (Fig. 3B) . 236
We then evaluated the heritability (ℎ 2 ) of the abundances of individual taxa using a mixed 237 model. This was done for up to 29 phyla, 53 classes, 86 orders, 116 families, 148 genera and 238 4,240 OTUs per data series. Heritabilities were estimated using a mixed model implemented 239
with GEMMA (Zhou & Stephens, 2012) . The model included random polygenic and residual 240 error effects. Kinship coefficients (to constrain the polygenic effect) were computed from 241 whole-genome SNP data, also using GEMMA (Zhou & Stephens, 2012) . To obtain unbiased 242 ℎ 2 estimates and associated p-values, we repeated the analysis 1,000 times with abundances 243 randomly permuted within litter. The average ℎ 2 across the permutations was then 244 subtracted from the ℎ 2 obtained with the unpermuted data to yield conservative estimates 245 of ℎ 2 . Their p-values were estimated as the proportion of permutations yielding an equally 246 high or higher ℎ 2 estimate. Analyses were conducted for the 12 measured data series. P-247 values were ≤ 0.05 for 4,219 (=14%) of the 30,127 realized tests, hence above random 248 expectations (Suppl. Table 4 ). The correlation between F6 and F7 ℎ 2 estimates (or their 249 log(1/p) values) were positive and highly significant (p ≤ 1.1x10 -30 and 1.07x10 -17 , respectively) 250 for D240 fecal samples, hence supporting genuine genetic effects at least for this trait (Suppl. 251 There was no significant correlation between pig and human values when considering all taxa. 257
It is noteworthy, however, that the taxon found to be the most heritable in human, the family 258
Christensenellaceae (Goodrich et al., 2014&2016) , was also amongst the most heritable in 259 pigs. At least seven other taxa were found to be heritable in both human and porcine adult 260 feces (Fig. 3D) . Table 5 ). Both models were implemented with the GenABEL R package (Aulchenko et al., 2007 ) 300 and included sex, batch and the three first genomic principal components as fixed covariates. 301 P-values were further adjusted for residual stratification by genomic control. We obtained 302 1,527 signals encompassing at least three variants with value ≤ 5 × 10 −8 (the standard 303 genome-wide significance threshold). To evaluate whether this number exceeded 304 expectations assuming that all tests were null hypotheses, we performed two analyses. In 305 the first we chose one of the largest (hence best powered) data series (day 240 feces in the 306 F7 generation) and repeated all GWAS on a dataset with permuted (within litter) genotype 307 vectors. The number of microbiota QTL (mQTL) signals detected with the real dataset was 308 221, while the number detected with the permuted dataset was 152, hence suggesting a true 309 discovery rate of ~30%. In the second, we collected -for each of the 1,527 signals with p-310 value(s) ≤ 5 × 10 −8 described above (corresponding each to a lead SNP x taxon x trait x 311 cohort x model combination) -the p-values for the same SNP x taxon x model combination, 312 yet for all other trait x other cohort combinations. Thus, we would typically collect ~5-7 p-313 values for each such signal. We reasoned that if the initial signals included a sufficient 314
proportion of true positives, the lead SNPs would have similar effects in at least some of the 315 traits in the other cohort and the collected p-values concomitantly shifted to low values. The 316
corresponding distribution of p-values was examined by means of a QQ-plot, and was 317 compared with the distribution obtained with an equivalent number of randomly selected mQTL in our data (Fig. 4A ). Of note, the average (F6 and F7) number of genome-wide 322 significant mQTL was positively correlated with the average (F6 and F7) taxon's heritability, 323 particularly for D240 fecal samples (p = 5.2x10 -6 ) (Suppl. Fig. 3B ). 324
To identify the corresponding mQTL yet properly accounting for the large number of realized 325 tests before declaring experiment-wide significance and simultaneously provide confirmation 326
in an independent cohort, we performed meta-analyses (across traits) separately in the F6 327 and F7 generations for the 1,527 above-mentioned signals. We designed an empirical meta-328 analysis approach that accounts for phenotypic correlation across traits if it exists (cfr. 329
Methods). The discovery threshold was set at 0.05 10 6 ×1,527×2×2 = 8.2 × 10 −12 hence corrected 330 for the size of the genome, the number of tested taxa (because genome-wide significant), the 331 two used statistical models, and the two studied cohorts (F6 and F7). There was no need to 332 correct for the number of traits as the meta-analysis generated one statistic for all traits. The 333 confirmation threshold was set at 0.05/ where was the number of signals exceeding the 334 discovery threshold in at least one cohort. Thus, we searched for signals (defined as lead SNP 335
x taxon x method combinations) that would exceed the discovery threshold in either the F6 336 or F7 cohort and the confirmation threshold in the other. We identified seven signals 337 exceeding the discovery threshold in at least one cohort, hence setting . For six of those, as well as genus p-75-a5 to which OTU-476 is assigned (Suppl. standing out as being highly significant in both F6 and F7 (Fig. 4E ) All lead SNPs of the meta-analyses conducted in the F6 and F7 generation map to the 3' end 398 of the porcine acetyl-galactosaminyl transferase gene that is orthologous to the gene 399 underlying the ABO blood group in human (Fig. 4B ). This is a strong candidate gene known to 400 modulate interactions with several pathogens (Cooling et al., 2015) , although not directly 401 known to affect intestinal microbiota composition in healthy humans (Davenport et al., 2016) . 402
Four non-synonymous ABO SNPs (R37G, A48P, S60R, G66C) segregated in the F6 and F7 403 generation but none of these were in high LD with any of the lead SNPs ( 2 ≤ 0.18 Fig. 4D ). In these samples (D120, D240, CC, CM), AO genotype explained 440 on average 7.9%, 3.2% and 6.6% of the variance in abundance for OTU-476, OTU-327 and 441 genus p-75-a5 (Suppl. Fig. 4E ). Of note, the abundance of OTU-476 and OTU-327 was shown 442 to be highest in cecal content where they account on average for respectively ~0.92% and 443
~0.02% of reads in AA/AO animals, and for 0.47% and 0.003% of reads in OO animals ( Fig. 5  444 and Suppl. Fig. 4F ). 445
The ABO locus is known in humans to be under strong balancing selection that has GalNAc ( 1-3 linkage) to a variety of glycan substrates sharing a Fuc 1-2Gal 1-4GlcNAc or 515
Fuc 1-2Gal 1-3GlcNAc (H antigen) extremity (Cooling, 2015) . In the gut, these include the 516 heavily glycosylated secreted and transmembrane mucins constituting the cecal mucus. repressors (agaR, gntR), for a total of six essential pathway constituents (Suppl. Table 7) . 548
Genes involved in the utilization of specific sugars (including GalNAc) tend to cluster and form 549 operons of potentially coregulated genes (regulons) that support all or most of the essential 550 TR/CP steps. The steps that are not encoded by the operon may be complemented in trans 551 by genes encoding enzymes that are often less substrate-specific (Lawrence, 1999; Koonin, 552 2009). We used GhostKOALA (Kanehisa et al., 2016) to search for orthologues of the 24 genes 553 in the two OTU476-like genomes and 3,111 MAGs. We generated two scores to evaluate the 554 capacity of bacterial species to utilize GalNAc. The first (pathway score) counted the number 555 of essential steps in GalNAc utilization (out of six) that could be accomplished by the set of 556 orthologues detected in the genome (cfr. Methods), irrespective of their map position. The 557 second (regulon score) counted the number of essential GalNAc utilization steps that could 558 be fulfilled by orthologues that were clustered in the genome, i.e. forming a potential operon. 559
Following Ravcheev and Tiele (2017), we used agaS as anchor gene to establish the regulon 560 score, i.e. we counted how many essential steps in GalNAc utilization (out of six) were covered 561 by genes located in the vicinity of agaS. 562
The first striking observation was that at least one orthologue of agaS was found in the two 563 (=100%) OTU476-like strains (4-15-1 and 4-8-110), in 31% of Erysipelotrichaceae MAGs 564 (n=248), yet in only 3.0% of other MAGs (n=2,863). The second, was that both scores 565 (pathway and regulon score) were very significantly higher for Erysipelotrichaceae than for for Erysipelotrichaceae and non-Erysipelotrichaceae MAGs combined (ppathway=2.2e-3 and 568 pregulon=1.2e-5) (Fig. 6B ). These comparisons accounted for variation in the MAGs' completion 569 score, number of contigs and predicted size of the corresponding genomes (see Methods) . agaR, respectively (Fig. 6C) . 581
We further showed that adding GalNAc in the culture medium indeed enhances the growth 582 of the two isolated OTU-476 like strains, indicating that these can indeed utilize GalNAc as 583 carbon source (Fig. 6D) . and F7 generation is remarkably uniform and close to expectations (i.e. 12.5%) both at 679 genome-wide and chromosome-wide level (Fig. 1C) , suggesting comparable levels of genetic 680 diversity across the entire genome. This does not preclude that more granular examination the same traits across the F6 and F7 generation, yet marked compositional differences 685 between traits (Fig. 2B ). Even at family-level, some taxa are found to be nearly trait-specific Pasteurellaceae, the firmicutes Clostridiaceae, Peptostreptococcaceae, Bacillaceae, 688
Leuconostocaceae, and the actinobacteria Microbacteriaceae are at least ten times more 689 abundant in ileal than in any other sample type. Amongst those, Leuconostocaceae are nearly 690 digesta-specific, while Pseudomonadaceae are nearly mucosa-specific. The Bacteroidetes 691
Odoribacteraceae and Rikenellaceae were found to be at least ten times more abundant in 692 day 25 feces than in any other sample type. The firmicutes Christensenellaceae were nearly 693 ten times more abundant in feces (irrespective of age) than in any other sample type. This To evaluate the importance of host genetics in determining gut microbiota composition we 698 first examined the relationship between genetic relatedness and microbiota dissimilarity. It 699 is worth re-emphasizing that food and environment was very standardized in this experiment 700 when compared to typical human studies. Genetic relatedness between individuals was 701 measured using genome-wide SNP information while microbiota dissimilarity was measured 702 using Bray-Curtis distance (f.i. Rothschild et al., 2017). We relied on two approaches to 703 mitigate confounding of genetic and environmental effects. In the first we restricted the 704 analyses to full-sibs raised in the same environment, i.e. we confronted genetic similarity and 705 microbiota dissimilarity of litter-mates. In the second we confronted genetic similarity and 706 microbiota dissimilarity across generations (F6 and F7), yet avoiding parent-offspring pairs. 707
Both approaches supported an effect of genetics on microbiota composition manifested by 708 significant negative correlations between genetic similarity and microbiota dissimilarity for 709 (some) individual traits as well as when combining information across traits (Fig. 3A&B) . 710
Regressing squared phenotypic difference on genetic distance is an established way to 711 estimate local and global heritability (Haseman & Elston, 1972; Visscher et al., 2006 ). Yet, 712
Bray-Curtis distance is peculiar in that the phenotypes between which a "difference" is 713 measured are not defined per se (Bray and Curtis, 1957) . To nevertheless evaluate to what 714 degree of heritability the observed negative correlations might correspond, we simulated quantitative traits with various degrees of heritability in the actual pedigrees and examined 716 the distribution of ensuing correlations between phenotypic distance (absolute value) and 717 genetic distance. These analyses indicated that (in the studied, genetically highly divergent, 718 population) the heritability of microbiota composition may be of the order of ~0.80 within 719 litter, and ~0.20 in the overall population (Suppl. Fig. 7) . That the heritability is higher within 720 litter than in the overall population is expected as the environment is obviously more 721 homogeneous within than across litters and generations. Strikingly the impact of genetics was 722 strongest for fecal samples at day 240 in all analyses. This may be in agreement with the 723 observation that, in human, microbiota composition stabilizes with age (Aleman & Valenzano, 724 2019). Yet, why heritability should be higher in feces than for ileal and cecal content and 725 mucosa remains unclear. Sample types with higher alpha-diversity may be more resilient and 726 hence more heritable. 727
We also measured the heritability of the abundance of individual taxa. As before, we only 728 extracted within-litter information to mitigate confounding between environment and 729 genetics. Convincing evidence for a genuine influence of host genetics on taxa abundance 730 was the observation of a significant correlation between heritability estimates in the F6 and 731 F7 generation for fecal samples at day 240 and -to a lesser extend -cecum content (Suppl. 732 Larger windows will increasingly be contaminated with recombinant A-O 766 haplotypes blurring the sought signal. Indeed, for windows ≥ 20-Kb or more, the gene tree 767 corresponds to the species tree, while for windows ≤ 15-Kb the tree sorts animals by AA vs 768 OO genotype (Suppl. Fig. 8 Fig. 8 ). This suggests that the O 772 allele may be older than the divergence of the Phacochoerus and Sus A alleles, i.e. > 10 MYA. 773
It will be interesting to study larger numbers of warthog to see whether the same 2. behavior. No significant effects were observed when accounting for multiple testing (Suppl. 792 Fig. 9 ), including those pertaining to immunity and disease resistance. 793
It is noteworthy that the old age of the "O" allele must have contributed to the remarkable 794 mapping resolution (≤ 3 Kb) that was achieved in this study. In total, 42 variants were in 795 near perfect LD ( 2 ≥ 0.9) with the 2.3 Kb deletion in the F0 generation, spanning 2,298 bp 796
(1,522 on the proximal side, and 762 on the distal side of the 2.3 Kb deletion). This 2.3 Kb 797 span is lower than genome-wide expectations (17th percentile), presumably due to the 798 numerous cross-overs that have accrued since the birth of the 2.3 Kb deletion that occurred 799 in the distant past (Fig. 1E ). Yet the number of informative variants within this small segment 800 is higher than genome-wide average of (57% percentile) also probably due at least in part to 801 the accumulation of numerous mutations since the remote time of coalescence of the A and 802 O alleles (Fig. 1D) . 803
The chromosome 1 QTL was the only signal that exceeded experiment-wide discovery and 804 conformation thresholds. QQ-plots obtained after removing chromosome 1 variants (272.8-805
273.1Mb interval) did not show convincing evidence for residual inflation of log(1/p) values 806 (Suppl. Fig. 3D ). This suggests that the residual heritability most likely has a highly polygenic 807 architecture, as becoming increasingly apparent for most complex traits. and 327) and genus p-75-a5, but affected at least some other Erysipelotrichaceae as well. As 811 mentioned above, effects of ABO genotype on host-pathogen interactions are usually 812 interpreted in the context of adhesion or immune response. Yet an alternative mechanism 813 is by altering the source of carbon upon which intestinal bacteria feed. Small and large 814 intestine are amongst the tissues in which the ABO gene is the most strongly expressed (Suppl. 815 least two predictions. The first is that the intestinal GalNAc content should be higher in AA 822 than in OO pigs and this was indeed shown to be the case (Fig. 6A) . The second is that the 823 bacteria affected by the mQTL should be able to use GalNAc. We isolated two strains with 824 16S rRNA sequences that were near-identical to those of the OTU strain (OTU476) that was 825 most affected by the mQTL, and sequenced their complete genome. We showed that it 826 contained the orthologues of eight genes known to be essential for GalNAc import (AgaPTS: 827 agaE, agaF, agaV, agaW) and catalysis of the first four GalNAc-specific degradation steps 828 (deacetylation: nagA; demanination/isomerisation: agaS; kination: fruK; aldolase: gatZ) 829 hence the five key steps in GalNAc utilization. Importantly, the eight genes clustered in a 50 830
Kb chromosome segment (Fig. 6C) . We generated 3,111 porcine intestinal MAGs from 831 metagenomic shotgun data for comparison. None of these would harbor a GalNAc gene 832 cluster encoding more than four of the five key steps. One catalytic function was always 833 provided in trans, whether GalNAc-6-P deacetylase, tagatose-6-P kinase, or tagatose-1,6-PP 834 aldolase. This finding clearly revealed the unique status of OTU476 with regards to GalNAc 835 utilization. Also consistent with the QTL findings, Erysipelotrichaceae MAGs were strongly 836 enriched in clustered GalNAc TR/CP orthologues when compared to MAGs assigned to other 837 bacterial species. Finally, the growth of the two isolated OTU476-like strains was shown to 838 increase when fed with increasing concentrations of GalNAc (Fig. 6D) . 839
Amongst the 3,111 studied MAGs, 15 harbored a gene cluster able to sustain four of the five 840 steps, while the fifth enzyme was encoded somewhere else in the genome. These included unidentified Firmicutes (Fig. 6C) . Gene order within the corresponding GalNAc clusters was 844 supported by observation in two or more independent MAGs and/or by the fact that all 845 concerned genes resided on the same contig. For all 15, the orthologue needed to fulfill the 846 fifth enzymatic reaction (2x tagatose-1-P kinase, 2x GalNAc deacetylase, 1x tagatose-1,6-PP 847 aldolase) was found somewhere else in the genome, allowing us to assume that all these 848 species are able to utilize GalNAc. Fifty additional MAGs contained the orthologues needed 849 to accomplish the five key steps in GalNAc utilization albeit without evidence for a similar 850 degree of clustering (either because the genes are indeed not clustered in the corresponding 851 genomes or because they were segregated across distinct sequence contigs). It is reasonable 852 to assume that several of those bacteria are also able to utilize GalNAc as carbon source. Why 853 then would the chromosome 1 mQTL only affect a small subset of Erysipelotrichaceae species? 854
We first reasoned that OTU476-like strains might be more dependent on GalNAc availability 855 than other species, for instance because they can't utilize alternative, common 856 monosaccharides as carbon source. To test this hypothesis, we searched for KEGG numbers 857 that would commonly occur in other genomes, yet were absent in the OTU476-like strains. 858
We performed this analysis for all MAGs, as well as separately for the MAGs that were 859 predicted to be able to use GalNAc (cfr. above). There was no convincing evidence that 860 By regulating expression of their GalNAc operon in response to ambient GalNAc availability, 894 species like E. Coli may fair equally well in the gut of AA/AO as in that of OO pigs, hence not 895 be affected by the mQTL. It is worth noting that the A allele is dominant with regards to 896 OTU476, OTU327 and p-75-a5 abundance (Fig. 5B) , suggesting that the additional increase in 897
GalNAc concentrations in AA (vs AO) animals does not further benefit these taxa. 898
We examined the effect of ABO blood group on the abundance of ~75 OTUs assigned to 899
Erysipelotrichaceae in human gut samples. Although we could not rigorously test this for all 900
OTUs (as some human V1-V2 and V5-V6 data could not directly be compared with porcine V3-901 V4 data) none of the OTUs detected in human samples were as closely related to the pig OTU-902 476, OTU-327 or p75.a5 as these were to each other. We found no evidence for an effect of 903 ABO blood group on the abundance of any of these OTUs. What underlies the difference 904 between pigs and humans is unclear. Either strains susceptible to ABO genotype are not 905 present at sufficient frequency in human feces, or the carbohydrate composition of human intestinal content makes these strains less sensitive to variations in GalNAc concentrations. 907
It is worth noting that the studied human samples were intestinal biopsies collected after a 908 standard gut cleansing procedure. The abundance of the genus p-75-a5 was recently found 909 to differ significantly between African subsistence categories and to be highest in pastoralists 910 similarity (from SNP genotype data) and microbiota dissimilarity (Bray Curtis distance 918 computed from 16S rRNA data) both within litter as well as between generations, supporting 919 an effect of host genetics and intestinal microbiota composition (Fig. 3A&B) . We took care in 920 these analyses to mitigate effects of litter on both genetic and microbiota distance metrics, 921 (as these may inflate statistical significance) by applying permutations tests. Regressing 922 squared phenotypic difference on genetic distance is a standard way to estimate local or 923 global heritability (Haseman & Elston, 1972; Visscher et al., 2006) . It can be shown that 924 −̂/(2̂2) estimates the narrow sense heritability ĥ 2 . In these, ̂ is the least square 925 regression coefficient and ̂2 an estimate of the phenotypic variance. In our analyses, and (microbiota composition) was deduced from the coincidence between the red and green lines. 947
As an example, the heritability of microbiome composition for data series F7-D120 was 948 assumed to be close to 0.8. We proceeded in the same way for the across generation analysis 949 (Suppl. Fig. 7C ). We used the actual measures of kinship across the F6 and F7 generations 950 computed with GEMMA. We standardized them (mean 0 and SD 1), and then scaled them 951 such that the values for an individual with itself would center on 1, and for full-sibs on 0.5. 952
Breeding values and environmental effects were sampled using mvrnorm and rnorm as above. Both models included sex, slaughter batch (21 for F6, 23 for F7) and the three first genomic 1155 principal components as fixed covariates. GWAS were conducted separately for each taxon x 1156 data series combination and p-values concomitantly adjusted for residual stratification by 1157 genomic control. P-values were combined across traits and/or taxa using a z-score. P-values 1158 were converted to signed z-values using the inverse of the standard normal distribution and 1159 summed to give a "z-score". Z-scores were initially calculated using METAL (Willer et al., correlation that exists between the phenotypic values of a given cohort across traits we also computed the genome-wide (i.e. across all tested SNPs) average (Z ̅ ) and standard deviation (σ Z ) specific annealing temperature for each set of primers and 72°C for 2 min), and 72°C for 10 1194 minutes on a PE 9700 thermal cycler (Applied Biosystem, USA). bins (metagenomic assembly genomes, MAGs) generated by the two binning algorithms were 1319 evaluated for quality and combined to form a MAG set using the bin_refinement module in estimate the completeness and contamination of each MAG (Parks et al., 2015) . The MAGs with completeness <50% and contamination >5% were filtered out. Non-redundant MAGs were then clustered to 97% identity level OTUs using the DNACLUST program (Ghodsi et al., from cecum samples, and contributed to ABO genotyping. YZ participated in the preparation