Metagenomic Analysis Identifies Sex-Related Cecal Microbial Gene Functions and Bacterial Taxa in the Quail

Background: Japanese quail (Coturnix japonica) are important and widely distributed poultry in China. Researchers continue to pursue genetic selection for heavier quail. The intestinal microbiota plays a substantial role in growth promotion; however, the mechanisms involved in growth promotion remain unclear. Results: We generated 107.3 Gb of cecal microbiome data from ten Japanese quail, providing a series of quail gut microbial gene catalogs (1.25 million genes). We identified a total of 606 main microbial species from 1,033,311 annotated genes distributed among the ten quail. Seventeen microbial species from the genera Anaerobiospirillum, Alistipes, Barnesiella, and Butyricimonas differed significantly in their abundances between the female and male gut microbiotas. Most of the functional gut microbial genes were involved in metabolism, primarily in carbohydrate transport and metabolism, as well as some active carbohydrate-degrading enzymes. We also identified 308 antibiotic-resistance genes (ARGs) from the phyla Bacteroidetes, Firmicutes and Euryarchaeota. Studies of the differential gene functions between sexes indicated that abundances of the gut microbes that produce carbohydrate-active enzymes varied between female and male quail. Bacteroidetes was the predominant ARG-containing phylum in female quail; Euryarchaeota was the predominant ARG-containing phylum in male quail. Conclusion: This article provides the first description of the gene catalog of the cecal bacteria in Japanese quail as well as insights into the bacterial taxa and predictive metagenomic functions between male and female quail to provide a better understanding of the microbial genes in the quail ceca.


INTRODUCTION
Japanese quail (Coturnix japonica), from the order Galliformes and family Phasianidae, are an important poultry species for egg and meat production and are widely distributed in China. Quail have a short maturation period, and female quail can lay their first eggs at 35 days old. The laying rate can reach 50% at 45 days, and the weight of annual egg production is an average of 20-25 times that of the female quail's body weight. Quail make ideal animal models owing to their small body size and short generation intervals. Quail byproducts also have commercial applications, for example, feathers in duvet manufacturing. Therefore, some researchers consider quail breeding to be the future of twenty-first century poultry breeding (1,2).
The gut microbiota plays an important role in production and disease resistance in many animals, especially in digestion and nutrient absorption, contributing to feed-related traits (3)(4)(5)(6)(7)(8)(9)(10). 16S rRNA genes were analyzed to profile the cecal bacterial communities of chicken; these analyses showed that the male animals had high abundances of Bacteroides, whereas female animals were enriched with Clostridium and Shigella. Microbiome-wide association analyses in the ilea and ceca of Japanese quail have shown that several quantitative feed-related traits, including feed or nutrient efficiency, feed intake, body weight gain, feed per gain ratio, and phosphorus utilization were associated with microbiota features at both the bacterial genus and operational taxonomic unit levels as characterized by 16S rRNA amplicon sequencing (11)(12)(13)(14)(15)(16). Beneficial groups of bacteria (i.e., probiotics) in the gut help balance the intestinal flora and improve the body's disease resistance via competitive inhibition of acid substance production. Bacillus and some subspecies of Enterococcus are considered probiotics in chickens and Japanese quail (17,18). However, to our knowledge, the reference gene catalogs of the gut microbiome, including the intestinal bacterial compositions and the functional capacity of the gut microbiome, have rarely been applied or reported for Japanese quail. This could hinder growth promotion in quail.
Female quail exhibit greater growth potential in the breast, wings and back than do male quail for both yellow and red laying Japanese quail, but no relevant research is available on the gut microbiotas of adult quail in China (1,2). Variations in growth rates between the sexes might be due to differences in their gut microbiomes because the gut microbiome has key effects on nutrient digestion, absorption, and metabolism in quail (19,20). However, knowledge of how the gut microbiome varies between sexes in quail is minimal. This study was conducted to examine how nutrient digestion and absorption by the gut microbiome differ between male and female quail and whether any gut bacteria are beneficial to quail.
Here, we used shotgun metagenomic sequencing to analyze the microbiomes from cecal samples of ten 70-day-old Japanese quail. At this age, the cecal microbiota should be mature and stable (21). We compared the gut microbiotas between female and male quail to study the gut microbial ecosystem and the differences in the gut microbiotas between the sexes in this economically important species.

