Caecal Microbiota Composition of Experimental Laying Hens Differs According to Genetic Line and Vaccination to IBV

Interactions between the gut microbiota and the immune system may be involved in the vaccine response. In the present study, we studied the correlations between caecal microbiota composition and the immune response in six experimental laying hen lines harboring different haplotypes at the Major Histocompatibility complex (MHC), 7 weeks after their rst vaccination against the infectious bronchitis virus (IBV). Two lines were previously considered as high responders (HR) to IBV vaccination and two other ones as low responders (LR). We explored to what extent the gut microbiota could be related to this variability through the characterization of caecal bacterial communities with a 16S rRNA gene amplicon sequencing approach, one week after an IBV infectious challenge.


Background
Microbes are becoming increasingly resistant to the antibiotics used to treat human and animal infections. This is why the use of antibiotics as growth promoters (AGP) in animal production is forbidden in a growing number of countries. It was completely banned in Europe in 2006 (Regulation 1831/2003/EC), and in the USA in 2013 with full effect in 2017 (1). In chicken ocks, the ban of antibiotics was often followed by the observation of more frequent gastro-intestinal disorders such as necrotic enteritis (2). Moreover, genetic selection on production traits such as growth rate, feed e ciency or meat and carcass yield may have indirectly impacted animal robustness, since immune phenotypes or disease resistance phenotypes were generally not included in selection programs by breeders (3,4). For instance in broiler breeding, selection for an increased body weight led to a decrease in the relative weight of primary and secondary immune organs (5). It may also compromise the host immune response to Escherichia coli vaccination (6) or increase the susceptibility to Marek's disease (7). Potential trade-offs between productivity and immunity are often hypothesized to explain such results (8).
In order to avoid economic losses due to infections by pathogens and to ensure animal health in the absence of antibiotics, reinforced vaccination and the inclusion of probiotics in animal feeding were among the rst approaches considered (9). In this context, increasing the e ciency of vaccines may also be considered. Vaccines are still suboptimal for many endemic chicken diseases due to the variety of geographical areas, types of production, and poultry species and breeds, which all in uence the outcome of the infection (10)(11)(12).
Several approaches have been considered to enhance the animal immune response to vaccination, and some may depend on the intestinal microbiota (13). Some of them consider the early manipulation of intestinal innate lymphoid cells (ILCs) (14). Furthermore, many new vaccine designs speci cally target the antigen-presenting cells (APC) and use bacterial components as adjuvants through ligand-based APC targeting. With this technology, protective antigens are selectively delivered to APCs by coupling them to ligands from bacterial origin such as lipoproteins, peptidoglycan, LPS or lipoteichoic acids which are recognized by pattern recognition receptors (PRRs) from APC. Some applications are even considered in poultry breeding. For instance, TLR ligands have been used as adjuvant in chickens immunized with inactivated avian in uenza IAV H5N2 co-administrated with Salmonella enterica agellin (11,15).
In parallel, the early modulation of the intestinal microbiota is already promoted to improve intestinal health of chickens (16,17). Besides, it is also increasingly considered as a promising strategy to enhance the human vaccine response (13,18). For instance, babies that were given a fermented infant formula with a high amount of Bi dobacteria longum subspecies infantis subsequently displayed a higher abundance of this bacterial probiotic in their intestines, and higher IgA titers after poliovirus vaccination (19). The abundances of this same bacterium in the intestines are also associated with a higher T cell response after tetanus vaccination (20). Other associations between speci c intestinal bacteria and response to speci c vaccines have been reported in human for vaccines against Rotavirus (21), Salmonellaenterica Typhi (22), or Shigella dysenteria (23). In chickens, studies report the effects of microbiota modulatory compounds, and of probiotics, such as Anaerosporobacter mobilis or Lactobacillus reuteri, on the outcome of vaccination against Campylobacter jejuni (24,25). At least one study reports the effect of variations in the commensal gut microbiota on adaptive immune responses in vaccinated chickens, with gut microbiota promoting higher titers of in uenza virus speci c IgM and IgG, and higher interferon gamma expression (26).
While investigating the contribution of intestinal microbiota on the vaccine response is of growing interest, it is also likely there is an effect of vaccination on the gut microbiota in chicken. A few recent studies report the effect of speci c vaccination against Salmonella on their gut microbiota (27).
Considering vaccination against parasites, vaccines against Eimeria that are composed with attenuated live parasites to induce a mild infection, which can explain that theydiminish the microbial diversity and Firmicutes relative abundance in broilers (28,29). The effects of the vaccination on the gut microbiota remain to be studied for a larger set of pathogens, and for different types of vaccine. Nevertheless, whenever a vaccine impacts the gut microbiota, the outcome on the immune system homeostasis and e.g. competitive exclusion, could also be questioned.
In addition to the intestinal microbiota, a parameter usually not considered when developing vaccination strategies is the potential in uence of animal genetics. Adjusting vaccines to the host genetics, or genetically selecting animals displaying a better immune response to vaccines, could be a relevant way of improving the vaccine response in farm conditions. Strong associations occur between the MHC (Major Histocompatibility Complex), which is a polymorphic region comprising genes coding for molecules involved in antigen recognition through cell surface molecules on e.g. antigen presenting cells, and the responses to many infectious diseases in chickens (30). This suggest that genetic variations at the MHC might also control part of the variation in the vaccine responses, as for example regarding the Newcastle disease (31).
In order to evaluate simultaneously the effect of gut microbiota composition and genetic variations on the host response to vaccination, we used several experimental laying hen lines homozygous at the MHC for different haplotypes and studied their response to the IBV (Infectious Bronchitis Virus) attenuated vaccine. IBV is a single-stranded RNA and enveloped Coronavirus causing a respiratory disease in infected chickens. The high rate of mutations of this virus is leading to the existence of many strains with different virulence levels. The morbidity through respiratory symptoms in infected husbandry is high.
Despite vaccination schedules, chickens still display suboptimal levels of protection worldwide. In layer hens, IBV leads to economic losses with both a drop of egg production and a deterioration of egg quality. Typically, egg laying declines range between 3% and 10%, but reductions of up to 50% have also been observed (32).
In the present study, we thus evaluated the association between the caecal microbiota composition and the vaccine response against IBV in layer lines with different MHC haplotypes. We followed, one week after an infectious challenge with IBV, a control group of animals that were not vaccinated and thus displayed an in ammatory response to the challenge, and a vaccinated group of animals that rather displayed a sub-clinical response to the infection. We studied the caecal microbiota composition at the end of the experiment, in order to evaluate the potential association of variations in gut microbiota composition with variations in immune phenotypes following the vaccination against IBV in these lines.
We thus wonder if a part of the MHC line divergence regarding immune response to IBV vaccination (33) could be due to variations in caecal microbiota composition potentially contributing to differences in immune phenotypes. More especi caly, we wanted to decipher whether: At 9 weeks of age, the 96 animals were slaughtered by cervical dislocation in conformity with the Danish regulation for experimental animals. After slaughter, caecal contents were gently removed from the caecal bags so as not to collect the mucosa and immediately frozen in liquid nitrogen before being stored at -80 °C. As described previously by Larsen et al. (33), EDTA-preserved blood were collected at the same time from the jugular vein for immunophenotyping of both non-vaccinated and vaccinated chickens.
DNA extraction, ampli cation and sequencing DNA was extracted from individual caecal contents using a standardized protocol adapted from Godon et al (35). DNA quantities were measured using a Qubit analyser, and DNA qualities using a Nanodrop analyser by measuring ratios of absorbances at 230, 260 and 280 nm. We diluted DNA samples to reach a concentration of 15ng/µL before PCR ampli cation. A 464 bp fragment targeting the hypervariable V3-V4 region from bacterial 16S rRNA gene was rst ampli ed with the following primers: 460-F-iII-16S (adapter-CTTTCCCTACACGACGCTCTTCCGATCT-343-ACGGRAGGCAGCAG-357), and 460-R-iII-16S (adapter-GGAGTTCAGACGTGTGCTCTTCCGATCT-781-TACCAGGGTATCTAATCCT-806), in accordance with primer adapter previously reported by (36). The rst PCR reaction was conducted with 1 µL of genomic DNA and 0.5 µL of each primer (10 µM), 0.5 µL of dNTP mix (10 mM), 2.5 µL of 10X MTP Taq buffer mix (10 mM), 0.25 µL of MTP Taq DNA Polymerase (SIGMA-ALDRICH) and water qsp 25 µL. PCR conditions were: initial denaturation at 94°C for 1 min, followed by 30 cycles of 94°C for 1 min, annealing at 65°C for 1 min, extension at 72°C for 1 min, then a nal elongation step at 72°C for 10 min. We checked by agarose gel electrophoreses that each PCR generated a unique product size at the expected length. Puri cations, second PCR with the Illumina adapters, and sequencing using a 2x300 bp MiSeq Illumina sequencer were performed by the GeT-PlaGe platform (INRAE, Toulouse). Samples from three animals could not be extracted, leading to a set of 93 animals. Raw sequences are available through the accession PRJNA658073 on the NCBI SRA (Sequence Read Archive). We veri ed on electrophoresis gel that after the rst pcr, the two pcr control samples did not produce any visible ampli cation, and that their sequencing depth was far below the mean sequencing depth of the caecal samples (9 582 and 1 643 raw sequences in both pcr control, against an average sequence of 131 789 in the 93 samples). Negative control sequences are available in the additional le (Additional File 1).

