Effects of the DMRT1 genotype on the body weight and gut microbiota in the broiler chicken.

Intestinal microbiota is a critical determinant of growth and risk of metabolic diseases. Our previous studies showed that the locus rs16775833 within the DMRT1 gene is significantly associated with variation in the population structure of the gut microbiota, which is involved in determining the BW of the chicken. To assess the accuracy of correlation of rs16775833 located in the DMRT1 gene on microbial population and BW in birds, 2 genotypes GG and TT in the rs16775833 were identified in Chinese Yellow broiler breeders. We found that BW in the TT genotype group was significantly higher than in the GG genotype group at 7 and 13 wk of age in 777 female chickens. A full-length 16S rRNA sequencing approach was used to further evaluate the fecal bacterial composition of female broilers in 11 TT genotype chickens with high weight (HW-TT) and 11 GG genotype chickens with low weight (LW-GG) at 91 D of age. Partial least squares discriminant analysis revealed that the microbiota of the HW-TT and LW-GG females were clearly separated into 2 clusters. Furthermore, we identified 13 significantly different (P < 0.05) microbes at the genus level and 17 significantly different (P < 0.05) species between the HW-TT and LW-GG groups. Our data show that rs16775833 can modulate the microbial community structure and is associated with the BW of birds. To our knowledge, this is the first time that DMRT1 has been identified as a specific host factor, which is not only involved in sex determination but also has an effect on microbial function that might regulate animal growth.

Introduction

Gut microbes have a close symbiotic relationship with their host. Interactions between the microbial population and host genetics affect the immunologic, nutritional, and physiological status of the host (Bonder et al., 2016). Chickens are representative bird species, which have been a useful model for biological research on growth, development, immunology, and evolution (Vainio and Imhof, 1995). It is estimated that 10–100 trillion microbial cells reside in the gut of broilers, and this microbiota accounts for the largest and most diverse population (McKenna et al., 2008).

Genome-wide association studies (GWAS) have identified multiple loci that contribute to the gut microbiota population. The GWAS approach has been used to identify multiple genetic loci, including the VDR gene, which are significantly associated with the overall microbial variation and variations at the individual taxa level (Wang et al., 2016). Approximately 10% of the overall variation in the gut microbiome can be explained by 4 genetic loci, which are located in the genes SLC9A8, CBEP4, TNFSF4, and SP140 (Ruhlemann et al., 2018). Microbial QTL plots showed the association of a functional LCT SNP (rs4988235) with the abundance of Bifidobacterium (Bonder et al., 2016). In a previous study, we showed that the locus rs16775833 within the DMRT gene cluster accounted for approximately 21% of the variation in the population structure of the gut microbiota by GWAS (Ji et al., 2019).

In this study, we further provide evidence for the effect of rs16775833 variation on the gut microbiota and its association with the BW of birds.

Materials and methods

Sample Collection

A total of 777 Chinese Yellow broiler breeders were obtained from a local hatchery (Lingnan, Guangdong Wiz Agricultural Science & Technology Co. Ltd., Guangzhou, China). Before 40 D of age, each hatch was maintained in a group cage (70 cm wide, 70 cm deep, and 35 cm high). To eliminate shared microbiome influence, from 40 to 91 D, each chicken was reared in an individual cage with an individual feed trough (28 cm wide, 40 cm deep, and 41 cm). The house was equipped with nipple drinkers and was continuously lit. All chickens were allowed ad libitum feeding a starter diet (200 g/kg CP and 2,900 kcal ME/kg) from hatch to 35 D of age and followed by feeding a grower diet (180 g/kg CP and 2,950 kcal ME/kg). At 91 D of age, 22 broilers were identified based on their DMRT1 genotype and BW from population, including 11 low-weight (LW) females with genotype GG (LW-GG) and 11 high-weight (HW) females with genotype TT (HW-TT). Venous blood samples were collected and immediately stored at −80°C. Cecal contents were transferred into separate sterile tubes, then kept on ice in a biosafety cabinet, and immediately stored at −80°C. All chickens used were cared for and used in the present study in accordance with the regulations of the Institute of Animal Science, Guangdong Academy of Agricultural Sciences (No. GAAS-IAS-2019-31).

Genotyping

