The animal experiment was conducted in accordance with Dutch law, and the Dutch Central Authority for Scientific Procedures on Animals (CCD) approved the experiment under license number AVD1040020173948. Furthermore, animal management and experimental procedures were approved by the Animal Care and Use Committee of Wageningen University & Research (Wageningen, The Netherlands).
Experimental design and procedures
To investigate the effect of ROM on early life GI development and immune function in pigs, 33 sows (Hypor Libra, Boxmeer, The Netherlands) and their litters (Maxter × Hypor Libra sow) were used at the Swine Research Centre (Trouw Nutrition, Sint Anthonis, The Netherlands). All sows got inseminated from a single boar with the aim of reducing genetic variation in the litters. After parturition, the piglets immediately received an ear tag and an intramuscular iron injection, and their birth weight and sex were recorded. Piglet tails were docked as corrective procedure to prevent tail biting. To lower the risk of infection pressure in the stable, Calcium carbonate powder (Power-Cal®, Power-Cal, Sint-Oedenrode, The Netherlands) was spread to all pen floors. One day after parturition, 192 female piglets were selected and cross-fostered to minimize possible confounding (e.g., birth weight, sow parity and litter size). Of these 192 female piglets, a total of 96 piglets were randomly assigned to either treatment (ROM) or control (Ctrl) groups, while groups were matched in terms of sow parity, time of birth and body weight (BW). Six female piglets per pen with 16 pens (and sows) and 48 piglets per treatment were distributed over four farrowing rooms. Males and female piglets that were not included in the experiment were also equally divided over pens and housed together with experimental piglets.
The study lasted 70 days starting at the birth of piglets, and eight pigs from ROM and control groups were euthanized and dissected on day 27, i.e., one day prior to weaning, and on days 44 and 70 (Fig. 1). During the suckling phase, all piglets were housed together with their littermates and respective biological or foster mother (pen = 5.21 m² of which 0.7 m² closed with tender slats in the piglet area and stainless-steel slats for the sow crate area) and adequate enrichment. The farrowing room was equipped with a computer-controlled climate system that was set to thermoneutral for the sow area, while piglets had a nesting area with floor heating and lamps. Starting three days before weaning (day 25), piglets had access to creep feed (weaner diet; Tab. S1) to familiarize them with solid feed prior to weaning. On day 28, the piglets were weaned, and a subset of experimental pigs (48 pigs, 3 per pen) were randomly selected and moved to corresponding nursery pens. In all pens solid feed and water were available ad libitum. A weaner diet was provided from days 29 to 44, and a nursery diet from days 45 to 70 (Tab.S1).
The intervention with ROM started on day 2 and ended on day 44 (Fig. 1). The product ROM contained around 40% mycelium of A. subrufescens, and was provided by Selko (Trouw Nutrition, Tilburg, The Netherlands). The material was freshly prepared and thoroughly mixed with tap water immediately before feeding, with the concentration of 100 mg/ml during the first week, followed by a concentration of 200 mg/ml during the remaining intervention period. During the suckling phase, piglets received ROM every other day as an oral drench (starting with 100 mg/day and doubling each week till 800 mg on day 28). An equal volume of tap water was used for the placebo (Ctrl group) using disposable syringes (Discardit II, BD). During the post-weaning phase, the dosage was increased to 1000 mg/day per administration from day 29 until day 44, after which the dietary intervention was stopped. To assess the effect of ROM on immune function, an oral live attenuated vaccine (Salmoporc® STM (lot number 0270617), IDT Biologika GmbH, Dessau-Rosslau, Germany) of attenuated S. Typh was orally administered on day 21 and day 45 according to manufacturer’s instructions. The vaccine suspension was freshly prepared according to manufacturer’s instructions before each oral administration. In addition, BW of piglets was measured at several time points (days 1, 14, 26, 35, 43, 59 and 69). The feed intake of each pen during nursery period was also recorded, and the average daily feed intake (ADF) and feed conversion ratio (FCR) were calculated.
Faeces, gut luminal digesta sampling and microbiota profiling
We collected faeces and luminal digesta during both the pre-weaning and post-weaning phases to evaluate the effect of ROM intervention on microbial composition. Fresh rectal faeces were sampled in cryotubes by using a wetted (with sterilized H2O) Puritan PurFlock Ultra cotton swab (Puritan, ME, USA) at eight time points (days 4, 8, 14, 26, 35, 43, 59 and 69) (Fig. 1) and immediately placed in a box with dry ice, followed by stored at -80 ℃ until further analysis. Furthermore, at each dissection day (days 27, 44 and 70) (Fig. 1), approximately 1 g of homogenized jejunal, ileal and caecal digesta per pig were collected in a sterile cryogenic vial, snap-frozen in a box with dry ice and afterwards stored at -80 ℃ until further analysis. The remainder of the digesta from each GI tract segment was mixed with sterilized H2O, and the pH was measured using a pH meter (ProLine). The detailed method of sample collection on each dissection day can be found elsewhere [15].
The microbial composition was analysed by sequencing of PCR-amplified and barcoded 16S ribosomal RNA (rRNA) gene fragments on the Illumina NovaSeq 6000 S2 PE150 XP platform. Faeces sampled from piglets that were dissected and all luminal digesta collected were selected for total DNA isolation with repeated bead beating for mechanical cell disruption [16]. In brief, faeces (~50 mg) or luminal digesta (~100 mg) was weighed into a screw cap tube containing sterilized 0.1 mm zirconia beads (0.25 g) and three 2.5 mm glass beads, then 300 μL Stool Transport and Recovery (STAR) buffer (Roche Diagnostics, USA) was added to the tube and followed by homogenization using a bead beater (5.5 ms, 3 × 60 s; Precellys 24, Bertin Techonologies, Montigny-le-Bretonneux, France). Afterwards samples were incubated at 95 ℃ with shaking at 300 rpm for 15 min, followed by centrifugation for 5 min at 16,100 × g at 4 ℃. The supernatant was collected, and the pellets obtained after centrifugation were resuspended in 200 μL STAR buffer and subjected to the same procedure as described above. The two supernatants were pooled, and 250 μL was used for the subsequent DNA purification using the automated Maxwell® 16 Research Instrument (Promega, Madison, USA) as previously described [15]. The purified DNA was eluted with 50 μL nuclease-free water (Qiagen) and the concentration was measured using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA). Purified DNA was adjusted to a concentration of 20 ng/ μL and was used as a template for PCR amplification with primers 515F (5’-GTGYCAGCMGCCGCGGTAA) and 806R (5’-GGACTACHVGGGTWTCTAAT) [17], targeting the V4 region of the bacterial and archaeal 16S rRNA gene. DNA from faecal samples was amplified in triplicate in 35 μL reactions that contained 25.5 μL nuclease free water (Promega, Madison, WI, USA), 7 μL 5 × HF buffer (Thermo Fisher Scientific, Vilnius, Lithuania), 0.7 μL of 10 mM dNTPs (Thermo Fisher Scientific), 0.35 μL DNA polymerase (2 U/ μL) (Thermo Fisher Scientific), 0.7 μL of 10 μM sample-specific barcode-tagged primers and 0.7 μL purified DNA (20 ng/ μL) template. PCR reactions with luminal digesta DNA were performed in duplicate in a total volume of 50 μL, and the template DNA was at the same concentration (20 ng/ μL) as described previously [15]. The detailed PCR program can be found elsewhere [18], and was used with minor modifications, i.e., with an annealing temperature at 50 ℃ for all samples, and with 30 cycles for jejunal digesta (as opposed to the 25 cycles used for other samples). Replicate PCR products were pooled and purified by using the CleanPCR kit (Clean NA, The Netherlands). Concentrations of purified DNA amplicons were determined with the Qubit BR dsDNA assay kit (Invitrogen by Thermo Fisher Scientific, Eugene, OR, USA). Finally, equimolar amounts of purified PCR products were pooled into libraries and sent for Illumina Hiseq sequencing with a read length of 2 x 150 bp ((Eurofins Genomics, Ebersberg, Germany). In each library, two artificial mock communities, biological replicates of random samples and a blank (water) were included for quality control. Raw sequence data was first processed using NG-Tax 2.0 [19] with default settings and assigned to amplicon sequence variants (ASVs) using the Silva132 reference dataset [20]. ASVs with a relative abundance lower than 0.1% in a given sample were excluded on a per-sample basis.
Gut mucosal scrapings and microfluidic qPCR for gene expression analysis
Mucosal scrapings from ileal and caecal epithelia were collected on two dissection days (days 27 and 44) (Fig. 1) to investigate the effect of ROM supplementation on pig GI tract mucosal immunity. Around 5 cm long ileal and caecal segments adjacent and proximal to the GI tract segments sampled for microbiota analysis were excised and then longitudinally cut open. Digesta were removed and segments were carefully rinsed with Phosphate Buffered Saline (PBS) without removing the mucus layer. Thereafter, the mucosa excluding the muscularis layer was removed by scraping with a scalpel, added to a snap-lock tube containing RNA later (Qiagen, Hilden, Germany), and immediately placed in a box with dry ice. Samples were shipped to the laboratory on dry ice and then stored at -80 ℃ until further analysis. For RNA isolation, a 3 × 3 mm piece of tissue for each animal was cut off, added to a 2 ml Eppendorf tube containing 500 μL RLT lysis buffer (Qiagen, Hilden, Germany), and homogenized using a Turrax (IKA-Werke GmbH, Staufen, Germany) for 90 s. Subsequently, 100 μL of homogenized tissue suspension was added to a new 2 ml Eppendorf tube containing 600 μL fresh cold RLT buffer and homogenized by pipetting 10 times. Total RNA was isolated using the RNeasy Mini Kit (Qiagen, Hilden, Germany) according to manufacturer’s instructions, including on-column DNase digestion for 15 min. Qubit BR RNA Assay Kit (Thermo-Fisher, Scientific) was first used for preliminary determination of RNA concentration, degradation, and contamination. A Qsep 100 bioanalyzer (GC Biotech, Waddinxveen, The Netherlands) was then used to determine RNA quality and integrity. RNA was immediately stored at -80 °C after isolation until cDNA synthesis. Reverse transcription was performed using the QuantiTect Reverse Transcription Kit (Qiagen, Hilden, Germany) following manufacturer’s instructions. cDNA samples were stored at -20 °C until further use.
We selected 96 genes for microfluidic qPCR analysis based on their functional relevance to GI tract mucosal barrier function and various parts of the mucosal immune response. qPCR was performed in a BioMark HD Reader and the 96.96 Dynamic Array (Fluidigm, CA, USA). Primer pairs (n=96) were designed using Primer3 version 0.4.0 (https://bioinfo.ut.ee/primer3-0.4.0/) for selected genes (Tab. S2). Distribution and length of introns and exons were identified using Ensembl (https://www.ensembl.org/index.html), and interspecies-variation was checked if needed, using BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Primers were synthesized at Sigma Aldrich. Pre-amplification of cDNA was performed using TaqMan PreAmp Master Mix (Applied Biosystems, Waltham, MA, United States) before qPCR. cDNA was pre-amplified in a 10 µL reaction mixture, which contained 3 µL TaqMan PreAmp Master Mix, 2.5 µL 200 nM mix of each of the 96 primer pairs, 2 µL low-EDTA TE-buffer (Panreac Applichem, Darmstadt, Germany) and 2.5 µL diluted cDNA (1:10 in low-EDTA TE-buffer). The cycling conditions were as follow: 95 °C for 10 min, followed by 19 cycles of 95 °C for 15 s and then 60 °C for 4 min. After pre-amplification, cDNA was treated with 16 U Exonuclease I (New England Biolabs, Ipswich, MA, USA) at 37 °C for 30 min, followed by 80 °C for 15 min, and then was diluted 1:10 in low-EDTA TE-buffer for microfluidic qPCR analysis. Microfluidic high-throughput qPCR was performed on a BioMark HD real-time instrument (Fluidigm) as previously described [21].
Blood, serum, and mesenteric lymph node sampling and corresponding analysis
Blood and tissue samples were collected both pre-weaning and post-weaning to evaluate the effect of ROM intervention on the immune response. Blood samples were obtained from the jugular vein using Sodium Heparin tubes or Serum Gel tubes (S-monovette®, Sarstedt, Germany) on days 14, 26, 35, 43 and 69 (Fig. 1). Blood samples collected in Sodium Heparin tubes were used for Flow Cytometry analysis and kept at room temperature (RT) until use. Blood samples collected in Serum Gel tubes were centrifuged at 2000 × g for 10 min to collect serum, which was stored at -20 ℃ until use. The ileocecal Mesenteric Lymph Node (MLN) was collected after euthanasia and exsanguination on days 27, 44 and 70 (Fig. 1) and stored in ice-cold RPMI 1640 Medium (with GlutaMAX™ supplement, Gibco®) that contained 1% L-Glutamine (Gibco®) and 10% foetal calf serum (FCS, Gibco®).
Serum samples were collected to measure vaccine-specific antibodies against S. Typh, for which an in-house ELISA was optimized. A detailed procedure can be found elsewhere [15]. In brief, S. Typh bacteria were recovered from the vaccine Salmoporc®, which contains live attenuated S. Typh. A freshly prepared 100 μL S. Typh bacteria suspension (2 × 108 cells/mL) was used to coat medium-binding 96 well plates (clear, flat bottom, Greiner Bio-One), and incubated overnight at 4 °C. Then, the bacterial suspension was removed, and bacteria still attached to the plate were fixed with 4% paraformaldehyde for 2 h at RT. The plates were blocked with a blocking solution consisting of 5% milk powder (ELK, FrieslandCampina, Amersfoort, The Netherlands) in demineralized water overnight at RT, and stored at 4 °C until usage. Before adding serum samples, blocked plates were washed with PBS/Tween20 (0.05%). Serum samples were diluted 250x (IgG) and 50x (IgA, IgM) in blocking solution, then 100μl of diluted serum samples were added to the plates and incubated for 1 h at RT. After incubation, the plates were washed two times, and 100 μL of 50,000 times diluted (in blocking solution) horseradish peroxidase (HRP) conjugated goat anti-Porcine IgG, IgM, or IgA (Novus Biologicals) was added. After 30 min, plates were washed five times and incubated with 100 μL of 3,3’,5,5’- tetramethylbenzidine (TMB) substrate solution (Enhanced K-Blue®, Neogen) for 15 min. The reaction was stopped with 100 μL of 2% HCl, and the optical density (OD) of the plates was measured at a wavelength of 450 nm (Multi-Mode Microplate Reader FilterMax F5).
To measure the effects of ROM on MLN immune cells, MLN cells were restimulated with 10, 1, or 0.1 µg/mL of lipopolysaccharide (LPS, serotype O55:B5/L2880, Sigma-Aldrich) for 24 h. Next, the levels of cytokine production (IL-10) were quantified by ELISA as previously described [22]. Identification of Natural Killer (NK) cells in isolated PBMCs was performed by Flow Cytometry, followed by analysis in FlowJoTM software (Version 10) [15]. The gating strategy for NK cell identification was designed according to previous studies [23-25].
Statistical analysis
For the growth performance data, BW was analysed with a linear mixed model for repeated measures analyses.
For microbiota data, ASV read counts were first transformed to relative abundance, and all statistical analyses were performed in R 3.6.1 [26]. Alpha diversity, with metrics of observed species, phylogenetic diversity, Shannon diversity and inverse Simpson (InvSimpson), was determined at ASV level using packages picante [27] and microbiome [28], followed by a Linear Mixed-Effects Model that was applied to assess whether alpha diversity was significantly different between ROM and Ctrl piglets over time. Beta diversity was estimated at ASV level based on weighted UniFrac [29] and unweighted UniFrac distance [30], and the results were further plotted by unconstrained principal coordinate analysis (PCoA) using the phyloseq R package [31]. Both metrics take the phylogenetic relationships among ASVs into account, with unweighted UniFrac only considering the presence or absence of ASVs, whereas weighted UniFrac takes the relative abundance of ASVs into account. The significance of differences in beta diversity was computed by permutational multivariate analysis of variance (PERMANOVA), using the Adonis’ function in the vegan package [32]. If there were significant differences in overall microbial composition between two groups, the linear discriminant analysis effect size (LEfSe) method [33] was used to identify biomarkers characterizing differences between ROM and Ctrl groups. A nonparametric Wilcoxon rank sum test was applied to assess whether the inter-individual variation of weighted and unweighted UniFrac distances was significantly different between both groups at each time point. In order to test for differences in relative abundance of microbiota at genus level between ROM and Ctrl groups over time, a generalized additive model for location, scale, and shape with a zero-inflated beta family (GAMLSS-BEZI) was used, as implemented with the taxa.compare function in the metamicrobiomeR package [34]. Significance was assessed with a false discovery rate correction for all multiple testing according to the procedure by Benjamini–Hochberg, with a threshold of 0.05 [35]. Differences were considered significant if adjusted p ≤ 0.05, or as a trend if the adjusted p value was below 0.1 but above 0.05 (0.05< p ≤ 0.10).
For microfluidic qPCR data of gut mucosal gene expression, the method for data pre-processing was as described previously [36]. In brief, data was first pre-processed, normalized and relatively quantified per time point and tissue type, using GeneEx5 (MultiD, Göteborg, Sweden). The algorithms geNorm [37] and NormFinder [38] were then used for reference gene selection per tissue origin (ileum or caecum). Primer pairs that resulted in inconsistent replicates were excluded from the dataset. The Cq values were converted to relative quantities for each primer assay, and relative quantities were Log2-transformed before downstream statistical analysis. Redundancy analysis (RDA) was performed to assess multivariate effects of environmental variables (treatment, time, pH and BW) on ileal and caecal mucosa gene expression, as implemented by the capscale function in the vegan package [32]. Missing values of pH and BW were estimated using the K-nearest neighbour algorithm from the vim package [39]. The statistical significance was estimated by an ANOVA-like permutation test with 999 permutations using the anova.cca function in the vegan package [32]. Differential gene expression between ROM and Ctrl groups per tissue and time point was tested using a linear model framework with an empirical Bayes moderated t-test, and resulting p values were adjusted for multiple testing using the Benjamini-Hochberg false discovery rate (FDR) from the limma R package [40]. –log10 adjusted p value versus log2 fold-change in expression for each pairwise comparison was plotted using EnhancedVolcano R package [41]. Differences were expressed as significantly down- or up-regulated if adjusted p ≤ 0.05 with a log2 fold change > 0.5 or tendencies if 0.05 < p ≤ 0.1 with a log2 fold change > 0.5.
To determine the strength of association between gut mucosal gene expression and corresponding relative abundance of microbial groups at genus level, the microbial composition data was filtered with relative abundance and prevalence thresholds of > 0.1% and 40% of the samples, respectively. Correlation analysis was performed according to Pearson’s product moment (95% CI), and the reported p values were corrected for multiple comparisons according to the Benjamini–Hochberg method.
For immunological analyses, a Linear Mixed Model was used to assess the immunological effects over time and the interaction between time and treatment, using R statistical software (version 3.6.2). Additionally, unpaired Student’s t-tests were performed to assess the differences between the treatment groups per time point. Normality of data (Shapiro-Wilk test) and homogeneity of variances (Levene’s test) were checked prior to statistical testing. Skewness values between −2 to +2 were considered acceptable [42]. Extreme outliers (indicated by R) were removed from the analysis. When the normality was not met, data were log-transformed, but presented as untransformed means. Results with an adjusted p-value below 0.05 were considered statistically significant and results between 0.05 and 0.1 were considered a trend.