Influence of the digestive/reproductive tract microbiota on chicken egg production beyond host genetics

, Abstract Background: The microbiota of the digestive and reproductive systems has a prominent role in animal health and performance, but the extent of its contribution is difﬁcult to determine. In chickens, the effect of host genetics on the reproductive and digestive tract microbiota is unclear, and the means by which digestive/reproductive microbiomes help improve egg production in chicken are unknown. Results: To gain insight into this, we examined genomes from 128 chickens reared under identical conditions and described their digestive (crop, gizzard and small intestine) and reproductive tract (vagina, uterus and isthmus) microbiota. Although the diversity, composition and predicted function of the digestive and reproductive tract microbiota exhibited notable microbiota variation substantially between different parts, host genetics had limited effects on the reproductive and digestive tract microbial community. The digestive and reproductive tract microbiota had a significant effect on egg production (accounting for 52.31% - 98.86% of the variance), after correcting for host genetic effects; in particular, the uterus and isthmus microbiota accounted for an average of 93.59% and 98.86%, respectively, of variance in egg production. We further identiﬁed four reproductive tract microbial species which were related to immune system, Bacteroides fragilis , Bacteroides salanitronis , Bacteroides barnesiae and Clostridium leptum , that were signiﬁcantly positively correlated with egg production. Chickens with a lower abundance of these species had produced signiﬁcantly fewer eggs at 300 days of age (37.13 vs. 113.75) than those with a higher abundance of these microorganisms. These taxa indicate potential roles play in promoting reproductive performance. Especially uterus and isthmus tract microbiota were major factors in regulating the chicken egg production. Conclusions: Host genetics has limited effect on digestive/reproductive microbiome composition. The distinct site-associated chicken microbiome may be determined by the differences of their physical function. These ﬁndings may help design strategies for controlling and altering the digestive/reproductive tract microbiota in chickens to improve egg production. (LDA > 4) and taxa whose abundances differed significantly among the six sites. phylum; o, NS, significant. g Number of species identified in all chickens at the six different sites.