Genomic DNA was extracted from vein blood samples through saturated phenol–chloroform extraction method and then was stored immediately at −80°C before further analysis. The DMRT1 sequence was retrieved from the GenBank database to design PCR primers based on DNA Select program of the DNASTAR package (version 6.0). The primers were designed to amplify rs16775833 for T/G phenotype. The sequence of forward primer was 5′-CTTTTCATTTCTTAAAAGGAAACTG-3′ and reverse primer was 5′-GCAACATCTTGCCTTCACCA-3′. The PCR performed in a total volume of 25 μL, and the mixture included 2 μL primers, 1-X reaction buffer, 0.3 μL DNA polymerase, and 2 μL DNA 2.5 mM dNTP mixture. The PCR cycling conditions were initial denaturation at 95°C for 5 min, followed by 5 cycles consisting of denaturation at 94°C for 30 s, annealing at 53°C for 30 s, and extension at 72°C for 30 s and 30 cycles of denaturation at 94°C for 30 s, annealing at 60°C for 30 s, and extension at 72°C for 45 s; this was followed by a final extension at 72°C for 10 min. The PCR product (6 μL) was digested in 10 μL reaction buffer with 0.6 μL Taal at 65°C for 3 h and analyzed using 1.5% agarose gel electrophoresis.

DNA Extraction

Genomic DNA was extracted from the feces using the QIAamp DNA stool mini kit (Qiagen, Venlo, The Netherlands), following the manufacturer's protocol, and was stored at −80°C before further analysis. The quality and concentration of the extracted DNA samples were measured using agarose gel electrophoresis and the NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), respectively.

Pyrosequencing of 16S rDNA Amplicon Pyrosequencing

PCR amplification of the nearly full-length microbe 16S rRNA genes was performed using the primer with forward 27F (5′-AGAGTTTGATCMTGGCTCAG-3′) and reverse 1492R (5′-ACCTTGTTACGACTT-3′). The extracted DNA for two-step PCR, with sample-specific 16-bp barcodes in the second step of PCR, was incorporated into the forward and reverse primers for multiplex sequencing. In both the steps of PCR, the reaction mixture included Q5 High-Fidelity GC buffer (5X) for 5 μL, Q5 reaction buffer (5X) for 5 μL, Q5 High-Fidelity DNA Polymerase (5U/μL) for 0.25 μL, DNA template for 2 μL (2.5 mM) of dNTPs for 2 μL, each of the forward and reverse primers for1 μL (10 μM), and ddH2O for 8.75 μL. Thermal cycling conditions contained initial denaturation at 98°C for 2 min, followed by 25/10 cycles (for first and second PCR amplification steps, respectively) consisting of denaturation of 30 s at 98°C, annealing of 30 s at 55°C, and extension of 90 s at 72°C, with a final extension of 5 min at 72°C. PCR amplicons were purified with Agencourt AMPure Beads (Beckman Coulter, Indianapolis, IN, USA) and quantified using a PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA). Amplicons were pooled in equal amounts after the individual quantification step, and single-molecule real-time sequencing was conducted by the PacBio Sequel platform at Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China). PacBio circular consensus sequencing (CCS) reads were derived from the multiple alignments of subreads to decrease the error rate. In CCS, a ligated circular DNA template is read by the DNA polymerase multiple times, which can effectively create a consensus sequence from multiple reads of a single molecule.

Raw sequences were filtered for a minimum of 3 passes and predicted accuracy for a minimum of 99% (minfullpass = 3, minPredictedAccuracy = 99) through the PacBio SMRT Link portal (version 5.0.1.9585). The predicted accuracy is defined as the threshold lower than which a CCS is acted as noise, was 99%. The files created by the PacBio platform were then trimming to remove sequences with a length longer than 2,000 bp using amplicon size.

Sequence Analysis

The Quantitative Insights Into Microbial Ecology (QIIME, v1.8.0) pipeline was used to process the pyrosequencing data, as previously described (Caporaso et al., 2010). Raw sequence reads are available on the National Center for Biotechnology Information website and have an SRA accession number, SRP200433. Original sequencing reads that exactly matched the barcodes were assigned to respective samples and then identified as valid sequences. The remaining high-quality sequences were assigned to operational taxonomic units (OTU) at 97% sequence identity through UCLUST after chimera detection (Edgar 2010). A representative sequence was picked for each OTU using default parameters. The OTU taxonomic classification was conducted by BLAST searching the representative sequences set against the National Center for Biotechnology Information 16S ribosomal RNA Database using the best hit (Altschul et al., 1997). An OTU table was generated to record OTU information (abundance and taxonomy of OTU) in each sample. The OTU containing <0.001% of the total sequences across all samples were removed. To minimize the difference in sequencing depth across samples, an average, rounded rarefied OTU table was created by averaging 100 evenly resampled OTU subsets lesser than 90% of minimum sequencing depth for further analysis.

Data Analysis