Bioinformatic analyses
We used the FROGS pipeline (v3.0) (37), to assemble read pairs, trim the 460-F-iII-16S and 460-R-iII-16S primers, and lter sequences to keep only amplicons with a length ranging from 300 to 490 bp. Unique sequences with their associated abundances in each sample were clustered with FROGS using Swarm (v 1.4.1) (38), rst with an aggregation distance of 1 and secondly of 3. FROGS removes chimera thanks to Vsearch (v2.6.0) (39). Rare OTUs were ltered out (abundance <0.005% of total abundance) (40) and the remaining sequences were annotated by NCBI Blast+ (2.6.0) (41) with an alignment on SILVA database restricted to 16S references (version 132 from 2017 update) (42). We nally obtained an abundance table. We produced a phylogenetic tree using FROGS Tree tool with Mafft (v7.310) (43) and Fasttree (2.1.10) (44). We rare ed counts for α-and β-diversity analyses at 10,000 counts per sample. The unrare ed OTU table, it corresponding taxonomic classi cations, and the metadata are available in the Additional Files 2, 3, and 4 respectively.

Microbiota taxonomic composition and diversity
Bio-statistical analyses were performed using R 3.5.2 and the Phyloseq 1.26.1 package. From the distribution of OTU counts, we described taxonomic composition at phylum, family and genus levels. To identify excluding or co-abundant genera, we calculated Spearman correlations between the abundance of each possible pair of genera on the overall dataset with the Hmisc 4.2.0 package. We represented on a network the moderate and high associations present in at least 20% of the sampled chicken using a cutoff of 0.5 for the Spearman coe cient, and 0.05 for the adjusted p-value.