Sample Collection
Samples were collected from the cecal contents of ten 70-dayold Japanese quail (one sample per quail), immediately frozen in liquid nitrogen, and stored at −80 • C until DNA extraction. We used five healthy adult female quail and five healthy adult male quail. All birds were offered the same diet and housed under similar environmental conditions. Supplementary Table 1 provides the details of the samples.
DNA was extracted using the Magen DNA Stool Kit per the manufacturer's protocol (Magen, Guangdong, China), using 200 mg of feces per sample. The DNA samples were tested using two methods. First, the degree of DNA degradation was determined on 1% agarose gels, and second, the DNA concentration was measured using a Qubit R dsDNA Assay Kit in a Qubit R 2.0 Fluorometer (Life Technologies, CA, USA). The optical density was between 1.8 and 2.0. The DNA contents were used to construct a sequencing library.

Library Construction
One microgram of DNA per sample was used as the input material to prepare the DNA samples. Sequencing libraries were generated using NEBNext R Ultra TM DNA Library Prep Kit for Illumina (NEB, USA), and recommendations and index codes were added to attribute sequences to each sample. Briefly, the DNA sample was fragmented by sonication to 350 bp, then the DNA fragments were end-paired, poly-A-tailed, and ligated with the full-length adaptor for Illumina sequencing with further PCR amplification. PCR products were then purified (AMPure XP system), and libraries were analyzed for size distribution using an Agilent 2100 Bioanalyzer and quantified using real-time PCR.