The chicken (Gallus gallus domesticus) is a commonly raised livestock animal with important economic value, contributing to global food production for humans by providing meat and eggs, as these products have a very high nutritional value and are relatively inexpensive [1]. Improvement of egg laying performance, such as egg production up to 46 weeks of age, egg weight and age at first laying, is urgently needed in the poultry industry to meet the tremendous global demand [2]. With the modern revolution in advanced technologies, genome-wide association studies (GWASs) have been widely used to identify single nucleotide polymorphism (SNP) genotypes and candidate genes associated with quantitative traits in birds [3,4]. To date, in egg production and egg quality traits of chickens, about 3,135 quantitative trait loci (QTLs) have been reported, among which 315 QTLs on 19 different chromosomes (1-11, 13, 17, 19-20, 22-24 autosomes and Z sex chromosome) were associated with egg number and are listed in the Animal QTLdb (https://www.animalgenome.org/cgi-bin/QTLdb/GG) [5]. There is plenty of GWAS evidence showing that promising mutants located within genes, such as APOO, RGS1, PRDX4, and CLSPN, are associated with reproductive traits [6][7][8].
Egg production is an important economic trait in the poultry industry [9], and it is a polygenic inheritance trait with low to moderate heritability values ranging from 0.05 to 0.44, depending on the period involved [10,11]. Therefore, other effective approaches for modulating egg production in laying hens are urgently required. Among them, the role of the microbiota is a relatively new field of study that is gaining importance. Egg production may be regulated by the gut and reproductive tract microbiota. In humans, some studies have been conducted to explore the relationship between microbes and reproduction. For example, the dynamics of microbiome interactions with the host during pregnancy leading to preterm birth were investigated [12] and temporal changes in the vaginal microbiome associated with full-term pregnancies were identified [13]. Manipulation of the gut and reproductive tract microbiome can alter the outcome of hormone-related pathologies, including pregnancy and infertility [14][15][16]. An abnormal vaginal microbiota may predispose individuals to ascending colonization of the genital tract, microbial invasion of the amniotic cavity and fetal damage [17,18].
The avian digestive and reproductive tracts house complex bacterial communities that are believed to play crucial roles in the overall health and egg production traits of birds [19][20][21][22]. Most of the studies have mainly focused on broilers to measure the composition of commensal and pathogenic bacteria in the chicken gut, likely due to the commercial importance of these birds in the meat industry [23,24]. It has been identified that Methanobrevibacter and Mucispirillum schaedleri play crucial roles in chicken fat accumulation [25]. Typically, the chicken gut and reproductive tracts are colonized by bacterial phyla mainly Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria and Fusobacteria [26,27], that are spatially organized within specific gut/reproductive compartments. In addition, Lactobacillus species were found to be as the keystone species residing in the chicken oviduct. The functions of these microorganisms in the reproductive system and their effects on egg production are worthy of study.
More recently, a metagenomic association study showed that the gut microbial communities of chickens with high egg production performance differ significantly from those of chickens with low egg production performance [28].
Several synergistic factors, such as the environment and diet dominate over host genetics in determining gut microbiota composition [29][30][31]. A comparative study of gut microbial diversity among parrot species, indicated the potential role of host ancestry in shaping gut microbiomes [32]. In birds, genetics specialists applied GWASs to demonstrate the relationship between host genotype and gut/reproductive microbiota population [33][34][35][36][37][38]. In a previous study, 14 QTLs were identified as exerting a strong effect on C. leptum and Lactobacillus and on related candidate genes involved in anti-inflammatory responses and in the motility of the digestive system [39]. On the other hand, recent studies have suggested that host genetics have a limited impact on the gut microbiota composition of humans [40] and chickens [25], and microorganisms play crucial roles in chicken fat accumulation [25], and microbes have important effects on host traits independent of host genetics. Therefore, we objectively realized that the colonization by some microorganisms is influenced by host genetics and that the microbiota of various sites can exert a great impact on different phenotypes of the host.
As suggested above, we speculated that the microbial complement of the digestive-reproductive system might be an important aspect of chicken egg production. To what extent host genetics affect the microbes of the digestive and reproductive system and the influence of these microbes on egg production are unknown. In this study, we performed whole-genome sequencing on 128 green-shell laying hens and 16S rRNA gene sequencing on 768 samples from digestive (crop, gizzard, small intestine) and reproductive (vagina, uterus, isthmus) sites (i) to characterize the digestive and reproductive tract microbiota of hens, (ii) to establish a correlation between host genetics and microbial diversity of the digestive/reproductive tract, and (iii) to identify the contribution of key microorganisms to egg production. The findings of this study will provide insights into the microbial communities in different digestive and reproductive tracts and their independent contributions to egg production.

Experimental Animals and Sample collection
The study was conducted on a flock of 128 green shell laying chickens reared on an experimental poultry farm at Sichuan Agricultural University in Ya'an, China. All birds were hatched on the same day and housed in individual pens. Feed intake was controlled daily according to standard farm husbandry practices and water was provided ad libitum. The number of eggs produced for the first 300 days of life was recorded daily for each bird. At the age of 300 days, 2 mL of whole blood was collected from the wing vein using venipuncture and stored at -20°C. Subsequently, each bird was culled by cervical dislocation followed by decapitation. After laparotomy, one centimeter of fresh tissues were collected from different reproductive tract (vagina, uterus and isthmus) and digestive tract sites (crop, gizzard and small intestine), as shown in Fig. 1 Finally, low-quality reads (> 50% bases with phred quality < 5) were also removed. Consequently, we generated 57,609,211 high-quality reads for further analysis (Additional file 1: Table S1). The quality of clean data was assessed by the length distribution of reads, quality distribution, and error rate distribution. More than 93.85% of the high-quality reads had lengths of 250 -260 nt (Additional file 2: Figure S1a). Data with a quality score greater than 20 accounted for 88.14% of all the effective bases (Additional file 2: Figures S1b, c). The error rate of the sequencing reads was relatively high in the ending position (Additional file 2: Figure S1d).

OTU cluster and species annotation
The remaining high-quality sequences were used to generate OTUs (operational taxonomic units) by

Alpha Diversity
Alpha diversity is used to analyze the complexity of species diversity for a sample based on the normalized OTUs through 6 indices (i.e., observed OTUs, Chao1, Shannon, Simpson, ACE, Good's coverage), using the package QIIME2 (Quantitative Insights into Microbial Ecology). Of these indices, two were selected to identify community richness (i.e., Chao1 and ACE); two were used to identify community diversity (i.e., Shannon and Simpson); and was used to characterize sequencing depth (i.e., Good's coverage). The differences in alpha diversity indices among the six tracts were calculated with the Wilcoxon rank-sum test using R software (Ver 2.15.3).

Beta Diversity
Beta diversity was used to evaluate differences in samples. Beta diversity in Bray-Curtis (BC), weighted and unweighted UniFrac distances was calculated by QIIME software. The Bray-Curtis ordination provides position values along an ordination axis and distances from the axis for samples of Cluster analysis was preceded by principal component analysis (PCA), which was applied to reduce the dimension of the original variables using the FactoMineR package and ggplot2 package in R software.
Principal coordinate analysis (PCoA) was performed to obtain principal coordinates and visualize complex, multidimensional data. A distance matrix of weighted or unweighted UniFrac distances among samples obtained previously was transformed to a new set of orthogonal axes, by which the maximum variation factor was demonstrated by the first principal coordinate, the second maximum variation factor was demonstrated by the second principal coordinate, and so on. PCoA was performed using the WGCNA package, stat packages and ggplot2 package in R software.
Unweighted Pair-group Method with Arithmetic Means (UPGMA) Clustering was performed as a type of hierarchical clustering method to interpret the distance matrix using average linkage and was conducted by QIIME software.

Prediction of the functional profiles of microbial communities
The functions of the microorganisms present in the microbial communities detected at each of the six anatomical sites (uterus, isthmus, vagina, gizzard, crop and small intestine) were predicted using PICRUSt2 [49]. We used the Wilcoxon rank-sum test to investigate the differences in pathways among the six tracts. These p-values were adjusted using the Benjamini-Hochberg method by the false discovery rate (FDR) with the p.adjust function in R.

Community differences analysis
The different sites were statistically compared using ANOSIM (also named permutational MANOVA) [50] with 10,000 permutations based on Bray-Curtis (BC) ordination to evaluate the reasonability of the division of groups.

Between-group variation analysis
The high-dimensional biomarkers were discovered by LEfSe (linear discriminant analysis (LDA) Effect Size) with the parameter 'LDA score > 4' [51], for identifying characteristics of abundance and related classes (e.g., genes, metabolites or taxa).

Whole-genome sequencing and quality filtering
After qualified host DNA samples were tested, a Covaris ultrasonic fragmentation instrument was used to randomly fragment the DNA, and then a sequencing library was prepared using a TruSeq Nano DNA bases having phred quality < 5; and putative PCR duplicates generated in the library construction process), which mainly resulted from base-calling duplicates and adaptor contamination, were removed using in house script [52]. Consequently, we obtained 1.30 Tb (~10.15-fold per individual) of high-quality paired-end reads, including 95.13% and 88.98% nucleotides with phred quality ≥ Q20 (with an accuracy of 99.0%) and ≥Q30 (with an accuracy of 99.9%), respectively (Additional file 1: Table S2).

