Influence of genetic background and dietary oleic acid on gut microbiota composition in Duroc and Iberian pigs

Background Phenotypic variability for productive and meat quality traits has been largely studied in Iberian pigs, especially in genetic selection and nutritional experiments. Complex interactions among genetic background, diet composition and gut microbiota hinder the correct assessment of each factor’s contribution on phenotypes. In order to disentangle these interactions, we evaluated changes in gut microbiota composition comparing 48 Iberian and Duroc pigs fed diets with different energy source (standard diet with carbohydrates vs sunflower oil-enriched diet with high oleic acid content). Results A higher richness was observed for Iberian pigs (p < 0.05) and compositional analysis was applied for beta-diversity, differential abundance and pairwise log-ratio analyses. We found significant differences in overall microbiota composition between breeds, and also between diets inside breeds, to a lesser extent. Differential abundance analysis revealed that Duroc animals have more proportion of Actinobacteria and Prevotella, while Iberian replace those microorganisms with other more variable taxa. According to dietary differences, high-oleic fed animals were richer in Prevotella. We also found microbial ratios capable of separating animals by breeds and diets, mostly related to Actinobacteria. Conclusion This study reveals that both genetic background and diet composition might have a relevant impact in gut microbiota composition. The application of compositional data analysis has facilitated the identification of microorganisms and ratios as possibly related to metabolic changes due to genetic background and, to a lower extent, to dietary changes. This may lead to a relevant progress in the knowledge of interactions between pig genetics, environment and gut microbiota.

resources available 2 . The combination of their natural high tendency to fat accumulation and the acorn and grass feeding in the late fattening phase ("Montanera"), confers their final products a unique fatty acid (FA) profile composition and a 11 higher intramuscular fat (IMF) content than commercial breeds. 12 The interest in this high-quality market has led to the search and legal regulation of strategies for the improvement of its 13 efficiency and productivity. In this sense, one of the most common practices nowadays is crossbreeding with Duroc, 14 resulting in a negative impact on final product composition and quality 3 . Crossbreeding consequences have been largely 15 studied through different points of viewsuch as phenotypic traits, meat quality, genetic markers or gene expression 4-6 . 16 As the high quality of Iberian products has become the flag for the Iberian pig market to enter international markets, 17 producers keep an interest in improving their products' organoleptic attributes through FA profile and IMF content 18 optimization 7,8 , especially in the recent more intensive production systems. Several strategies such as genetic selection, in Iberian pig genotypes subjected to different diets 4,9-11 . 23 Besides, with the advent of new mass-sequencing techniques, alternative ways have been opened for animal breeding and 24 for in-depth studies of the molecular basis of phenotypic variability. Microbiota studies have hugely benefitted from these 25 NGS techniques, as culture-free sequencing allows obtaining a lot more information from microbiome communities. Gut 26 microbiota has become a key focus in animal breeding industry, as well as in human health, since gut microbiota and host 27 are linked in a bidirectional way; first, diet composition and other environmental elements, as well as animal genetics, 28 might modify microbial population structure 12 ; second, microbiota has been reported to have an impact on animal 29 physiology, on multiple animal productive traits and also on the final products' quality and composition 13,14 . This intricate 30 correlation network makes gut microbiota composition a key factor for understanding the relationships between genotype, 31 phenotype and environment. As a consequence, it might result obvious that different genetic backgrounds can promote 32 different microbiota population structures as it has been proven in other studies 15 . For this reason, the two main breeds 33 employed in Iberian pig production (Iberian and Duroc) might have differences in gut microbiota composition possibly 34 affecting relevant phenotypic traits. Also, dietary modification or supplementation strategies such as FA addition might 35 have an impact on gut microbiota composition, potentially contributing to differences in animal traits, with special 36 relevance on final product composition and quality. Understanding all these relationships and their implications is crucial 37 to comprehend their global effect on animal phenotype.

38
In this study we compared the gut microbiota composition in pigs from Iberian and Duroc breeds, fed two different 39 isocaloric and isoproteic diets with different energy source (oleic acid vs carbohydrates), in order to dissect the diet and genetic background effects on the gut microbiota composition, as well as to deepen in the metabolic differences between 41 breeds and diet groups, potentially mediated by gut microbiota. 42