DNA Sequencing and DNA Assembly
Index-coded samples were clustered using a cBot Cluster Generation System per the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq platform, and paired-end reads were generated.
After preprocessing the raw data, clean data were obtained using subsequent analyses, which included removing reads that contained low-quality bases (default quality threshold value ≤38) above a certain portion (default length: 40 bp), removing reads in which the N base had reached a certain percentage (default length: 10 bp), and removing reads that overlapped by more than a certain portion using Adapter (default length: 15 bp). To remove the effects of host pollution, clean data were analyzed using Basic Local Alignment Search Tool (BLAST) against the Japanese quail genome database (INSDC: LSZS01000000), which defaults to using Bowtie2.2.4 software (Bowtie2.2.4, http:// bowtiebio.sourceforge.net/bowtie2/index.shtml) to filter reads of host origin using the following parameters (22,23): -end-to-end, -sensitive, -I 200 and -X 400.

Construction of the Gene Catalog and Abundance Analysis
To construct a comprehensive Japanese quail gut microbial gene catalog, all reads were assembled de novo from each sample into longer contigs, which were assembled and analyzed (24) using SOAPdenovo software (V2.04, http://soap.genomics.org. cn/soapdenovo.html), with the parameters -d 1, -M 3, -R, -u, -F and -K 55 (25)(26)(27)(28). The assembled scaffolds were then broken from the N junction to acquire scaffolds not containing N (called scaftigs) (26,29,30). The clean data were compared with each respective scaffold using Bowtie2.2.4 to acquire unused pairedend reads with the parameters -end-to-end, -sensitive, -I 200 and -X 400 (26). Fragments shorter than 500 bp were filtered from all scaftigs because of the lengths of the assembly sequences.
The assembled scaftigs (≥500 bp) were used to predict the open reading frames (ORFs) using MetaGeneMark (V2.10, http://topaz.gatech.edu/GeneMark/) software, then the sequences with coding region lengths <100 nt were filtered from the predicted results using the default parameters (26,(30)(31)(32)(33). CD-HIT software (34, 35) (V4.5.8, http://www.bioinformatics. org/cd-hit) was used to obtain the unique initial gene catalog (the genes referred to here are the nucleotide sequences coded by unique and continuous genes) (36), with the parameters -c 0.95, -G 0, -aS 0.9, -g 1, and -d 0 (33,36). Each gene's abundance was calculated from the number of mapped reads and the gene length. The formula used was r represents the number of reads mapped to the genes, and L represents gene length (31). The gene numbers were analyzed from the abundance of each gene in each sample in the gene catalog; the analyses included the basic statistical information, core-pan gene analysis, correlation analysis of the samples and a Venn diagram.

Taxonomic Assignment of Genes and Construction of Taxonomy and Relative Abundance Profiles
DIAMOND [28] software (V0.9.9) was used to analyze the unigenes via BLAST against the bacterial, fungal, archaeal and viral sequences, which were extracted from the NR database (Version: 2018-01-02, https://www.ncbi.nlm.nih.gov/) of the National Center for Biotechnology Information with the parameter settings blastp and -e 1e −5 . Because each sequence could have multiple aligned results, sequences were chosen for alignment if the E-value was ≤ the smallest E-value ×10 (45). The least common ancestors (LCA) algorithm was also applied to the system classification using MEGAN software (46) to confirm the species annotation.
The number of genes and the abundance information for each sample were obtained using the LCA annotation results and the gene abundance for each taxonomic hierarchy (kingdom, phylum, class, order, family, genus, and species). The species abundance in one sample equaled the sum of the gene abundances annotated for the species. The gene number of a species in a sample equaled the number of genes with abundances >0.
These analyses were based on the abundance results for each taxonomic hierarchy, including the relative abundance and abundance cluster heat maps. Differences between two groups were tested via analysis of similarity (ANOSIM; R vegan package, version 2.15.3). Metastats and linear discriminant analysis effect size (LEfSe) analysis were used to detect different species between groups. Permutation tests between groups were conducted in Metastats for each taxonomic level to obtain the p-value, then the Benjamini and Hochberg false discovery rate procedure was used to correct the p-value and acquire a q-value (47). LEfSe analysis was conducted using LEfSe software, with a default LDA score of 3. Important species were screened out via MeanDecreaseAccuracy and MeanDecreaseGini, and the receiver operating characteristic curve was plotted to cross validate each model (default: 10 times).

Microbial Genome Sequencing and de novo Metagenome Assembly
We sequenced 107.3 Gb of high-quality data, with an average of 10.73 Gb per sample (SRA accession: SUB9360165) after removing low-quality reads, N reads and adaptor sequences. This included 99.5 Gb of no-host data. Table 1 lists the raw, clean and preprocessed data. The best assembly with different k-mer sizes was chosen on the basis of contig N50 and the mapping rate. The lengths of the contig N50 in each sample ranged from 1,662 to 4,217 bp. Table 2 summarizes the assembly results. The length distribution of the contigs ranged from 500 to 3,000 bp, with most of the length distribution falling between 500 and 1,000 bp (Supplementary Figure 1).

Gut Microbiome Gene Catalog Prediction and Taxonomic Annotation
Based on the scaftigs, 3,378,594 ORFs were constructed, and core-pan gene analysis showed that the curve approached saturation (Supplementary Figure 2). We obtained 1,247,092 ORFs after reducing the data redundancy; 48.76% of these were complete ORFs, and 608,045 genes were obtained. The female quail contained 303,162 common ORFs; the male quail contained 401,051 common ORFs (Supplementary Figure 3).
The gene abundances revealed several dominant microbial species in both groups at different levels. The top five were assigned to Bacteroidetes, Firmicutes, Proteobacteria, Spirochaetes, and Fusobacteria ( Figure 1A).

Gene Function Analysis
To obtain gene function information, the genes were analyzed separately via BLAST against the KEGG, EggNOG, CAZy, and CARD databases using DIAMOND software. For the functional annotation, 804,714 genes were annotated in the KEGG database, followed by 791,180 in EggNOG, 46,317 in CAZy, and 496 in CARD (Supplemental File 5). Of the 496 genes in CARD, 308 were for antibiotic resistance.
The KEGG database profile on the relative abundances of the gene functions suggested that the main activities of the genes identified were associated with metabolism, genetic information processing and environmental information at the first classification level. Genes related to the three terms of replication, recombination and repair, amino acid transport and

EggNOG
The BLAST results based on EggNOG showed that more unigenes were in the classes of amino acid transport and metabolism (55,478 genes) and carbohydrate transport and metabolism (54,316 genes), while a cluster of 43,338 genes was assigned to inorganic ion transport and metabolism.

ANOSIM and LEfSe Analysis
To examine the similarities between the male and female quail, we analyzed the differences in relative gene abundances and their corresponding functions in the ten quail. CAZy results showed differences between the groups ( Figure 4A); however, the differences between the measurement data were not statistically significant. The results suggested that the abundances of gut microbes that produce carbohydrate-active enzymes may differ between female and male quail.
To investigate the difference in functional capacities of gut microbiome between female and male quail, we performed a LEfSe analysis, the significant difference was yielded based on the result of Cazy analysis. The genes of GH31 annotation was significantly different between female and male quail ( Figure 4B).

Analysis of Antibiotic-Resistance Genes and Bacteria
In contrast to the above results, fewer genes were related in the CARD database, which matched 496 genes; 308 of these were antibiotic-resistance genes (ARGs), with an abundance of 2.98 × 10 −4 . Of these genes, 110 were found in the ten quail (Supplemental File 7). The most abundant genes included tetQ, tetW, ermF, adeF, tetW/N/W, lnuC, sul2, tet40, tetO, and APH3-IIIa. Species information was collected for each ARG from each quail. These ARGs were harbored by diverse phyla, including Bacteroidetes, Firmicutes, and Euryarchaeota ( Figure 5).
Analysis of the ARG distributions showed that the female quail had 147 common genes and 24 unique genes, and the male quail had 124 common genes and ten unique genes (Supplemental Figure 4 and Supplemental File 7). Bacteroidetes were the predominant ARG hosts in the female quail; Euryarchaeota were the dominant ARG hosts in the male quail. Tenericutes and Kiritimatiellaeota were detected only in the female quail ( Figure 5).

DISCUSSION
We conducted an Illumina short-read-based gene-annotation sequencing analysis of metagenomic DNA taken from the ceca of ten Japanese quail, including five female quail and five male quail, to obtain a catalog of bacterial gene communities in adult quail. To our knowledge, this is the first analysis of the gut microbiota gene compositions for this species. The results from the gene annotation against the MicroNR database revealed several abundant groups at the phylum level, including Bacteroidetes, Firmicutes, Proteobacteria, Spirochaetes and Fusobacteria. 16S rRNA gene sequencing revealed similar results for the ileal microbiotas of this species (15). Additionally, Tenericutes was a major phylum in cecal samples from Japanese quail (12); however, the abundance of Firmicutes was higher than that of Bacteroidetes (12,15), both of which likely dominated the cecal microbiota (11,12,15). Other studies have shown that higher abundances of Firmicutes could lead to more efficient absorption of calories from the digesta (48,49), which may be related to a higher abundance of Firmicutes in the ileum and a lower abundance in the cecum. Bacteroides spp. play key roles in immunoregulatory FIGURE 4 | (A) Analysis of similarity (ANOSIM) between female (F) and male (M) Japanese quail (Coturnix japonica) based on the CAZy database. The horizontal axis represents the grouping information; the vertical axis indicates the distance data, and "Between" is the combined information for both groups. The median line for "Between" being higher than those of the other two groups suggests that the groupings were appropriate. The R-value was between −1 and 1 and >0, indicating a significant difference between the groups. The P-value represents the reliability of the statistical analysis, and P < 0.05 indicates statistical significance; (B) Different functional capacities of gut microbiome based on Carbohydrate-Active Enzymes (Cazy) analysis between female (F) and male (M) Japanese quail (Coturnix japonica).
abilities and help provide valuable nutrients and energy to the host by breaking down food (50,51). Bacteroidetes may also contribute to disease resistance in the host. Moreover, Bacteroides can degrade indigestible fiber in chicken guts (20), and the enriched Bacteroides from the quail ceca was consistent with substantial microbial fermentation in the hindgut.
Our gene function prediction results showed that the microbial species in the quail ceca contributed to metabolic functions involving production of energy and nutrients by digesting food, showing that the cecal microbiota is an important energy supplier for the host. The GH31 might contribute to metabolic functions, since the functional capacity analysis of gut microbiome found the GH31 annotation was significantly different between female and male group. These results reflected the important role of the cecal microbiota in some feedrelated traits, such as phosphorus utilization, daily gain, feed intake and feed per gain ratio reported in Japanese quail as assessed by mixed linear models (15). Our data revealed genes belonging to subspecies of Bacillus and Enterococcus, which are considered probiotics in chickens and Japanese quail (17,18). Bacillus exerts positive effects on several feedrelated traits.
Our analysis identified ∼3.4 million comprehensive Japanese quail gut microbial genes. Animal gut microbiotas, including those of humans (31,32,52), dogs (53), monkeys (54), mice (55), rats (56,57), chickens (58), and pigs (25,59), are receiving increasing attention worldwide for their important contributions to host nutrition and health. The size of the Japanese quail gene catalog was smaller than those of chickens (9.04 million genes), humans (9.9 million genes) and pigs (7.7 million genes) (58), possibly owing to the limited samples from only the lumens of the ceca, which do not contain gene information on the microbiome in the foregut. From our data, means of 88.74, 0.59, 0.0364, and 0.203% of the genes were annotated to bacteria, viruses, eukaryotes and archaea, respectively. No viral genes have been found in chicken gut microbiotas (58). Metagenomic sequencing identified 308 ARGs, of which, the dominant genes were tetQ, tetW, ermF, adeF, tetW/N/W, lnuC, sul2, tet40, tetO, and APH3-IIIa. Some ARGs, such as tetM, vanX and bla, existed before the advent of antibiotic use (60). The ARG abundances in the quail (P = 2.98 × 10 −4 ) in our study were lower than those reported in the pigs (61) and were mainly distributed in Bacteroidetes, Firmicutes and Euryarchaeota. High-throughput next-generation sequencing enables detecting the presence of broad-spectrum antibiotic resistance in animal gut fecal resistomes. To date, numerous studies on ARGs in regard to existing and emerging antibioticresistance threats have been reported for humans, sheep, chickens, pigs, and cows (62)(63)(64). Because of the intensive link between the spread of ARGs and human health, increasing knowledge of ARGs is important for addressing potential threats to human health.
Our results showed that female quail had relatively high microbial species abundances, whereas male quail had a higher microbiome diversity (Chao1 index) (12). The differences between the male and female quail were analyzed, and 17 species that were clustered into four genera (Anaerobiospirillum, Alistipes, Barnesiella, and Butyricimonas) from four families (Succinivibrionaceae, Rikenellaceae, Barnesiellaceae, and Odoribacteraceae) were differentially abundant. Compared with the families Lactobacillaceae and Catabacteriaceae, which were differentially abundant in the large intestines of male and female Japanese quail (12), our results revealed significant differences within the quail. The differences in these microbiomes may be related to the differences in growth potential between male and female quail. Bacterial communities also differ significantly between male and female broiler chicken ceca. Male chicken ceca were enriched with Bacteroides, whereas female chicken ceca were enriched with Clostridium and Shigella (20). These detected bacterial species varied between chickens (65) and quail. However, some of our results are preliminary, and the functions of these gut microbes require further study.
In summary, our study provides the first reference gene catalog for the Japanese quail gut microbiome, which will be an important addition to animal gut metagenomics. Metagenomic analysis will contribute to future studies on the differences in the mechanisms of feed-related quantitative traits and the gut microbiome in quail. Our results also help explain why female quail exhibit greater muscle growth potential than do male quail.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: SRA accession: SUB9360165.

ETHICS STATEMENT
The animal study was reviewed and approved by Animal Care Committee of Nanchang Normal University (Nanchang, China).

AUTHOR CONTRIBUTIONS
J-EM, X-WX, and J-GX carried out the design and wrote the manuscript. J-SG, Y-FL, QX, and Y-BY helped to revise the draft manuscript. JL, X-NZ, Y-WT, and W-TS help to collect the samples. X-TT, Z-FW, and MZ help to analyze the data. C-YZ sorted some of the tables. X-QZ provided some effective suggestions. Y-SR designed the study and revised the manuscript.