Read mapping, genomic variant calling and annotation
The remaining high-quality reads of each birds were aligned to the chicken reference genome (Gallus_gallus-6.0 Ensembl release 94, http://asia.ensembl.org/Gallus_gallus/Info/Index/) using the Burrows-Wheeler Alignment tool (BWA ver 0.7.15) [53] with the command 'mem -t 10 -k 32'. BAM alignment files were then generated using SAMtools (ver 0.1.19). Additionally, we improved the alignment performance through filtering the alignment reads with mismatches ≤ 5 and mapping quality = 0. After sorting by SAMtools, the sorted bam file was marked in duplicate using the command "MarkDuplicates" in the package picard (ver 1.119).
Subsequently, we performed gVCF calling according to the best practices using the Genome Analysis Toolkit (GATK, ver v3.7) [54] with the HaplotypeCaller-based method, and then population SNP calling by merging all gVCFs with the commands "CombineGVCFs".
To obtain high-credibility SNPs, we applied hard filter command 'VariantFiltration' to exclude variants were further filtered when more than three SNPs clustered within a 10-bp window were removed [55]. Finally, we used vcftools (ver 0.1.15) [56] to obtain biallelic variants with the following parameters: sample call rate > 90%, SNP call rate > 95%, minor allele frequencies > 1% and Hardy-Weinberg equilibrium p-value < 10 −5 [57]. Ultimately, a total of 10,819,101 high-credibility SNPs in 128 chickens were retained (Additional file 1: Table S3). We further annotated SNPs using the package ANNOVAR [58], the SNPs were classified into different genomic regions (i.e., exonic, intronic, splice sites, upstream and downstream around gene regions and intergenic).

Construction of microbial relationship and host genetic relatedness matrices
The OTUs identified at each anatomical site were normalized to a zero mean and unit variance. We then constructed a microbial relationship matrix (MRM) [36] using an R script based on the following equation: where represents the tested microbial relationship in tract t between chickens i and j; and are the abundance of OTU o in tract t in chickens i and j, respectively; ̅̅̅̅ is the average relative abundance of OTU o in tract t in the population; 2 is the variance in the abundance of OTU o in tract t; and N T is the total number of OTUs in tract t used for the computation of relatedness. The high-quality SNPs were further used to detect independent markers using PLINK [59], with the following parameters: 50 kb window size, 10 SNPs per step and 0.2 as a squared Pearson's correlation coefficient ( 2 ). All 10,809,968 SNPs were used to compute the principal components (PCs) and GRM [60] using GCTA (ver 1.91.1) [61]: where ℎ is the tested genetic relationship between chickens i and j; and represent the number of reference alleles in chickens i and j, respectively; the ̅ is the frequency of the reference allele in the population; and N is the number of variants.

Heritability analysis
To estimate the effects of host genetics on the digestive/reproductive tract microbiota at different anatomical sites, we first computed the correlation between the GRM and BC distances at each site using both Pearson's and Spearman's correlation-based Mantel tests with 10,000 permutations. The correlation between the GRM and the MRM was also computed. To estimate the correlation between the GRM and the reproductive/digestive tract microbiota community, we computed the heritability at the OTU, genus and species levels. OTU abundance information was normalized using a standard sequence number corresponding to the sample with the least sequences. The microorganisms that were present in < 60% but ≥ 20% of the samples were dichotomized as present or absent [62], and the microorganisms that were detected in < 20% of the samples from each site were excluded from the analysis.

Genetic and microbial parameters of egg production
As the chickens examined in this study had no pedigree information, we computed the SNP-based heritability of the egg production phenotype (EN300) instead, using the following model: where y is an observed value (EN300); c is a vector of fixed covariates with the corresponding design matrix K; e is the residual effect; and g is a vector of aggregate effects of all SNPs with an ~N(0, 2 ), where G and 2 are the GRM and polygenetic variance (overall SNP effects), respectively. The top five host genetic PCs were considered covariates in the model to account for the calculated population stratification, as described above. The likelihood ratio test p-value was calculated to examine the significance of the association between SNPs and EN300.
The fraction of the EN300 variance explained by the microbial variance was calculated as 2 = with a univariate linear mixed model (LMM), which was performed using GEMMA [63]. The LMM was calculated as follows: where the model parameters are the same as those described in model [A], except , which is the random effect of the digestive/reproductive tract microbiota in location s following the multinomial distribution ~N(0,M 2 )), and M is the MRM. We then used the MRM in GCTA to calculate 2 .
The genome-wide significance threshold was 10e-6. We then extracted these SNPs with significant effects on EN300 and calculated the PCs using PLINK. The first two PCs and the top five host genetic PCs were then used as covariates in model [B] to account for host genetics.
Identification of microbiota constituents related to egg production EN300 between two groups (the 20% lowest-and 20% highest-ranked chickens regarding EN300) was then compared using the Wilcoxon rank-sum test. Microorganisms with P < 0.05 and FDR < 0.05 were retained. Furthermore, we calculated the Spearman's and Pearson's correlations between EN300 and the abundance of each microbiota constituent at the genera, OTU and species levels. The presence of a microorganism was considered to correlate with EN300 if P < 0.05, as determined using the psych package in R with the p-value adjusted using the BH method. The overlapping microorganisms obtained from the Wilcoxon rank-sum test and Spearman's and Pearson's correlations were considered to have a potential relationship with EN300. We subsequently characterized the amount of each EN300-associated microbes at each site.