44
Experimental design and sampling 45 The experiment was carried out at the facilities of the Instituto Tecnológico Agrario de Castilla y Leon (ITACYL) Pig 46 Test Center (Hontalbilla, Segovia, Spain). A total of forty-eight castrated male purebred pigs, 19 Duroc (DU) and 29 Iberian (Torbiscal strain) 16 (IB), were raised in the same commercial farm, weaned at 28 days old and transferred to the 48 experimental facility one month after weaning. The animals were kept under identical management conditions, housed in 49 batches of 4 pigs/pen (1 m 2 pig -1 ), with a concrete floor and straw bedding. Temperature was controlled at a mean of 23  Raw reads were pre-processed using Prinseq-lite tool 18 and in-house software. After trimming low quality bases in each 74 read, sequences shorter than 50 bp and with an average quality score below 30 on a window of 20 bases were discarded.

75
After joining forward and reverse files, a total amount of 11,599,350 reads were processed using Quantitative Insights 76 Into Microbial Ecology (QIIME) version 1.9.1 19 into operational taxonomic units (OTUs), following the de novo OTU-77 picking approach, with a clustering threshold of 97% identity. Chimeras were filtered using USEARCH algorithm (v. 6.1) 78 20 . SILVA reference database (release 132) 21 was used for taxonomic assignment and chimera filtering. An OTU 79 abundance threshold of 0.005% was established for final quality-filtering, as described in Bokulich et al. 22  transformation. This process moves the data to an additive scale, which makes the performance of multivariate hypothesis 101 testing easier. The analysis was repeated glomming the dataset to OTU, genus, family and phylum levels. R packages 102 zCompositions 26 and easyCODA 27 were used for this transformation.
Beta-diversity: Beta-diversity was computed using a compositional Principal Component Analysis (i.e., PCA after CLR transformation of the data). For a better visualization of the behaviour of both breed and diet effects, a mixed "Breed-105 Diet" factor was considered as PCA grouping variable. The distance matrix was built using Aitchison distance measure, 106 which is the Euclidean distance after CLR transformation. Individual beta-diversities were calculated as the distance to 107 centroid for groups from mixed "Breed-Diet" effect (using vegan::betadisper). Homogeneity of dispersions respect to the 108 group centroids was measured through ANOVA-like permutation test (vegan::permutest) and post-hoc Tukey HSD test.

109
For testing the significance of multivariate effects in our a priori groups, non-parametric permutational multivariate 110 analysis of variance (PERMANOVA) 28,29 was performed using adonis2 function from vegan package in R 30  118 normalization by estimation of size factors as the median ratio of counts 33 and negative binomial GLM fitting with Wald 119 significance tests were performed. Differences between breeds and between diets were analysed in separate models 120 correcting by the other factor. Interactions between factors were also evaluated with an additional model ( Table 1). P-121 values were corrected for multiple testing through Benjamini-Hochberg method. Significantly differentially abundant 122 (DA) OTUs were considered with a false discovery rate (FDR) cut-off of 0.05 and a fold-change (FC) higher than 1.5 or 123 lower than -1.5 (i.e., |log2FC| > 0.59).

126
Input for ALDEx2 analysis was prepared as the CLR-transformed posterior distribution of the data generated by 128 Pairwise log-ratios analysis 135 Closured data were also used for another approach, using pairwise log-ratios (p-LR) instead of centred log-ratios 27

170
The overall composition is summarised in Table 2. Our final OTU table was composed by 47 samples, as we pruned one 171 sample with less than 30k reads. As mentioned in methods, we applied an additional filter removing those reads not 172 assigned to Bacteria kingdom. After this step, 97.44% of total reads were kept, representing 1,246 OTUs. From this 173 bacteria-classified fraction, 99.93% reads were assigned to phylum level. A total of 17 phyla were identified. Prevalence 174 and abundance of phyla are represented in Figure 1. Total OTU number, % of sequences assigned to a known taxon at different ranks and unique taxa found at each rank. 1 Proportion of total reads classified to phylum and lower taxonomic levels are calculated relative to the total reads after filtering OTUs not assigned to Bacteria kingdom.  (i.e., genus and OTU levels). Centroid position also revealed that separation between breeds is higher than separation between diets. At OTU level, centroid distance is maximum for both breed and diet groups. Individual beta-diversities 195 were calculated as the Euclidean distance to centroids for groups from mixed "Breed-Diet" effect. Average individual 196 dispersion between groups was not significantly different, except for OTU-glommed dataset. At OTU level, DU-C 197 samples had an average beta-diversity index significantly lower compared to DU-O, IB-C and IB-O groups, although p-198 value in DU-C vs IB-C contrast (same diet in different breeds) was close to 0.05. (Figure 4).

199
PERMANOVA analysis confirmed that samples from different breeds are significantly different. Differences between 200 diets (corrected by breed effect) are also significant at OTU, genus and family level, although these differences blur at 201 phylum level (Table 3). Nonetheless, the differences between groups are moderate, as roughly 10% of variance from our 202 data can be explained by Breed effect, and only 5% by Diet effect. 203

