Metagenomic sequencing and gene prediction
Metagenomic sequencing generated 291 billion raw reads corresponding to 123 Gbp of raw data from the cecum content of the 24 broilers. Perform low-quality read filtering on each sample data to obtain clean reads (Supplementary Table S1). After assembling and predict sequencing data, display an information about reads, length of open reading frame (ORF), N50 etc. of each broiler is shown in Table 1, and the gene length distribution of the non-redundant gene set is shown in histogram (Fig. 1a). We identified 7.11 million non-redundant genes, with an average ORF length for 551 bp. Moreover, we investigated the difference in gene number among six group, and found that the number of non-redundant genes in higher in FKXHC compared to other five groups, but the difference is not significant as shown in the scatter plot (Fig. 1b). The numbers of common and unique genes among the six groups are shown in (Fig. 1c), and 198778 genes were found to be common in both regions.
Table 1. Sequenced reads analysis and assembly statistics
Sample
|
Contig Num.
|
Total Len. (bp)
|
N50 (bp)
|
GC (%)
|
Mapped (%)
|
FKXHC1
|
392222
|
555133757
|
1777
|
48.05
|
94.65
|
FKXHC2
|
433139
|
628311260
|
1884
|
49.23
|
93.86
|
FKXHC3
|
262011
|
458872158
|
2874
|
51.35
|
95.88
|
FKXHC4
|
407587
|
608707582
|
2009
|
48.44
|
94.35
|
HZHXC1
|
237593
|
361955020
|
2072
|
48.56
|
88.32
|
HZHXC2
|
296745
|
484357564
|
2374
|
48.34
|
95.15
|
HZHXC3
|
308913
|
464894041
|
2005
|
48.71
|
93.88
|
HZHXC4
|
295257
|
476343903
|
2362
|
51.55
|
94.65
|
QYMC1
|
219833
|
401060089
|
3327
|
48.47
|
94.28
|
QYMC2
|
239798
|
355326828
|
1943
|
48.01
|
88.95
|
QYMC3
|
321474
|
525704855
|
2419
|
48.74
|
94.64
|
QYMC4
|
263335
|
426235271
|
2359
|
50.54
|
93.82
|
XYHXC1
|
394324
|
549891671
|
1734
|
47.23
|
92.43
|
XYHXC2
|
379211
|
540337865
|
1832
|
48.05
|
94.05
|
XYHXC3
|
401311
|
571101588
|
1812
|
50.01
|
92.6
|
XYHXC4
|
385463
|
546185646
|
1823
|
48.97
|
92.38
|
YSC1
|
237314
|
226163115
|
893
|
46.98
|
77.19
|
YSC2
|
375498
|
508359646
|
1646
|
49.05
|
91.58
|
YSC3
|
338497
|
498029370
|
1936
|
47.65
|
93.23
|
YSC4
|
370584
|
528209061
|
1811
|
47.8
|
91.82
|
ZSSLC1
|
282670
|
442677074
|
2215
|
48.64
|
94.43
|
ZSSLC2
|
348071
|
512786581
|
1914
|
48.06
|
95.04
|
ZSSLC3
|
315508
|
464274920
|
1903
|
48.85
|
94.97
|
ZSSLC4
|
348665
|
518056690
|
1978
|
47.58
|
95.03
|
Microbial Taxonomy Annotation
The scaftigs were aligned against the non-redundant protein (NR) database to get the microbial population composition and relative abundance information. Phylum and genus level distributions for individual samples are shown in Fig. 2a and b. At the phylum level, Bacteroidetes and Firmicutes were two dominant phyla in all samples, accounting for 75-98% of the total bacterial community. At the genus level, 3072 genera were identified across all group (Supplementary Table 2). The Bacteroides and Prevotella were dominant genus in QYMC group (13.9% and 15.5%, respectively), and the dominant genus of the other five groups is only Bacteroides. Interestingly, we compared the abundance of the top 30 abundance genus in the six groups, and found that only the QYMC group had a significantly higher level of Prevotella than the other five groups (Fig. 2d and e). To distinguish between the variations of host gut microbiomes, the α-diversity and β-diversity of microbial communities among different chicken groups were also evaluated. The chao1 index (the abundance within samples) demonstrated the number of genera in XYHXC group was higher than in other groups, whereas the difference is not obvious. The shannon index (the diversity within samples) of the cecum microflora was more diverse in HZHXC and ZSSLC group (Supplementary Figure S1). In addition, the results of Principle Coordination Analysis (PCoA) and Non-Metric Multi-Dimensional Scaling analysis (NMDS) for dimension reduction analysis based on the binary_jaccard distance revealed the dissimilarity of bacterial communities among all groups (Fig. 3a and b). The closer the sample points on the scatter chart, the greater the similarity of the composition. Analysis of similarity (ANOSIM) is a nonparametric test that assesses whether variation between groups is significantly greater than variation within groups, which helps to evaluate the reasonability of the division of groups (R=0.747, P=0.001) (Fig. 3c).
Functional annotation
We then predicted the functional and metabolic pathways of the broilers’ gut microbiome non-redundant (NR) genes sequenced based on the following functional databases: Kyoto Encyclopedia of Genes and Genomes (KEGG), Evolutionary genealogy of genes, Non-supervised Orthologous Groups (eggNOG), and Carbohydrate-Active enzymes Database (CAZy). From KEGG pathway results, we found that a large number of pathways belonged to metabolism groups (77.0%), including carbohydrate metabolism (10.3%), Amino acid metabolism (8.6%), nucleotide metabolism (6.2%), metabolism of cofactors and vitamins (6.1%), energy metabolism (4.4%), glycan biosynthesis and metabolism (2.7%), and lipid metabolism (2.4%) (Fig. 4a). The eggNOG analysis showed most of the gene functions remained unclear (18%), while the known functions are relatively more abundant in replication, recombination and repair (6.7%); translation, ribosomal structure and biogenesis (5.3%); cell wall (4.8%) and amino acid transport and metabolism (4.6%) (Fig. 4b). Meanwhile, the CAZy database analysis exposed that the most enzymes are classified as glycoside hydrolases (44.9%), followed by classified by glycosyltransferases (22.5%), polysaccharide lyases (13.5%) and carbohydrate esterases (11.8%) (Fig. 4c).
We further compared α-diversity in the individuals at the KEGG orthology (KO) level, Cluster of Orthologous Groups (COG) level and family level of CAZy levels (Supplementary Figure S2). There was no significant difference between the chao1 and shannon indies at the KO level and COG level, except that chao1 index at the COG level was significantly different between the XYHXC and QYMC group (Supplementary Figure S2). Moreover, The Shannon index shown that the diversity is similar in all groups on the KO and COG levels, and XYHXC group was significantly higher diversity than other groups on the family level of CAZy.
Antibiotic resistance gene (ARG) annotation
We used the annotation results compare with the CARD, and a total of 989 ARGs were identified in all samples as shown in Table S2. There were some differences in ARGs numbers in the six breeds of broiler. Especially, the number of ARGs in HZXHC and QYMC group were lower than in other groups ((Fig. 5a). Starting from the relative abundance table of AGRs, the ARGs results of the maximum abundance ranking top 20 were shown as histogram (Fig. 5b), and the ARGs abundance in ppm and relative percentage in each broiler was calculated. The results showed that the abundance and percentage of tet(Q) was higher than other ARGs in all samples. Identified ARGs were further categorized on the basis of their resistance profiles, and each ARG was annotated with information of resistance type. The overall presence of these resistance types described as a circos plot (Fig. 6a), which indicated the most predominant types of ARGs in the six groups of chicken are tetracyclines, multidrug and aminoglycoside, and found that the KFXHC group is the highest proportion of resistance type among the six groups. Interestingly, the highest abundance of ARG types in the XYHXC group was also found by comparing the dominant types of ARG in 6 groups of broilers (Fig. 6b), the result analysis implied that Xinghua chicken is more resistant to antibiotics than the other five yellow-feathered broilers.