Diversity and composition of the reproductive and digestive tract microbiota
To describe the microbial landscape of the chicken digestive and reproductive tracts, we profiled the  Table S1). De novo clustering after singleton removal produced 46,480 operational taxonomic units (OTUs) at an identity cutoff of 97%. We then discarded OTUs that were not found in at least 20% of these samples, resulting in a total of 6,776 OTUs.
For alpha diversity analysis, the Good's coverage ranged from 96.3% to 99.6% with an average of 98.7%, indicating that the measurement depth met the requirements (Additional file 2: Figure S2a).
Analysis of the other five indices (observed OTUs, Chao1, Shannon, Simpson and ACE indices) indicated that all pairwise comparisons showed significant differences (P < 0.001) except for two comparisons (uterus vs. isthmus was not significant for all indices; and small intestine vs. vagina was not significant for observed OTUs, Shannon and Simpson indices) (Additional file 2: Figures S2b-g); the microbiota alpha diversity in the uterus and isthmus was higher than that in other sites, indicating that these sites contained more microbial taxa than the vagina and the digestive tracts (  Only genera with an abundance of >1.0% in any of the sites are described (Additional file 1: Table S5).
Although similar phyla dominated the microbiota at all six sites (Firmicutes, Proteobacteria and Cyanobacteria, accounting for 71.45% to 97.86% of all OTUs), there were some noticeable differences (Additional file 2: Figure S2g) Figure S2g). Interestingly, the vagina had the highest abundance of Fusobacteria (11.51%) among these six sites.
At the genus level, the crop (73.74%) and small intestine (30.70%) were dominated by Lactobacillus, and the crop exhibited the lowest alpha diversity (Fig. 1a). Gallibacterium (2.14%) was found only in the crop. Lactobacillus (24.76%) was found abundant in the gizzard. In the reproductive tract samples, Lactobacillus (7.24%-9.27%) were also dominant, and bacteria such as Exiguobacterium  Table S5). Furthermore, for the majority of genus, the relative abundance at one site was not associated with that at other sites. Among these 2,475 Spearman correlation coefficient, 2,113 (85.37%) of them were not significantly correlated (Additional file 2: Figure S3). A limited number of genera belonging to the Proteobacteria and Firmicutes phyla showed significant (P < 0.001) and positive correlations between the crop and gizzard, the gizzard and small intestine or the three reproductive tract sites.

Predicted functions of the microbiomes and site-associated microorganism identification
To determine whether the microbiomes identified as described above have specific functions within their distinct microenvironmental niches. We analyzed the functional capacity of the reproductive and digestive tract microbiota using PICRUSt2. Our results revealed that the top 50 KEGG pathways enriched at each site were mostly the same (shared 36 pathways), and twelve (33.33%) of the 36 shared pathways were associated primarily with metabolism (Additional file 2: Figure S4). Seven pathways were identified that were only enriched at a single site: 3 and 3 pathways for the crop and gizzard respectively, 'Riboflavin metabolism' for the vagina. The abundance of the OTUs involved in these pathways differed significantly (adjusted p-values < 0.001) among the six sites (Additional file 1: Table   S6). For example, the crop microbiota was significantly enriched in 'D-Alanine metabolism' and 'Nucleotide excision repair', the gizzard microbial community had important roles in 'Sulfur relay system' and 'Starch and sucrose metabolism', and the small intestinal microbial community had an important role in 'Valine, leucine and isoleucine biosynthesis', as indicated by the moderate row z-scores for each pathway. Moreover, 'Propanoate metabolism' and 'Bacterial chemotaxis' were overrepresented in the vagina. Meanwhile, 'Bacterial secretion system' was overrepresented only in the uterus and isthmus (Fig. 2a).
We then used LEfSe to explore site-associated bacterial features (Additional file 2: Figure S5), and the relative abundances of these features was visualized on a heat map (Fig. 2b). LEfSe analysis confirmed most of the observations described above. For example, the isthmus and uterus had similar abundances of similar microorganisms; Lactobacillus and Gallibacterium were classified as crop-associated bacteria, whereas unidentified_Chloroplast was a gizzard-associated bacterium. Of note, Helicobacter and unidentified Erysipelotrichaceae, which were associated with the small intestine, showed the highest abundance among the six sites (Fig.  2b).
Unidentified_Erysipelotrichaceae also showed high abundance in three reproductive sites (1.83%-2.40%, Additional file 1: Table S5). Romboutsia, Fusobacterium, Turicibacter and Clostridium_sensu_stricto_1 were vagina-associated bacteria. The bacteria associated with isthmus and uterus both showed higher abundance at those sites than in other sites.  Table S6). The heatmap is color coded based on row z-scores of means of KEGG pathway abundances.
b Heat map showing the 65 site-associated bacterial taxa identified by LEfSe (LDA > 4) and taxa whose abundances differed significantly among the six sites. p, phylum; c, class; o, order; f, family; g, genus.