241
At phylum level diet differences were not clear inside each breed, as it was observed in beta-diversity analysis.
At family level, three ratios were capable of separating animals in the four Breed-Diet groups with high precision. As subset, denovo283234 (Treponema 2), which was negatively associated to backfat MUFA and oleic acid content, while, a total of 87 individual OTUs were found as significantly associated to phenotypic traits in IB subset (Figure 9). Most of them belonged to Alloprevotella and other Prevotella groups and were associated to multiple traits. Ruminococcaceae through fermentation 36 and their abundance is highly related to high-fiber long-term diets 37,38 . Their presence in fecal 294 microbiota is highly reported in human 39 and other animals such as ruminants, predominating in rumen microbiota [40][41][42] .

295
Multiple studies have shown the dominance of Prevotella in pig gut microbiota 43 as well as their relevance in enterotypes 296 44 or even immune response 45 . Alloprevotella genus is also one of the most representative in our data, with an average 297 RA of 4.2% in DU and 11.6% in IB animals. This genus is also known for its saccharolytic activity 46 , although it is less 298 studied than Prevotella due to its lower abundance in gut communities. Lachnospiraceae and Ruminococcaceae families 299 also comprehend well-known polysaccharide fermenters, which in addition participate in methanogenesis by producing 300 H2 and formate as substrates to archaeal communities 47 . Abundance balance of these three families in gut microbiota, 301 through a number of abundance ratios, might be relevant for host traits related with feed efficiency, metabolism or health.

302
For instance, Prevotella-to-Bacteroides ratio has been related with weight and fat loss in humans 48 , and Firmicutes-to-