Effect of vaccination and chicken line on the microbiota
We used the Vegan 2.5.5 package to evaluate the α-diversity by calculating the species richness, and the Chao1, Shannon, and Inverse Simpson indexes that we represented for each vaccine group, each line and each vaccine extreme responder group separately (HR versus LR). β-diversity indices were calculated using the Bray-Curtis (BC) distances. We described β-diversity in each line, and we represented it using a NMDS analysis (Non metrical Multidimensional Scaling). We rst depicted the effect of vaccination by gathering the 6 chicken lines, and then for each vaccine group. We represented the vaccination impact in the different lines, adding on the NMDS ordination the centroids displacements from the control group to the vaccinated in each line. To assess the effects of line and vaccination on α-diversity, we performed an ANOVA on α-indices. We used the vaccine status as a co-variable to assess the line effect, and the line as a co-variable to assess the vaccine-status effect. We did not correct statistical models for a "sex" effect, because only 6 out of 93 animals were males, and because we assessed with an NMDS analysis that the microbial community of the male birds was not dissimilar from the microbial community of the female birds (data not shown). Similarly, we did not correct the models for a putative pen effect, because in this experiment it was confounded with the vaccine status group (data not shown). Signi cance threshold was set at 0.05 on the p-value. To assess the line and then the vaccination effect on β-diversity, we used an Adonis model. Similarly, we used the MetagenomeSeq 1.24.1 package to compare OTU abundances between each couple of lines, and between each couple of vaccine status groups. This package enables the normalization of raw sequences counts through a cumulative-sum scaling (CSS) method to reduce biases due to uneven sequencing depth (45). We then used a zero-in ated Gaussian mixture model with vaccine status or line as main effect to identify, respectively, OTU with an abundance differing between each pair of lines, or rather between each vaccine status, with a threshold of 0.05 on the adjusted p-value (FDR). Such analyses were performed independently at the species and genus levels.

Association between microbiota and immune phenotypes
We searched for associations between individual caecal microbiotas (OTUs) and the immune phenotypes by calculating Spearman correlations between each immune parameter and the abundance of each OTU. To perform it, we used as starting point the normalized abundances from previous MetagenomeSeq normalization. As it appeared that both the vaccination and the line had a signi cant effect on the microbiota, we collected residuals from linear models implemented on those normalized OTU counts with lines and vaccine status as main effect. Similarly, we collected residuals from linear models implemented on each immune parameter with once again the vaccine status and line as main effects. Then, these two set of residuals were used to calculate Spearman correlation coe cients, and we set a cutoff of 0.05 on their adjusted p-values.
We rst performed this analysis on a subset of contrasted lines from the HR group (B2 and B21 lines) and from the LR group (B19 and B12 lines) to have a closer look at the vaccine e ciency. We then repeated this analysis on the whole dataset, and then on each vaccine status group separately. To analyse each vaccine status group separately, we collected residuals from linear models implemented on each normalized OTU count and each immune parameter with only the line as main effect. We thus assessed the vaccine effect from the comparison of the separate results. We represented the intermediate and high associations on a network, adjusting the cutoff on the Spearman correlation in order to keep at least 5 signi cant associations. Thus, we set a cutoff of 0.4 considering the subset on HR and LR group, and on the entire dataset, and a cutoff raised to 0.5 on the control and vaccinated groups separately that displayed more high associations.
We also searched for associations between the microbiota and the immune phenotypes regarding the effect of microbiota β-diversity calculated with BC distances on each immune parameter. We used an Adonis model performed on the entire dataset with line and vaccine status as main effects. In order to identify the lines that possibly contribute more to some of the observed effects, we then performed these analyses on each line separately by implementing the Adonis model only with the vaccine status as main effect. For the immune parameters that displayed the higher association with the microbiota, we represented their level according to the microbiota ordination on a NMDS plot.

Prediction of putative functions carried by OTUs and identi cation of differentially abundant functions between HR and LR group
We realized a metagenome inference from the caecal microbiota using the OTU abundance tables with the PICRUSt2 method implemented with the Integrated Microbial Genome (IMG) database from November 2017 (46). We used the QIIME2 plugin to assign sequences of the identi ed OTUs into a reference multiple-sequence alignment and then into the reference phylogeny to get representative OTUs. From the Kyoto Encyclopedia of Genes and Genomes (KEGG) database as genome functions reference, we thus obtained a prediction of KEGG Orthology (KO) abundances. Using the DESeq2 1.22.2 R package, KO abundances for each sample were then submitted to a model with the vaccine e ciency group (HR vs LR) as main variable and the line as co-variable. Following proper data normalization, the DESeq2 method, originally developed for RNAseq data analyses, performs similarly to or better than many other algorithms developed speci cally for microbiome data (47,48).
We then performed a Wald test to identify sets of KO functions signi cantly more abundant in either the HR or the LR group, with an adjusted p-value threshold of 0.05 through Benjamini-Hochberg method. We then used the clusterPro ler 3.10.1 R package to obtain the pathways, KEGG identi ers, and abbreviations corresponding to these two sets of KEGG functions. In order to search for possible immune interactions, we only considered from the output some genes coding for epitopes with known or possible immunogenicity. From this prediction, we thus looked for differentially abundant expression of genes linked to lipopolysaccharide, peptidoglycan, ABC transporters, and N-glycan synthesis or transport.