Analyses of sequence data were largely performed using Quantitative Insights into Microbial Ecology and R packages (v3.2.0). Differences in the UniFrac distances for pairwise comparisons among groups were determined by Student t test and the Monte Carlo permutation test with 1,000 permutations and visualized through constructing box-and-whisker plots. The taxonomy abundances were statistically compared among samples or groups by Metastats at the phylum, class, order, family, genus, and species levels (White et al., 2009) and presented as violin plots. Partial least squares discriminant analysis (PLS-DA) was also recommended as a supervised model to reveal the microbiota variation among groups, through the “plsda” function in R package “mixOmics” (Chen and Jiang, 2014).

Association analyses of polymorphisms with BW in group with 777 broilers were evaluated by the GLM procedures of SAS 8.0 software (SAS Institute Inc., Cary, NC, USA) through the following model:where Yijk was the weights at 4, 7, and 13 wk of age, μ was the overall population average value, Hi was the fixed effect of hatch i, Fj was the effect of sire j, Gk was the fixed effect of genotype k, and eijk was the residual effect. Comparisons between 2 groups were determined using the Student t test, and the significance of differences was considered at P < 0.05.

Results

Identification of rs16775833 Genotypes in Birds

The rs16775833 SNP is located in DMRT1, which is located at 26185353 on Z chromosome in chicken (Version 5.0). The 2 genotypes, GG and TT, were observed. An amplicon of 227 bp was amplified using the primer pair (Figure 1A). On digestion with Taal, the amplicon from the GG genotype gave a single band of 227 bp, whereas in the case of TT genotype, 2 bands of 162 and 65 bp were obtained (Figure 1B).

Figure 1

Agarose gel picture showing the 2 genotypes of rs16775833. (A) PCR-amplified product of rs16775833. (B) The amplified product of rs16775833 digested with Taal. The SNP genotypes of rs16775833: left lanes (1)–(11) = genotype GG; right lanes (1)–(11) = genotype TT.

rs16775833 is Associated With the BW in Birds

To investigate the role of rs16775833 polymorphism on phenotypic traits of chickens, the GG and TT genotypes were correlated with the BW at different times (at 4, 7, and 13 wk of age). A total of 777 females were genotyped for rs16775833 in DMRT1, and 481 females were found to have the GG genotype, 128 female chickens had the TT genotype, and 168 female chickens had TG genotype.

The association analysis of the rs16775833 genotype with chicken BW was performed. At the age of 4 wk, there was no significant difference between the phenotypes of GG and TT with respect to the BW (P < 0.05). The BW of birds with the TT genotype was higher than those with the GG genotype at 7 and 13 wk of age (Figure 2). The association analysis indicated that rs16775833 was associated with BW in the chicken.

Figure 2

Association of GG and TT genotypes of rs4958843 with the BW of chickens. The BWof chickens with GG and TT genotypes at 4, 7 and 13 wk of age. Data are shown as means ± SEM. Comparisons between means were performed using t tests (∗P < 0.05). (n = 481, GG genotype group; n = 128, TT genotype group).

The Diversity of the Microbial Community Between the HW-TT and LW-GG Groups

Fecal samples from 22 birds which included 11 HW-TT females and 11 LW-GG females were collected. The average number of high-quality sequences from bacterial populations generated per sample was 11,026. The average length of sequence reads was 1,462 bp. The rarefaction curves were close to saturation. Based on 97% sequence similarity, all the sequences were clustered into 3,820 bacterial OTU (Supplementary Figure 1). The number of shared OTU between the HW-TT and LW-GG groups was 861, which is shown as a Venn diagram in Figure 3A. Partial least squares discriminant analysis showed a clear separation of samples between the HW-TT and LW-GG groups in a model based on variable importance in projection (Figure 3B). For α-diversity, the index of Chao, Shannon, and Simpson showed no significant difference between the 2 groups (Supplementary Figure 2). At the phylum level, the HW-TT and LW-GG groups were similar (P > 0.05) for fecal microbial community structure. At the phylum level, the percentage distribution of the microbiome community in HW-TT and LW-GG chickens was 64.2 and 61.1 for Bacteroidetes, 21.8 and 24.8 for Firmicutes, 3.4 and 5.7 for Verrucomicrobia, 3.3 and 2.8 for Proteobacteria, 1.7 and 1.1 for Synergistetes, 1.1 and 1.4 for Tenericutes, respectively, and less than 1 for other phyla (Figure 3C). These results demonstrated that chicken GG and TT genotypes can be a clear separation, but there was no significant difference at the phylum level.

Figure 3