Weak association between host genetics and digestive/reproductive tract microbial communities
Next we sought to determine the association between the microbiota and genetics in the same cohort of laying hens. To do this, we first generated 1.76 tera bases (Tb) of high quality sequences from 128 chickens (Additional file 1: Table S2), of which the quality of 96.39% of the bases was more than Q20.
Mapping to the chicken reference genome resulted in an average of 99.81% reads covering an average of 95.59% of the nongap genome with a ~10.15-fold average depth. A total of 10.82 Mb of highly credible single nucleotide polymorphisms (SNPs) with a density of 10.29 SNPs/kb were acquired from the 128 individuals. We then calculated the correlation between host genetics (using genetic relatedness matrix) and microbial beta diversity (based on Bray-Curtis (BC) distance) at all six sites. The results showed that there was no significant association (|r|< 0.033, P>0.05, Fig. 3a-f and Additional file 1: Table S7) between host genetics and microbiome composition. However, the correlation coefficients for adjacent anatomical sites of BC distances were larger than those for nonadjacent sites, indicating that microbiomes in adjacent sites were similar; for example, a positive correlation was detected between the uterus and isthmus (r=0.4260, P<0.0001, Spearman correlation). No significant association (r=0.0186, Spearman correlation) was found between the crop and isthmus microbial communities (Additional file 1: Table S7). We then estimated the association between the genetic relatedness matrix (GRM) and microbial relationship matrix (MRM), and obtained similar results: both Pearson's and Spearman's correlation coefficients of the GRM with any of the six MRMs were close to zero (Additional file 1: Table S7), suggesting that there is no connection between genetics and the microbiota.
Then, we regarded each microorganism as a quantitative trait (abundance as a phenotype) and computed the heritability of each microorganism at the species, genus and OTU levels. The microorganisms in more than 20% but less than 60% of the samples were thus analyzed qualitatively as dichotomous traits; microorganisms present in less than 20% of samples were excluded from the analysis (Fig. 3g, Additional file 2: Figure S6a). At the species and genus level, there was no statistical correlation (P < 0.05) between the presence of a SNP and the presence of a specific microbe (Additional file 2: Figure S6b Table S8). Most of these heritable bacteria detected at all three levels belonged to the Firmicutes phylum (Additional file 2: Figure S7a, c, e). Reproductive segments had more heritable bacterial phyla than digestive segments (Additional file 2: Figure S7b, d, f). The cumulative abundances of these heritable bacteria were only 0.22%, 4.14%, 1.46% and 1.61% (P < 0.05) in the small intestine, vagina, uterus and isthmus, respectively (Fig. 3h). Similar results were obtained at the genus and OTU levels (Additional file 2: Figure S6c). This result was consistent with the results described above showing that host genetics have a limited effect on shaping the reproductive and digestive tract microbial composition. h Cumulative relative abundance of heritable microbial species in each segment.

Heritability and microbiability of egg production
As there was no evidence of an association between host genetics and the digestive/reproductive tract microbiota, we next asked whether host genetics and the microbiota at different reproductive/digestive sites affect egg production. To do this, we determined the mean number of eggs that each hen had laid by 300 days of age (EN300). The mean EN300 was 75.32 (range: 24-129) (Additional file 2: Figure   S8), and the trait fit a normal distribution (Shapiro-Wilk normality test, p>0.05). A genetic relationship matrix (GRM) of pairs of samples was used as input for the restricted maximum likelihood analysis to estimate the heritability of EN300 explained by the whole genome SNPs. As a result, we found that the EN300 exhibited relatively low heritability (h 2 = 0.28) (Additional file 1: Table S9). The heritability of EN300 was within the scope of previous reports (Additional file 1: Table S10). This suggests that host genetics have a limited effect on egg production in chickens.
Then, we performed a genome-wide association study (GWAS) and identified 11 EN300 related SNPs ( Fig. 4a and Additional file 1: Table S11), five and four of which were located on chromosomes 18 and 3, respectively. All 11 SNPs were noncoding mutations: 7 were intergenic, and 4 were intronic.
One of the intronic SNPs was located in GRM5 (glutamate receptor, metabotropic 5), and three were located in ROCK2 (rho-associated, coiled-coil containing protein kinase 2). The lower left corner of The fraction of the EN300 variance explained by the microbial variance was calculated as microbiability (m 2 ). We used all 11 SNPs as additional covariates to correct for host genetic factors in the m 2 estimation. The estimated EN300 m 2 for the reproductive sites (0.923 for the vagina, 0.936 for the uterus and 0.989 for the isthmus) was higher than that for the digestive sites (0.867 for the crop, 0.873 for the gizzard and 0.523 for the small intestine) (Additional file 1: Table S12). Taken together, these results suggest that egg production chicken is primarily determined by the digestive and reproductive tract microbiota independent of host genetics.