Taxonomic composition
From the 93 caecal samples processed, 16S rRNA gene sequencing provided a total of 12,256,442 reads that after quality control led to 4,952,178 sequences with an average of 53,249 sequences per sample (range: from 32,560 to 76,819 sequences). We identi ed 487 OTUs from 5 phylum (Suppl Table S1 Table S1-c). The massively dominant phylum was Firmicutes (average abundance 96%), with essentially Lachnospiraceae (48.1%), Ruminococcaceae (42.1%), and in a lesser extent Bacilli (2.1%). It was followed by the much less abundant Proteobacteria (average abundance 3.9%). Over the 73 genera identi ed, Ruminococcus torques group (12.4%), Eisenbergiella (9.0%), and Ruminococcaceae UCG-014 (8.0%) were the three dominant ones, and 34 signi cant correlations were observed between their abundance (Suppl Table S2, Suppl Figure S2). Eisenbergiella was negatively associated with Ruminococcaceae UCG-014 (rho=-0.52), and Escherichia-Shigella, which is often described as pathogenic, was positively associated with Anaerostipes (rho=0.54), and with an unknown genus from the Gammaproteobacteria class (rho=0.95).  Table 3, Additional File 5), meaning that the microbial dissimilarity observed between the lines was not due to a higher bacterial richness in some of them.
Results from pairwise comparisons of BC distances between lines showed that the line B14 displayed the greatest dissimilarity with the other ones (all pairwise p-values < 0.0002), and the greatest homogeneity within its bacterial community (Table 1, Figure 2b-d). Then, the line B21 was the more dissimilar one, followed to a lower extent by B2 and B12. Furthermore, we observed a closer proximity between B12 and B14, and between B19 and B21, which did not match the lines clustering regarding their vaccine e ciency. The number of OTUs differing between a pair of lines were in accordance with their corresponding β-diversity pro les as it ranged from none (B12 -B15 and B2 -B21), to 126 (B14 -B21) (Suppl Table S4-a). Globally, the line B14 had more abundant OTUs from the Ruminococcus torques group, and fewer abundances of OTUs from Ruminococcaceae UCG-014 and Ruminiclostridium 5 genus (Suppl Table S5).
The B21 line is mostly dissimilar with B12 and B14, and it globally harbored a larger abundance of OTU from Ruminococcaceae UCG-014 and GCA-900066575 genera and a lower abundance of OTU from Flavonifractor, Escherichia-Shigella, and Oscillibacter genera (Suppl Table S5). Finally, the B12 line harbored a higher abundance of OTU from Escherichia-Shigella genus in comparison with the lines B2, B21 and B19 (Suppl Table S5).
The same pattern was observed considering the vaccinated and control groups separately: we did not observe a line effect on α-diversity in any group, and once again, the lines B14 and B21 displayed the more dissimilar microbiota in both groups (Suppl Table 4-b-c). However, we observed that vaccination brought the lines closer: they were less dissimilar in the vaccinated group than in the control group, especially for B21 (all p-value<0.006 in pairwise comparison on BC distance implying B21 in the control group, while only two comparisons displayed p-value<0.03 in the vaccinated group, with signi cance maintained with B14 and B12). Furthermore, vaccination also brought closer the microbiota inside each vaccine e ciency group: in the LR group, the Bray-Curtis dissimilarities observed between lines B19 and B12 in the control group (p-value=0.01) disappeared in the vaccinated group (p-value=0.11); in the HR group, the BC dissimilarities between lines B21 and B2 line in the control group (p-value=0.0008) disappeared with vaccination (p-value=0.23).
Vaccination effect on the caecal microbiota of IBV challenged chicken Vaccination was associated with a reduced species richness (p-value=7.0 x 10 -4 ) and Chao1 estimator (p-value=7.2 x 10 -4 ), but it did not had any signi cant effect neither on Shannon (p-value=0.24) nor on inverse Simpson (p-value=0.92) ( Table 2, Suppl Figure S4). Since Shannon index is more sensitive to the more abundant OTUs, we can hypothesize that the greater richness observed in the control group is due to its higher carriage of rare OTUs. As Lachnospiraceae and Ruminococcaceae were the dominant families observed, we wondered whether vaccination had also an effect on species Richness and Chao1 in each of them separately. In the Ruminococcaceae subgroup, vaccination signi cantly decreased OTU richness (p-value=1.1 x 10 -7 ) and Chao1 (p-value=9.8 x 10 -7 ), but it had no signi cant effect in the Lachnospiracea group.
Regarding β-diversity, vaccination had a signi cant effect on Bray-Curtis distances (p-value=1x10 -4 ) with vaccinated versus control samples forming two distinct clusters noticeable on NMDS representation (Figure 3-a). Furthermore, it strengthened microbiota homogeneity as it signi cantly reduced β-diversity which is measured through the mean average distance of each sample to the centroid of the corresponding group (p-value=4x10 -3 ). Finally, by looking at the centroids displacement on an NMDS from the control group to the vaccinated group, vaccination similarly disturbed the microbiota β-diversity in all lines (Figure3-b).
Differential analysis on OTU abundance carried considering all lines together showed that most of the OTUs impacted by the vaccination were from the Ruminococcaceae and Lachnospiraceae families, and to a lower extent from Clostridiales vadinBB60 group, and Enterobacteriaceae family (Suppl Table S7). Twenty-three out of 25 differentially abundant OTUs from the Ruminococcaceae UGG-014 group genus, and 5 out of 6 OTUs from Faecalibacterium genus were less abundant in the vaccinated group. On the contrary, 11 OTUs out of 12 differentially abundant OTUs from Eisenbergiella genus, 5 out of 7 OTUs from Ruminococcus torques group, 5 out of 6 OTUs from Ruminoclostridium 9 and Oscillibacter genera, all 3 OTUs from Escherichia-Shigella, Lachnospiraceae NK4A136 group, Subdoligranulum, and Lactobacillus genus, and both OTUs from Tyzerella 3 genus were more abundant in the vaccinated group.
Considering the lines separately, the trend for Ruminococcaceae was similar, as most of the OTUs signi cantly different in this family were less abundant in the vaccinated group (Suppl Table S8). Ruminococcaceae UCG-014, the most affected genus within this family, was especially impacted in the lines B2, B19, and B21. Finally, vaccination disturbed more the B2 line (39 OTUs differentially abundant), then the B19 line (22 OTUs), and to a lower extent the B15 and B21 lines (9 OTUs).
Association between vaccine e ciency and caecal microbiota Microbial differences between high and low vaccine response lines We compared the HR (B2 & B21 lines) and LR (B19 & B12 lines) groups to study the association between vaccine e ciency and caecal microbiota (n chicken =62, n OTU =487). There was no signi cant difference for the richness indices observed between these groups (Suppl Figure S5), which suggests that the better vaccine response from HR animals was not necessarily correlated to a higher caecal microbiota richness, however there was a signi cant β-diversity difference on Bray-Curtis distances (p-value=0.001, Suppl Figure S6-a). This microbial dissimilarity was also present in the control group (n vac =48, Adonis p- In order to search for potential bacterial species contributing to this observed dissimilarity, we searched for OTUs with different abundances between the HR and LR groups whatever the vaccine status (Table  3). We observed 7 bacterial species from the Lachnospiraceae family, and 4 bacterial species from the Ruminococcaceae family that were signi cantly more abundant in the HR group, against 4 OTUs from the Lachnospiraceae family, 5 from the Ruminococcaceae family, 1 from the De uviitaleaceae family, and 1 from the Escherichia-Shigella genus that were more abundant in the LR group. Among the most signi cant differences, two OTUs were more abundant in the HR group than in the LR group: one from the Tyzzeralla genus (abundance ratio HR/LR of 4.6, adjusted p-value= 0.003) and one from the Angelakisella genus (abundance ratio HR/LR of 7.6, adjusted p-value= 0.003). On the contrary, the genus CAG-352 and Escherichia-Shigella were both more abundant in the LR group (resp. abundance ratio LR/HR of 3.6 and 2.3, adjusted p-value of 3.39x10 -3 and 2.09x10 -2 ). The Cluster_143 from Flavonifractor genus was more abundant in the HR group only considering the vaccinated animals (adjusted p-value=0.03).

Association between microbiota and immune phenotypes in LR and HR groups
We also searched for associations between the microbiota and the immune phenotypes regarding the effect of the microbial β-diversity on each of 28 immune parameter with an PERMANOVA (Adonis) model performed between the LR and HR groups. We found 9 phenotypes that were signi cantly different between HR and LR bacterial communities considering Bray-Curtis distances ( Table 4). TCR γδ expression on TCRγδ T cells was the immune parameter with the most signi cant difference (Adonis p-value=2x10 -4 ), followed by monocyte counts (p-value=0.001), CD45 expression on heterophils (p-value=0.002), and  (Table S9). It is to notice that among those few high associations, two implied a negative correlation between OTUs from Escherichia-Shigella genus and TCR γδ expression on TCR γδ + T cells (Cluster 3: rho=-0.60, adjusted-p-value=3.5x10 -7 ; Cluster 544: rho=-0.55, adjusted-pvalue=4.3x10 -6 ).
The two immune parameters that displayed the greatest numbers of high or moderate associations with OTU abundance were also the immune parameters that had the most signi cant differences based on Adonis model performed on the HR and LR Bray-Curtis distance. We indeed observed 12 OTUs that shared signi cant moderate or high associations with TCR γδ expression on TCR γδ + T cells (5 being positive), and 10 with monocyte counts (4 being positive) (Table S9). Among those OTUs, Cluster 65 from Ruminococcacea UCG-014 genus displayed a negative association with frequency of TCR1+ CD8+ expression on TCRγδ T cells (rho=-0.44, adjusted-p-value=3.5x10 -4 ) and with CD4 expression on CD4 T cells (rho=-0.42, adjusted-p-value=8.0x10 -4 ), and was also signi cantly more abundant in the HR group than in the LR group (adjusted-p-value=2.7x10 -2 ) (Table 3). Similarly, Cluster 759 from Escherichia-Shigella genus that displayed a negative association with CD4 expression on CD4+ T cells (rho= -0.43, adjusted-p-value=5.0x10 -4 ) and was also more abundant in the LR group (adjusted-p-value=2.1x10 -4 ). About TCR γδ expression on TCRγδ T cells, among the 12 positive associations, 6 implied OTUs from the Escherichia-Shigella genus (Table S9). CD8 αß expression on TCRγδ T cells and CD4 expression on CD4 + T cells both displayed 6 associations with OTU abundance from the microbiota. Considering the monocytes counts, among the 6 negative associations observed, 2 involved OTUs from the Flavonifractor genus, that were both more abundant in the B12 line than the B21 line (Cluster 14, Cluster 625, adjusted p-values<0.04).

Prediction of microbial functions differentially abundant between HR and LR groups
Picrust2 enabled us to predict the presence of 4.509 functions in caecal bacterial communities of both HR and LR groups, with 268 KEGG genes predicted to be more abundant in the HR group, against 1451 in the LR group (Suppl Table S10). From the enrichment analysis it appears that in the HR group there was a signi cant enrichment in predicted functions related to quorum sensing (adjusted p-value=0.001), with 80 KEGG genes from this pathway observed to be more abundant in HR than LR animals (Suppl Table S11).
In the LR group there was a signi cant enrichment in predicted functions related to cationic antimicrobial peptide (CAMP) resistance (adjusted p-value=1.2x10 -9 ), two-component system (adjust-p-value=7.5x10 -8 ), Escherichia coli bio lm formation (adjusted-p-value=7.5x10 -8 ), lipopolysaccharide biosynthesis (adjusted-p-value=1.3x10 -7 ), bacterial secretion system (adjusted-p-value=1.6x10 -3 ), Vibrio cholerae bio lm formation (adjusted-p-value=5.1x10-3 ), phenylalanine metabolism (adjusted p-value=2x10 -2 ) and degradation of aromatic compounds (adjusted-p-value=4.9x10 -2 ). From this list of "LR" enriched gene functions, we then focused on lipopolysaccharide biosynthesis and bacterial secretion system pathways. We found some differentially abundant genes mostly for functions related to lipopolysaccharide biosynthesis, then to ABC transporters, and to a lower extent to peptidoglycan biosynthesis, and n-glycan biosynthesis (Suppl Table S10 Associations between immune phenotypes and microbiota according to the genetic line We searched for associations between individual caecal microbiota and immune phenotypes with ADONIS. We rst used this model between each immune parameter and the Bray-Curtis distances of the caecal microbiota from the whole population, considering the 6 lines together. Then, we repeated these analyses in each line separately to search if some of them could contribute more than others to the observed associations. We found that 5 immune parameters had a signi cant association with the caecal microbiota. As we Associations between immune phenotypes and microbiota according to the vaccine status Finally, we searched for the bacteria that could be involved in these associations through Spearman correlation analysis that we calculated between the level of each immune parameters and the relative abundance of each OTU. We rst calculated them on the entire population to search for association present whatever the vaccine status and the level of IBV infection ( Figure 5-a, Table S12). Then, we calculated it on the animal from vaccinated ( Figure 5-b, Table S13) and control group ( Figure 5-c, Table  S14) separately to identify association speci c of an acute infection with an in ammatory state group (control group) or of a rather sub-clinical infection (vaccinated group).
We only get few common associations between these three analyses. We observed negative association between TCR ϒδ expression on TCRγδ T cells and OTU from Escherichia-Shigella genus in the three datasets. Especially, Cluster 544 from this genus was observed in the three datasets whatever the level of infection. Such consistency strengthens the likelihood of a biological link between Escherichia-Shigella and TCR ϒδ expression on TCRγδ T cells.
All the other associations we observed on the vaccinated group and the control group implied quite different immune parameters and different OTUs.

Discussion
We report the combined effects of vaccination to IBV and of different MHC genotypes of experimental laying hens on their caecal microbiota in a context of IBV infection. Experimental IBV challenge enabled us to fully follow the capacity of birds to decrease the viral load and thus their level of protection and recovery as IBV-speci c immunoglobulin level per se is not a good predictor of protective immunity (49,50). Once established the respective effects of the line and of the vaccination on the caecal microbiota, we investigated the associations between the caecal microbiota composition and immune phenotypes of the challenged animals. A previous study allowed us to identify two extreme groups (LR/HR) among the 6 genetic lines, according to their response to vaccination (33). We used this classi cation to analyse the gut microbiota, by investigating whether the microbiota from chicken lines with different vaccine response levels could contribute to the immune response during an IBV infection.
To investigate it, we considered those 6 genetic lines of mainly females layers reared together.
Animals were studied one week after an intranasal viral infection, which corresponds to the acute phase of the disease in non-vaccinated animals, and to a sub-clinical or recovery phase of the disease in their vaccinated counterparts. To our knowledge, there are no published reports on the effect of IBV infection on the intestinal microbiota. Some elements on richness drop have been proposed based on research performed on mice infected with either Respiratory Syncytial Virus (RSV) or In uenza: gut microbiota of treated and control mice were dissimilar, and challenged mice had a reduced intestinal richness (51,52).

An association between genetic lines and caecal microbiota composition
We observed that while the microbiota of the six genetic lines did not differ for their α-diversity, the microbial communities were dissimilar in their global composition, as assessed by the comparisons of their Bray-Curtis distances and the differential analysis performed on their OTU abundances. In our experimental design, animals from the different lines were reared in the same pens and were thus pecking a litter that harbored faeces from all the lines (53,54). Those environmental conditions could have favored more microbial similarity between the lines, which gives more credit to the host genetic effect we observed. Intestinal microbial dissimilarities are already reported among different breeds in poultry, and also among different broiler lines that diverge on their digestive e ciency, or on their susceptibility to Salmonella Enteritidis and Campylobacter jejuni (55,56). These microbial dissimilarities between lines of animals reared in the same environment could have resulted from physiological particularities related to their immunological immune system (57), or from differences in their digestive anatomy (58-61). The Page 16/31 differences observed between lines could be due to variations in the MHC locus but not only: although they are all White Leghorn, the lines we studied are not isogenic. They may differ for genomic regions other than the MHC locus, which could also have an effect on the caecal microbiota composition. Progenitors of the studied animals were reared in the same environment, so that we can discard the hypothesis made in other studies of a selecting pressure exerted by the environment through yolk antibodies (53,62). This reinforces the hypothesis that genetic differences shaped the differences of caecal microbiota between lines through an action on immunity, because genetics studies have shown that some quantitative trait loci (QTL) involved in immunity are associated to the higher abundance of speci c bacterial species (63). The present lines differ for their MHC, which has a pivotal role in the adaptive immune response. Indeed, MHC may shape the intestinal microbiota through its in uence on the antibody response, on some CD4+ Th cell subsets, or on innate lymphoid cells that recognize commensal antigens. Its involvement in the gut enterocytes that secrete antimicrobial molecules able to in uence microbiota composition is also suspected (64-66). Furthermore, differences in both a nity and number of amino acids recognized in binding pockets between alleles of MHC class I molecules in birds is a new concept that could explain the association between MHC haplotypes and some disease resistances (67).
Their MHC motif can indeed be considered as "promiscuous" or "fastidious". Alleles leading to a lower expression of MHC class I molecules having a lower a nity but a greater ability to bind a large variety of peptides are identi ed as promiscuous (generalist) against a wide variety of infectious pathogens. At the opposite, fastidious alleles display a higher expression of MHC class I molecules that can bind to a reduced variety of peptides, which can support an advantage against some speci c pathogen. In a study on wild seabirds, but with limited sampling, an effect of the MHC on the gut microbiota was observed (68). In stickleback sh, 6 out of 14 MHC motifs from class II molecules were signi cantly associated with the abundance of 7 bacterial families, and shes with more diverse MHC motifs had less diverse gut microbiota(65). Even if it is di cult to compare effects of the MHC on the gut microbiota between poultry and sh, which harbor quite dissimilar bacterial communities and environments, we can notice that the Enterobacteriaceae family was prone to sex-speci c MHC effects in stickleback sh. In a similar way, we also observed in the present study that OTUs from this family displayed higher relative abundance in the B14 and B12 lines. Thus, despite the few studies, an effect of MHC genotype on the Enterobacteriaceae abundance is possible. Furthermore, it is a group of interest as many bacteria from this family are pathogenic and lead to intestinal in ammation (69). In mice, 189 OTUs displayed differential abundances between animals homozygous or heterozygous at the MHC locus (66). As we observed in the present study, OTUs from Lactobacillus, Anaerostipes, Ruminococcus, and Butyricicoccus genera were differentially abundant between different MHC groups(64). In poultry, and more speci cally in layers, only one similar study was performed on the caecal microbiota from two congenic inbred lines that differ by their MHC genotypes (70). As we observed in the present study, they did not noticed any signi cant difference in α-diversity metrics in the microbiota composition of the two lines. To our knowledge, there is no other publication about the association between MHC haplotypes and gut microbiota in poultry and the present study thus brings some promising elements to deeply investigate through a detailed characterization of the MHC region.

A clear association between vaccination and caecal microbiota composition
We observed that the bacterial richness was lower in vaccinated chickens upon challenge as compared to challenged unvaccinated chickens, especially in the Ruminococcaceae family. Moreover, the vaccination effect was similar in all the lines, with globally a reduced abundance of OTU from the Ruminococcacea UCG-014 and Faecalibacterium genera, and an increased abundance of OTU from the Eisenbergiella genus and to a lower extent of Ruminococcus torques group, Tyzzerella 3, and Ruminoclostridium 9 or 5 genera. Due to the IBV viral challenge, we can suppose that the vaccination effect we observed was mostly the consequence of a different level of infection between the control and vaccinated groups.
Animal clinical signs showed that animals from the control group were in the acute phase of the disease with an increased body temperature, a reduced appetite and a higher in ammation state, while the vaccinated animals rather displayed sub-clinical symptoms. Thus, we can not exclude that the differences we observed on the caecal microbiota between the control and the vaccinated group may have resulted from the difference of those clinical state. Independently of a disease outcome, there are still only a few reports on the effect of the vaccination on the microbiota, and the results are contradictory (51,71,72). Investigating the effect of vaccination on gut microbiota is thus of interest to improve vaccine design, and the present study provides promising results to also pursue it with unchallenged animals.
A set of correlations between immune phenotypes and caecal microbiota composition TCR ϒδ expression on TCRγδ T cells was the immune parameter with the greatest consistency between the different analyses performed. It was particularly associated with the caecal microbiota of the lines B21, B14 and B15. Moreover, Larsen et al. observed a large variation in blast transformation between MHC-B lines for γδ T cells (33). Thus, our results extend to the microbiota the phenotypic differences already observed between MHC-B lines. TCR ϒδ expression on TCRγδ T cells also shared negative associations with many OTUs from the Escherichia-Shigella genus, in both the control and the vaccinated group, whether restricting the analysis to extreme responders to vaccination or not. Gamma delta T cells are a subpopulation of T lymphocytes with TCR composed by ϒ and δ chains. In adults chickens, ϒδ T cells are the dominant population among intestinal epithelial lymphocytes (20 to 60% of the lymphocytes) (73), circulating T cells (20 to 50%), and splenocytes (30%) (74). They are also contributeing to the rst line of defense against pathogens (73). This population grows during infectious challenges with Salmonella enterica (75), and demonstrates in vitro a capacity for cytokine and interferon production (IFNϒ, TNFα, IL17). Their dependency on APC and MHC restriction for antigen recognition remains uncertain, and some ϒδ T cells recognizing directly bacterial and viral epitopes were identi ed (73). In mice intestinal lamina propria, γδ T cells are the major source of IL-17, which is gut protective as it promotes the repair of damaged intestinal epithelium through the regulation of occludin subcellular localization to cell junctions (76,77). Furthermore, some reports in mouse show that lipid antigens from E. coli are recognized by the hepatocytes CD1d, which alleviates hepatic γδ T-17 cell clones (78). In the present study, we observed that several OTU from Eschericha-Shigella genus were negatively correlated with TCR ϒδ expression on TCRγδ T cells, and were also more abundant in the line B12 that was characterized as a low responder to the IBV vaccination. We also observed that Lipid X, IVa, IIa and IIb, Lipid A disaccharide-synthase were predicted to be more abundant in microbiota from the low responder group than from the high responder group. Despite that no orthologous gene of CD1d are known in chicken genome, we can wonder if the CD1 paralogous genes are implied in a similar mechanism (79). Thus, the effect of bacteria from Escherichia-Shigella genus on poultry γδ T cells, and on systemic level of bacterial lipid antigens seems worth of being deeply studied.
A relation of vaccination e ciency with a set of OTUs and predicted functions Some OTUs were signi cantly differentially abundant between groups having high or low vaccine e ciency. Among the most signi cant differences, an OTU from the Tyzzeralla family, from the Angelakisella genus, and to a lower extent an OTU from the Flavonifractor family were more abundant in the HR than in the LR group, while an OTU from the CAG-352 genus was more abundant in the LR than in the HR group. The identi ed OTU from the Angelakisella genus was the only bacterial species observed in this genus, while we identi ed four OTUs from the Tyzerella genus, and 9 from Flavonifractor in the community. It remains di cult to interpret these higher and lower abundances of OTUs. To our knowledge, only few of the identi ed bacteria have already been linked to immunity. Intestinal bacteria from Flavonifractor genus are known to produce butyrate (80), which has anti-in ammatory effects and regulate the production of in ammatory cytokine through its action on intestinal immune cells (81). Flavonifractor bacteria are also polyphenol converting bacterium that degrades avonoids into phenolic acids and 3hydroxyphenylpropionic acid. This last product has been associated with an improved response to In uenza infection through an increase of type I interferon signaling (82). Furthermore, the abundance of Flavonifractor is also associated with a reduction in plasma C-reactive proteins (83), and to higher intestinal endocrine function through GLP-1 secretion (84). In addition to harboring a microbiota dissimilar from HR, the bacterial community from LR was also predicted to have a signi cantly higher abundance of genes involved in lipopolysaccharide (LPS) and peptidoglycan (PG) synthesis. Even if there are still no proof for coronaviruses, LPS and PG are known to enhance the virion stability of poliovirus and reovirus and thus increase their infectivity (85). We can thus wonder if there were different LPS and PG loads between the high and low responder group, that could have subsequently impacted the intestinal IBV viral load, and if there were difference in type I interferon production between the HR and LR animals.

Conclusion
In the present study, performed on inbred chicken lines with contrasted abilities to response to IBV vaccination, we observed during an IBV infection that the different lines harbored dissimilar microbiota. In addition to this role of the host genetics on the microbiota composition, we observed that the microbiota from the vaccinated chickens was highly dissimilar from the microbiota of the control group. We hypothesized that a bi-directional crosstalk may have occurred between immunity and the gut microbiota. Vaccinated animals, that had developed a protective immunity, displayed a different immune response upon challenge which may have dampened the shift of their gut microbiota during the infection. On the contrary, the control animals that encountered IBV for the rst time had a higher immune disturbance, which in conjunction with a reduced appetite and thirst, an increased body temperature, may have fostered their gut microbiota shift. Nevertheless, regardless of vaccine status and the ability to respond to IBV vaccination, we observed that some bacteria from the Escherichia-Shigella genus were always negatively correlated with TCR ϒδ expression on circulating TCRγδ T cells. Those bacteria may affect the γδ T cell differentiation or proliferation, and thus their subsequent cytokine production. Furthermore, we observe that OTU from the Tyzzeralla family, from the Angelakisella and the Flavonifractor genus, and that predicted genes from lipopolysaccharide and peptidoglycan biosynthesis were differentially abundant between the high and low responder to IBV vaccination.