(A) Venn diagram showing the distribution of shared operational taxonomic units between the HW-TT and LW-GG groups. (B) The partial least squares discriminant analysis of operational taxonomic units. (C) Relative abundance of microbes between the HW-TT and LW-GG groups at the phylum level. Abbreviations: HW-TT, TT genotype chickens with high weight; LW-GG, GG genotype chickens with low weight; PLS-DA, partial least squares discriminant analysis.

Bacterial Composition at the Genus Level

To further characterize the changes in the microbial population structure that were imposed by the rs16775833 genotype TT and GG, 16S rRNA data were classified taxonomically at the genus level. For statistical analysis to detect the effects of the host quantitative genotype, 187 genera were used. Of these, 13 differed significantly (P < 0.05) between the HW-TT and LW-GG birds (Figures 4A, 4B). Further analysis showed that the most significant difference genus (P < 0.05) between the HW and LW birds was Anaerobacterium (LW-GG > HW-TT). These results demonstrate that limitations of chicken GG and TT genotypes influence the population structure of the gut microbiome at the phylum level, but some differences are present in the microbiome population between HW-TT and LW-GG birds at the genus level.

Figure 4

(A) Relative abundance of microbes between the HW-TT and LW-GG groups at the genus level. (B) Abundance of the significantly different genera. Abbreviations: HW-TT, TT genotype chickens with high weight; LW-GG, GG genotype chickens with low weight.

Diversity of the Microbial Community at the Species Level

At the phylum level, 428 genera were used. Of these, 17 differed significantly (P < 0.05) between the HW-TT and LW-GG birds (Figures 5A, 5B). Further analysis showed that the Bacillus toyonensis was the most significantly different (P < 0.05) between the HW and LW birds (HW-TT > LW-GG). These results demonstrate that the chicken GG and TT genotypes of rs16775833 also affect the microbiome population between HW-TT and LW-GG of birds at the species level.

Figure 5

(A) Relative abundance of species between the HW-TT and LW-GG groups. (B) Abundance of the significantly different species. Abbreviations: HW-TT, TT genotype chickens with high weight; LW-GG, GG genotype chickens with low weight.

Discussion

Our previous study showed that the loci rs16775833 was associated with β-diversity, being located within the DMRT1 gene cluster, which might influence the BW of chickens (Ji et al., 2019). The DMRT1 gene family is widely known for its involvement in sex determination in many organisms (Lambeth et al., 2014, Ayers et al., 2015). To further assess the effect of DMRT1 phenotype on the BW and gut microbiota, TT and GG genotypes of rs16775833 were identified in birds.

We observed that 128 birds with the TT genotype had a higher BW than 481 ones with the GG genotype at 7 and 13 wk of age (Figure 2). The gut microbiome is a complex ecosystem that is functionally involved in various biological processes. Changes in the gut microbial composition are influenced by both the host genotype and environmental factors (Neish, 2009, David et al., 2014). Therefore, we further investigated the difference on the intestinal microbial communities between the HW-TT and LW-GG birds based on the 16S rRNA profiling.

To compare the microbial diversity across groups, PLS-DA and Venn diagrams were used. The Venn diagram, used to make qualitative comparisons between the HW-TT and LW-GG chickens of the same (71 D) age, indicated that 861 OTU were common to the HW-TT and LW-GG groups. Partial least squares discriminant analysis analysis further showed a clear separation of samples in HW-TT and LW-GG groups. This suggests that rs16775833 might be an important factor in determining the microbial composition of the gut in chickens. We further analyzed the differences at the phylum level but found was no significant differences between the 2 genotype groups.

Analysis at the species level showed significant differences in 17 species, among which Bacillus toyonensis (HW-TT > LW-GG) and Anaerobacterium chartisolvens (LW-GG > HW-TT) were the most significantly different; B. toyonensis is used as a feed additive and has been shown to have beneficial effects when administered to animals. It can act against enteric pathogens and stimulates Na-dependent glucose absorption in pigs (Breves et al., 2000, Kantas et al., 2015). In addition, B. toyonensis strains can be used to enhance the gain in BW, improve the feed conversion ratio, and significantly improves the livability of birds (Grela et al., 2009, Batkowska et al., 2015). The chickens harboring B. toyonensis were reported to have more abdominal fat, and their meat had significantly higher conductivity than the meat of chicken in the control group (Novak et al., 2011). The Anaerobacterium OTU were associated with low feed efficiency. Anaerobacterium is negatively correlated with total feed intake, total BW gain (Siegerstetter et al., 2018).

In summary, there are measurable differences in compositional microbiota between chicken GG and TT genotypes. The results of present study indicated that the DMRT1 genotype strongly influences the BW in birds and is associated with the change in intestinal microbiota. The results provide new insights into the biological function of the DMRT1 gene.