303
Bacteroidetes ratio has been repeatedly associated to fiber metabolism or obesity in human populations 37,49 and murine 304 models 50 , although there is controversy about these links 51 . Prevotellaceae, Lachnospiraceae and Ruminococcaceae also 305 compose the core microbiome in rumen, being shared by multiple ruminants in an important abundance 52 . Our data 306 support that this core microbiota might also be shared by other monogastric species along all digestive tract. Similar same IB-C pig population. However, they reported a range of 35 to 41% of Prevotella abundance within all colon regions.
These discrepancies must be attributed to different intestinal region, different genotype and diet of the total animal 310 population studied (colon from IB-C animals vs rectum from IB and DU fed C and O diets

327
Differential abundance analysis gives a more specific insight about these changes in microbiome population structure 328 among experimental groups. DESeq2 is a RNAseq-based, robust and well documented method broadly used in 329 microbiome DA analyses. However, in recent years the concern about compositional nature of these data has increased, 330 and some authors remark that RNAseq-based tools are not accounting for the compositional nature of microbiome data.

331
For this reason we also performed differential abundance analysis with ALDEx2 method, an ANOVA-based approach 332 which accounts for this intrinsic compositionality and has a similar sensitivity than other methods, also reducing false 333 positive rate near to zero 24

361
Bifidobacterium are well known anaerobic lactic acid producer bacteria 58 whose abundance in fecal microbiota has been 362 reported as negatively correlated with pig age, being more abundant at early life stages 59  degradation and deposition rates in Iberian pigs, compared to other commercial breeds 62 . Corynebacterium might be contributing to this differential protein synthesis and degradation balance, through a differential bioavailability of essential amino acids in the gut, thus, their abundance could be correlated to animal protein deposition.
Differential abundance of other groups might be related with richness and resilience of Iberian pigs. In this sense, the 375 lower abundance of multiple Prevotella OTUs in IB pigs seem to be related with the previously reported increase of 376 microbial diversity, as it happened when comparing C and O diets. This richness increment might likewise be related with 377 the lower abundance of Shigella-Escherichia group, as a higher α-diversity has been associated to the prevention in the 378 establishment of pathogenic microbes such as Enterobacteriaceae 37   carried out with the same animals have proven the existence of a differential response to O diet for several productive 405 traits and for adipose tissue transcriptome 11 . In summary, this differential response may exist for gut microbiome 406 composition and structure, but again might be not clearly noticed due to the low number of animals and to the relatively 407 high resemblance between feed formulations and to the intrinsic difficulty in evaluating and interpreting interaction 408 effects.

409
We also used the log-ratio analysis approach, assuming that ratios between OTUs can also be biologically relevant despite 410 of single association with OTU RAs. Thus, composition differences have been explored more deeply by means of p-LR 411 analysis, which detected those microorganism pairwise ratios whose abundances could separate animals by both breed 412 and diet. In order to check which microbes were mainly responsible of separation between groups at each ratio, 413 supplementary differential abundance analyses were made, using ALDEx2 (as this pipeline uses the same methods as for 414 RPART trees building) with phylum, family and genus-grouped database.

415
With this approach we identified some ratios that could be useful as differentiators between Duroc and Iberian phenotypes.

431
This type of approach might be interesting to detect microbial biomarkers, given their importance as tools for genetic 432 diagnosis or diet determination. In the case of Iberian pig industry this might be especially useful in order to address 433 frauds due to crossbreeding and dietary manipulation, which can importantly affect final products' quality 6,7 .
As an additional approach we related the individual OTU abundances with different fat composition traits, in order to deepen in the potential effects of the gut bacteria on host's metabolism. A former study with the same animals reported 438 significant differences in FA composition between C and O diets 67

459
As our experiment has focused on the bacterial subset of fecal microbiota, interactions with other microbial clades, 460 especially eukaryotes, must be studied for a better understanding of the real nature of microbial population changes. On 461 the other hand, metabolite variations between breeds and diets must be analyzed in order to stablish the true relationship 462 between microbiome and host phenotype. Mid-gut microbiota populations might also be interesting to observe dietary Our study reveals that genetic background has an important impact on pig gut microbiota composition, while dietary 469 changes have a samller and more variable effect, which could be highly dependent on fiber content. It also brings light to 470 the complexity of microbial relationships. We report a higher OTU richness in Iberian pigs, which highlights the 471 importance of genetic background (i.e., animal breed) in overall microbiome composition resemblance, and may suggest 472 a potential relationship between gut microbiota diversity and composition and Iberian pigs' resilience. Dietary 473 modifications had a small effect in microbiome composition, much less important than host genotype. Differential 474 abundance revealed a complex picture of microbe relationships, with a considerable effect of the breed and multiple 475 interactions between breed and diet. In spite of this complexity, we identified relevant DA taxa and taxa ratios as 476 potentially associated to the metabolic differences between breeds and to a lower extent for the diet influence. Finally, an 477 approach to the use of OTU pairwise ratios as a predictive tool has been performed, with some interesting results that 478 should be subject of further researches. Figure 1 -Phyla prevalence and abundance in dataset. Each dot represents an OTU, being the x-axis its average relative abundance and the y-axis the proportion of samples in which it is present. "Unclassified_Bacteria" is the denomination given to those OTUs with no phylum assignment.    contrasts. Red dots represent OTUs with an adjusted p-value below the FDR cut-off (0.05) and a FC value above 1.5 or below -1.5 for DESeq2 contrasts, or a median difference between CLR values above 1 or below -1 for ALDEx2 contrasts. OTUs with Log2(FC) < 0 or diff.btw < 0 are more abundant in Duroc and Control groups, respectively.  Normalized counts from DESeq2 algorithm and standard error of the mean (SEM) are represented in y-axis. Figure 8 -RPART trees from pairwise log-ratios at phylum, family, genus and OTU levels. Each decision level indicates the threshold value of one p-LR to classify a sample in one phenotypic group or another. Each node box displays the classification, the probability of each class at that node (i.e. the probability of the class conditioned on the node) and the percentage of observations used at that node. Class probabilities at each node are sorted as: DU-C, DU-O, IB-C and IB-O. Limma regression analysis, when no diet correction was applied. Edge colour: red = positive correlation; blue = negative correlation. Central nodes have more connections (i.e., significantly associated to more phenotypic traits) than peripheral nodes. Numbered nodes represent OTUs classified as genera with high relevance due to their abundance patterns in our dataset.

Supplementary information
Additional file 1 -Diet composition.
Nutrient and fatty acid composition for Control (C) and High-Oleic (O) diets.

Additional file 2 -Relative abundance of genera
Tables containing overall relative abundance of identified genera and relative abundance by Breed-Diet groups. Mean and standard deviation of RA are shown for each OTU, as well as taxonomic classification. RA_Rank column represents the rareness of each OTU (High_RA: RA ≥ 1%; Mid_RA: 1% > RA ≥ 0.1%; Low_RA: RA < 0.1%).

Additional file 3 -Differentially abundant OTUs
Tables containing differentially abundant OTUs for each factor (Breed and Diet) with both pipelines (DESeq2 and ALDEx2) and overlapping of both pipelines. Prevalence represents the proportion of samples in which each OTU is found. For DESeq2, normalized counts per group (Ncounts), fold change (log2FC) and adjusted p-value (padj) are shown. For ALDEx2, relative abundance per group (RA), median difference in CLR values between factor groups (diff.btw) and Welch's test adjusted p-value (we.eBH) are shown. In overlapping tables information of both pipelines is shown.