The microorganisms significantly associated with egg production
We then analyzed how well egg production can be inferred from the digestive/reproductive tract microbial communities in comparison with that from host genetics. A double-divergent Wilcoxon rank-sum test was conducted to identify egg production-related microbial species, genera and OTUs (p<0.05). Most of the microorganisms detected at the three different levels that were significantly associated with egg production belonged to the Firmicutes phylum (Additional file 2: Figure S9). Only microorganisms that exhibited a significant correlation of egg production with the relative abundance as determined by both the Pearson correlation coefficient (PCC) and Spearman correlation coefficient (SCC) were considered causal (p<0.05). In total, 39 OTUs, 26 genera and 24 species fulfilled these criteria (Fig. 5a, Additional file 2: Figure S10a). In addition, we identified candidate genera related to egg production by LEfSe using downloaded data for commercial laying hens with high egg production rates and local chicken breeds with low egg production rates at similar ages (~50 weeks). The results showed that the microorganisms associated with egg production overlapped with those identified in this study (all five genera related to egg production in feces were detected in our study, and 4 of 7 genera related to egg production in the cecum were candidate genera in our study) (Fig S10b). Most OTUS, genus and species present at the three digestive tract sites were negatively correlated (using PCC) with egg production (negative/positive: 19/6, 8/8 and 16/5 respectively), whereas most taxa at the three reproductive tract sites were positively correlated with egg production (positive/negative: 13/9, 11/10 and 13/4 respectively) (Fig. 5b, Additional file 2: Figure S10c, d). The microorganisms in the uterus were most strongly correlated with each other (Additional file 2: Figure S11), which implied a strong symbiotic/competitive relationship.  Figure S12). The 20% of birds with the highest Firmicutes bacterium ZOR0006 abundance exhibited a significantly higher EN300 than the 20% of birds with the lowest abundance of this microorganism (Fig. 6a) (P < 0.05) at reproductive sites.
Moreover, the 20% of birds with the highest Bacteroides fragilis, Bacteroides salanitronis, Bacteroides barnesiae and Clostridium leptum abundances showed a significantly lower EN300 than the 20% of birds with the lowest abundance in reproductive segments, with the exception of Clostridium leptum abundance in the vagina (Fig. 6a).
We then characterized the spatial distribution of these five EN300-associated microorganisms (Bacteroides fragilis, Bacteroides salanitronis, Bacteroides barnesiae, Clostridium leptum and Firmicutes bacterium ZOR0006) at multiple segments. Bacteroides fragilis was detected in almost all samples and accounted for 0.05% to 1.29% of the total abundance (Fig. 6b). Bacteroides salanitronis and Bacteroides barnesiae were detected at similar rates at the six sites and in most reproductive site samples; both accounted for the highest abundance in the vagina. Firmicutes bacterium ZOR0006 was also detected in most reproductive segment samples (74.22%-89.8%) and in half of the digestive site samples (48.44%-58.59%), accounting for 0.61% to 2.40% of the total abundance. The detection rates (28.13%-64.84%) and relative abundance of Clostridium leptum were much lower than those of the other microorganisms at all six sites, but accounted for the highest abundance in the isthmus and uterus (Fig. 6b).

Alpha diversity
In order to characterize the digestive and reproductive tracts microbiota of hens, we demonstrated the spatial heterogeneity of the microbiota throughout contiguous segments of the digestive and reproductive tracts in chickens (including the crop, gizzard, small intestine, vagina, uterus and isthmus).
Compared to the digestive system, the reproductive system, especially the upper reproductive tract were no differences between the isthmus and uterus. A recent study also showed no significant differences in the abundance of four microbes composition among different sections of the oviduct (infundibulum, magnum, isthmus, and uterus) [66]. Nonetheless, chicken digestive tract microbiome diversity was comparable with those of the reproductive tract, suggesting that the different function of chicken digestive and reproductive tract microbiome.

The chicken digestive and reproductive tract microbiome composition
This study provides a comprehensive view of the succession of the chicken digestive and reproductive tract microbiome. Such a study design allowed us to address the logical question regarding the chicken microbiome: How does the chicken microbiome change across different digestive and reproductive tracts?
Lactobacillus and Exiguobacterium were found at all sites due to their broad adaptability or beneficial functions. Lactobacillus are thought to play key protective roles by lowering the environmental pH through lactic acid production [67], so the high abundance of this genus in the digestive tracts corresponds to the low pH value of this organ system [68]. Lactobacillusor other lactic acid bacteria-dominant communities were detected less abundant in the uterus and isthmus (~14.03%) and the relative abundance was lower than that of the mature human vaginal microbiota (~50%) [69]. Consistent with a recent study in chickens [66], this finding may be related to the alkaline pH value in the chicken reproductive tract that maintains sperm motility [70,71], in contrast to the low pH in humans that prevents pathogen infection and growth [72][73][74]. Exiguobacterium is a genus of bacilli and a member of the low-GC phylum Firmicutes; it has been found in areas covering a wide range of temperatures, pH values, salinity and other conditions [75], and it also produces organic acids and hydrolytic enzymes [76], even having some potential function in the degradation of toxic substances [77]. Given that the reproductive tract and digestive tract share a common exit in the cloaca, we propose that microbial exchange occurs at this site, resulting in similar microbiotas at the distal end of both tracts. Thus, we found that the vagina as a reproductive organ that acquires microbe communities from the isthmus and uterus but also has diversity similar to that in the small intestine ( Fig. 1b). A functional capacity analysis of the six sites microbial communities revealed that the top 50 pathways enriched in each segment were mostly the same (shared/total: 36/50) and primarily associated with metabolism. Specifically, 'bacterial secretion system' and 'bacterial chemotaxis' are over represented in the reproductive tract (vagina, uterus and isthmus). Previous study found that successful bacterial pathogens have evolved versatile protein secretion systems that they use to promote their survival and fitness in response to different environmental challenges, and also to modulate host immunity [78].
Our findings confirm that the chicken digestive and reproductive tract microbiota is primarily determined by the physiological function of each compartment within these systems.

Site-associated microorganisms in chicken digestive and reproductive tract
This study also enabled us to address some other important biological questions regarding the chicken digestive and reproductive tract microbiome, such as: What are the site-associated microorganisms in diverse sites? Site-associated bacterial features in this study were identified by using LEfSe, a total of 65 site-associated bacterial taxa were identified among six sites.
Six Lactobacillaceae bacterial taxa were detected as crop-associated. In the chicken digestive system, the crop acts as a reservoir for the storage of food prior to its digestion, where food mixes with many beneficial Lactobacillus bacteria (73.84% at the genus level) that produce lactic acid before moving on to the proventriculus [79,80]. Next, a strong muscle known as the gizzard functions similar to teeth for chickens by grinding any remaining large food particles with the assistance of grit (upon kneading), releasing abundant unidentified_Chloroplast (2.94%) and mitochondria-like (2.37%) microbes from plant consumption. The small intestine had the most abundant microbes in the three digestive sites, which is the main area where further digestion occurs and fermentation begins. The overrepresented Paenibacillaceae species were found with an optimum growth at pH 6.0-7.0 [81].
Helicobacter specifically inhabits the chicken small intestine to maintain near neutral pH and a microaerophilic environment and may be involved in inflammation, metabolism and neutralization of gastric acid [82][83][84]; whether it is beneficial or harmful remains unknown.
In this study, we observed that some representative genera in the vagina, such as Romboutsia, Fusobacterium and Clostridium_sensu_stricto_1, accounted for more than 25% of the vaginal microbiota, which were in lower abundance in other sites, although the role of most of these dominant bacteria in vagina remains unknown. Previous study found that genus Romboutsia was negatively associated with the body weight, fasting serum glucose and insulin. Fusobacterium has been reckoned as a proinflammatory organism [85]. It has been observed that the genus Clostridium is capable of utilizing a broad range of organic substrates for efficient hydrogen production [86] and thus adjusting pH value of host. Gammaproteobacteria and Betaproteobacteria were detected as uterus-associated.
Significant changes in these bacterial taxa (Gammaproteobacteria, Betaproteobacteria) were detected in immune-suppressed adults emerging from Varroa -infested colonies [87]. Six Bacteroidetes bacterial taxa were detected as isthmus-associated, Bacteroides species are thought to play a vital role in the breakdown of polysaccharides into simpler compounds that related to host immunity [88]. Acidic pH [68] and alkaline pH [70,71] values were observed in the chicken digestive and reproductive sites respectively. Thus, we speculate that the reproductive microbiota plays a vital role in improving protective immunity, defending against invasion of the reproductive tract by harmful bacterial and adaptation to pH levels.
The predominant microbes identified at the different anatomical sites in this study largely correlated with the physiological functions carried out at those sites.

Host genetics has limited effect on chicken digestive/reproductive tract microbes
Recent studies have suggested that host genetics have an impact on the microbiota composition [25,34,89]. However, the extent to which host genetics shape the digestive/reproductive tract microbiota remains debatable. In order to explore the effect of host genetics on chicken microbes. We first performed a GWAS to obtain EN300-associated SNPs and introduced these SNPs as additional covariates in the model [B] to adjust the effects of host genetics. Our estimation of SNP-based heritability for EN300 was moderate (0.28), and the value was comparable to that obtained in previous estimates [10,11]. We then identified SNPs in two genes (GRM5 and ROCK2) that appear to be strongly associated with egg production. GRM5 was found to be potentially involved in dermatological diseases/conditions in chickens [90], but the function of ROCK2 is unknown. Further studies need to be conducted to explore the role of these two genes in regulating chicken egg production. The closer it was to the ovary, the larger m 2 value of the segment was observed, indicating that the microbiome of the isthmus was the most informative site with respect to egg production. Commercial egg producers have special interest in hen oviducts because pathological changes or disrupted activity directly affect egg production efficiency and ultimately decrease economic profitability [92]. A previous study demonstrated that vaginal microbiota of sows could affect endometritis [93], with endometritis, the proportion of the smallest follicles was only 60% [94]. The isthmus in chicken forms both inner and outer shell membranes, while calcification of the eggshell and its subsequent pigmentation and cuticle deposition occur in the uterus, followed by expulsion of the egg through the vagina [95]. Functionally, it is likely that the microbiota in the uterus and isthmus may affect the health of the oviduct and thus influence chicken egg production, seen as an increased abundance of pathways for bacterial motility proteins, the bacterial secretion system, and membrane and intracellular structural molecules. Thus, we found that the microbiome could influence egg production in chickens, especially the microbiome in the reproductive tract.
Egg production associated microbial species found in the chicken reproductive tract were related to immune function Our finding showed that the reproductive/digestive microbiota was significantly associated with egg production led us to investigate which microorganisms play key roles in this association. At all three levels (OTU, species and genus), we found that the isthmus and uterus contained more microorganisms significantly associated with EN300 than other sites. We observed that the highly abundant species Bacteroides fragilis, Bacteroides salanitronis, Bacteroides barnesiae and Clostridium leptum were significantly associated with a higher EN300 in the three reproductive sites. Meanwhile, a high abundance of Firmicutes bacterium ZOR0006 was significantly associated with a low EN300 in three reproductive sites and a high EN300 in three digestive sites. Among these five EN300-associated species, only one (Bacteroides barnesiae, h 2 =0.43) was affected by host genetics.
In chickens, the genus Bacteroides is found in the caecum and represents one of the predominant anaerobic genera [96]. Bacteroides species are thought to play an important role in the breakdown of polysaccharides into simpler compounds that are used by the animal host as well as the microorganisms themselves [88]. Intestinal anaerobic bacteria such as Bacteroides fragilis and Bacteroides salanitronis have been suggested to possess metabolic pathways for N-glycan production [97]. The symbiont Bacteroides fragilis exists in a commensal relationship with the host as it expresses a relatively large number of genes involved in polysaccharide metabolism, which benefits the host [98]. The surface of Bacteroides fragilis can produce polysaccharides; in particular, capsule capsular polysaccharide A (CPSA) has been shown to be a key mediator of mammalian immune system development [99,100]. Surprisingly, CPSA has also been shown to have protective effects on autoimmune disorder models, such as antibiotic-induced experimental encephalomyelitis [101]. It is thus likely that genus Bacteroides can regulate bird reproductive activity by mediating the avian immune system.
Clostridium leptume a major member of the Firmicutes phylum, can alter the gut microbiota in rats, especially in obese individuals [102]. In human infants, fecal levels of Clostridium leptum were found to be negatively correlated with proinflammatory marker levels [103]. Colonic colonization of Clostridium leptum is associated with the accumulation of regulatory T cells, which inhibits the development of inflammatory lesions [104]. The proliferation and activation of regulatory T cells is crucial to establishing and maintaining an appropriate level of immune tolerance. In addition, our results demonstrated that Clostridium leptum had a large range of associations with other uterus or isthmus microbiota constituents (but limited association was observed with digestive microbiota constituents) and was not influenced by host genetic regulation. It is thus likely that this microorganism serves as a stimulator of regulatory T cell production and inhibitor of inflammatory lesions and then could regulate and maintain immunologic tolerance and the microbiota composition of the reproductive tract (especially the uterus and isthmus). The function of Firmicutes bacterium ZOR0006 is unknown, the abundance of this species was negatively correlated with EN300. In summary, microbial species found in the chicken digestive/reproductive tract that are most strongly associated with egg production probably exert this effect by modulating immune function, and further studies need to be conducted to explore their function in detail.

Conclusions
In summary, our results indicated that the reproductive/digestive tract microbial community of chickens is largely independent of host genetics. However, the microbiota was highly adapted to the environment of a specific site. A small proportion of the variability in egg production was associated with the reproductive/digestive tract microbiota in chickens. In particular, the genus Bacteroides and the species Clostridium leptum and Firmicutes bacterium ZOR0006 were strongly associated with EN300. These findings provide new insight into the roles of the reproductive/digestive tract microbiota in complex traits, and may help contribute to the development of effective therapies for improving egg production.

Availability of supporting data
The sequencing data for this project have been deposited in the National Genomics Data Center under accession number CRA002196 (microbiome) and CRA002195 (Genome resequencing). The authors declare that all other relevant data supporting the findings of the study are available in this article and its Supplementary Information files, or from the corresponding author upon request.

Ethics approval and consent to participate
Chickens were managed according to the Institutional Animal Care and Use Committee of Sichuan Agricultural University under permit number DKY-2018102015.

Consent for publication
Not applicable