Mapping of Quantitative Trait Loci Controlling Egg-Quality and -Production Traits in Japanese Quail (Coturnix japonica) Using Restriction-Site Associated DNA Sequencing

This research was conducted to identify quantitative trait loci (QTL) associated with egg-related traits by constructing a genetic linkage map based on single nucleotide polymorphism (SNP) markers using restriction-site associated DNA sequencing (RAD-seq) in Japanese quail. A total of 138 F2 females were produced by full-sib mating of F1 birds derived from an intercross between a male of the large-sized strain with three females of the normal-sized strain. Eggs were investigated at two different stages: the beginning stage of egg-laying and at 12 weeks of age (second stage). Five eggs were analyzed for egg weight, lengths of the long and short axes, egg shell strength and weight, yolk weight and diameter, albumen weight, egg equator thickness, and yolk color (L*, a*, and b* values) at each stage. Moreover, the age at first egg, the cumulative number of eggs laid, and egg production rate were recorded. RAD-seq developed 118 SNP markers and mapped them to 13 linkage groups using the Map Manager QTX b20 software. Markers were spanned on 776.1 cM with an average spacing of 7.4 cM. Nine QTL were identified on chromosomes 2, 4, 6, 10, 12, and Z using the simple interval mapping method in the R/qtl package. The QTL detected affected 10 egg traits of egg weight, lengths of the long and short axes of egg, egg shell strength, yolk diameter and weight, albumen weight, and egg shell weight at the beginning stage, yellowness of the yolk color at the second stage, and age at first egg. This is the first report to perform a quail QTL analysis of egg-related traits using RAD-seq. These results highlight the effectiveness of RAD-seq associated with targeted QTL and the application of marker-assisted selection in the poultry industry, particularly in the Japanese quail.


Introduction
Tremendous progress has been made in poultry egg production and breeding programs over the past decades, leading to the fast-growing livestock industry. In the future, increased advances are required to enhance egg-quality traits in particular. Egg-related traits are complex traits controlled by multiple quantitative trait loci (QTL) and are influenced by genetic and environmental factors and their interactions [1]. QTL analysis can rate, egg weight, and egg shell weight. Recoquillay et al. [51] found four QTL for the mean egg weight located on chromosomes 1, 3, and 18, two QTL for the number of eggs laid positioned on chromosomes 3 and 18, and two QTL for the age at first egg placed on chromosomes 3 and 19 in a study of medium-density genetic map and QTL analyses of the Japanese quail, in which egg weight as an egg-quality trait and two egg production traits were investigated. Minvielle et al. [52] identified QTL for egg weight, egg number, and age at first egg on chromosome 6 and eggshell weight on chromosomes 1, 5, and 20 using microsatellite markers.
Next-generation sequencing [53] and the construction of high-resolution linkage maps based on SNP markers have been performed in the Japanese quail [6,51]. To the best of our knowledge, there are no published reports on QTL analysis of egg-related traits using RAD-seq in the Japanese quail. However, genetic analysis studies using the RAD-seq method in chickens [54,55], aquaculture [56][57][58], and mammals [59,60] have been reported previously. Hence, the aim of this study was to detect QTL affecting egg-quality and -production traits by constructing a genetic linkage map based on SNP markers using RAD-seq in the Japanese quail.

Experimental Birds
LS and NS quail strains were reared at the Research Farm of Hiroshima University, Japan. F 1 birds were obtained from an intercross of one LS male with three NS females. Subsequently, 138 F 2 females were produced from an intercross between three F 1 males and nine F 1 females. In addition to the above three parents and nine F 1 females, the phenotypic values of 97 parents and 16 F 1 females were used for comparison. Chick management and food supply were in accordance with a previous study [61]. Newly hatched chicks were leg-banded and weighed before being moved to heated brooders, where they were reared until 4 weeks of age. Thereafter, they were housed in individual cages (depth: 15 cm; width: 18 cm; height: 18 cm). Chicks were fed a standard chick diet (22% crude protein (CP); 2900 kcal metabolizable energy (ME) kg −1 ) ad libitum for ages 0 to 4 weeks, followed by a grower diet (17% CP and 2850 kcal ME kg −1 ) from 4 to 16 weeks of age. The birds were maintained under a 24-h lighting photoperiod for 4 weeks, followed by a 14 h:10 h light:dark cycle. They were reared according to the protocol described in the Guidelines for Proper Conduct of Animal Experiments [62].

Traits
A total of 263 female quails were used to assess egg-related traits. Blood samples from all birds were collected using the method described by Kabir et al. [63]. External and internal egg-related traits were measured at two different egg production stages, at the beginning of the egg production period (first stage) and 12 weeks of age (second stage). The two stages are indicated with the subscript letters 1 and 2 in the abbreviations of traits. Five eggs from each quail hen were evaluated at each stage. In general, a total number of 2630 eggs were measured from the parental, F 1 , and F 2 generations. A total of 1000, 250 and 1380 eggs were evaluated at the parental, F 1 , and F 2 generations, respectively. External and internal egg-related traits, including egg weight (EW), egg long axis (ELA), egg short axis (ESA), egg shell strength (ESS), egg shell weight (ESW), egg equator thickness (EET), yolk weight (YW), yolk diameter (YD), yolk color, lightness (L* value) (YC-L*), redness (a* value) (YC-a*), and yellowness (b* value) (YC-b*), and albumen weight (AW), as well as the age at first egg (AFE), total number of laid eggs from maturation up to 16 weeks of age (TLE), and egg production rate (EPR) were used for evaluation in this study. EW, ESW, YW, and AW were measured using a digital balance GX-2000 (A & D Company Ltd., Tokyo, Japan). The lengths of ELA and ESA, and YD were measured using a digital micrometer ABS Digi-Kanon (Nakamura Mfg. Co. Ltd., Tokyo, Japan). For the ESS, an egg shell strength meter (FHK Fujihira Industry Co. Ltd., Tokyo, Japan) was used. EST was measured using an eggshell thickness micrometer (Mitutoyo Industry Co. Ltd., Tokyo, Japan). Yolk color results were obtained using the international colorimetric system of CIELAB. The calibration was based on a black standard with L* = 0 and a white standard with L* = 100. The balanced CIELAB system was determined by three mutually perpendicular axes L*, a*, and b* defined by lightness, redness, and yellowness of yolk color, respectively. The coordinates a* and b* represent the regions of the spectrum with wavelengths corresponding to colors from green (−a) to red (+a) and from blue (−b) to yellow (+b), respectively. The scale of lightness L* represents an interval from 0 (black) to 100 (white). Thus, the complementary color system was based on the differences in the three elementary color pairs: red/green, yellow/blue, and black/white [64,65]. A Chroma meter CR-300 was used to measure the YC-L*, a*, and b* values (Konica Minolta Holdings Inc., Tokyo, Japan).

RAD Library Preparation and Sequencing
Collected blood samples from all birds were used for the extraction of genomic DNA using the phenol-chloroform method and Qiagen kit (DNeasy Blood & Tissue Kits-QIAGEN, Venlo, The Netherlands), according to the manufacturer's protocol. The quality of extracted DNA in each sample was measured using a Qubit 3.0 assay fluorometer (Thermo Fisher Scientific Inc. Waltham, MA, USA), and the DNA was finally adjusted to a concentration of 20 ng/µL for RAD-seq. RAD-seq was performed for the 4 parents and 12 F 1 and 138 F 2 individuals. The RAD library construction procedure followed the methodology previously described by Sakaguchi et al. [66]. The RAD library was sequenced using HiSeq 2500 (Illumina, San Diego, CA, USA) in a 50-bp single-end adapter using EcoRI and BglII. The RAD-seq read data were consigned in the DDBJ Sequence Read Archive (accession no. DRA011153). The RAD-seq reads were trimmed using the Trim_Galore program (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/, accessed on 16 March 2020), and the trimmed reads were mapped onto the Japanese quail genome (GCA_001577835.1 Coturnix japonica 2.-NCBI) using Bowtie2 [67]. The resulting binary sequence alignment/map format (BAM) files processed with SAMtools [68] were used for variant detection. Variant detection was initially performed for F 1 strains. The BAM files of F 1 strains were processed using SAMtools mpileup and varscan2 mpileup2cns [69] with default parameters but changed to min-coverage 5. The variant call format (vcf) files were merged with bcftools [70], and the merged vcf file was further screened using vcftools [71] with the following parameters: minDP 5, min-meanDP 5, maxDP 100, minalleles 2, and max-alleles 2. The screened sites, which were heterozygous for all F 1 strains, were summarized in the position list. With the position list, polymorphisms of all samples, including strains of P 1 , P 2 , F 1 , and F 2 , were called using samtools mpileup and varscan with the above parameters and merged using bcftools. After the polymorphism detection steps, only the GT fields were exported using vcftools and used for further analysis.
The RAD-seq reads that showed different alleles between the two parental strains and genotyped in more than 90% of the F 2 resource population were selected for mapping. The Map Manager QTX b20 application was used to construct the linkage map [72]. The Kosambi map function was used to calculate the genetic distance with linkage criterion p = 0.001 [73].

QTL Analysis
The scanone function of R/qtl [74] was used for QTL analysis using the simple interval mapping method [75]. The genome-wide significance (1% and 5%) and suggestive at (10%) threshold levels were determined by 1000 permutation tests [76]. The significant thresholds for the Z chromosome were estimated using the method described by Broman et al. [77]. The confidence interval from the logarithm of odds (LOD) drop-off method was calculated as 1.8 in this study [78]. The statistical analysis for the calculation of confidence interval, percentage of phenotypic variance, and additive and dominant effects of QTL was conducted using R/qtl based on a previous study [76].

Statistical Analysis
Egg-related data were adjusted for the effects of birth date and dams using the least squares method in R-Project for Statistical Analysis v. 3.6.1. [79]. Egg-related trait comparisons among parental, F 1 , and F 2 generations were performed using one-way analysis of variance (ANOVA), followed by Tukey's HSD test with JMP v. 11.0.2 (SAS Institute Inc., Tokyo, Japan).

Phenotypic Values
The estimated means and standard errors for egg-related traits are shown in Table 1. The LS birds presented the highest values for EW 1 , ELA 1 , ESA 1 , AW 1 , YW 1 , ESW 1 , EW 2 , ESA 2 , AW 2 , and ESW 2 traits. Significant differences were observed for the YC-a* 1 and ESA 2 traits between each generation. The highest values were detected for YC-a* 1 in NS, both stages of YC-L* in F 1 , and YC-b* 1 , YD 2 , and AFE in the F 2 generation. On the contrary, the lowest values of ESA 2 , YD 2 , and YW 2 for NS and YC-a* in both stages and ESS 2 for F 2 generation were observed. No significant difference was observed between the parents and the F 1 generation for ESA 2 , TLE, and EPR. In addition, no significant difference was observed for EW 1 , ELA 1 , ESA 1 , and YW 1 between the NS and both filial generations.

RAD Sequencing and SNP Markers
The Illumina HiSeq 2500 yielded 123,700,031 RAD-seq reads for F 2 birds. After removing the uninformative markers, a total of 25,631 SNP markers were identified. Four hundred twenty-five SNPs were excluded that did not fit the Chi-squared goodness-offit test (p < 0.05). Out of the remaining 25,206 SNPs, 9843 markers were discarded that genotyped at least 90% of all F 2 individuals. After stringent filtration, 14,659 SNPs with a low level of heterozygosity in the F 1 parents were removed. After quality-based filtering, 571 SNPs were discarded from further analyses to reduce missing and wrongly called SNPs. A total of 133 high-fidelity SNP with fixed genotypes in both parents were obtained after excluding those with a deviation for the Mendelian segregation pattern. The QTX Map manager detected 15 unlinked markers that were excluded. Finally, 118 SNP markers were found to be informative between the parental strains.

Linkage Map
A linkage map was successfully drawn using the selected markers in Figure 1. The summary statistics of the constructed genetic map shown in Table 2. Linkage maps contained 118 SNP markers arranged in 13 linkage groups. The chromosome Z in this study was divided into two linkage groups. Each linkage group contained 2 to 32 SNPs. The length of Linkage groups ranged from 0 to 205.3 cM. SNPs covered 776.1 cM of total genetic length with an average spacing of 7.4 cM.

QTL Detection
Nine main-effect QTL were detected on chromosomes 2, 4, 6, 10, 12, and Z (Table 3). Genome-wide significant and suggestive levels for simple interval mapping calculated were 2.963-4.562 in LOD score for the detected QTL. Genome-wide LOD plots for detected QTL are shown in Figures 2-5. The plot of a phenotype against the genotypes at a marker for the identified QTL was drawn in Supplementary Figures S1-S12. Two genome-wide significant and suggestive QTL were found on chromosome 2 for AFE and YC-b* 2 , respectively. QTL for AFE was located at 108.0 cM with 10.4% of phenotypic variance, and YC-b* 2 QTL positioned 123.5 cM with a phenotypic variance of 9.9%. A single genome-wide suggestive QTL was detected for ELA 1 located at 13.3 cM on chromosome 4. The detected QTL on chromosome 6 were located at 54.0 and 69.0 cM, which were responsible for YW 1 and ESS 1 , respectively. Discovered QTL on chromosome 10 affected EW 1 , ESA 1 , AW 1 , and YW 1 traits that ranged in the same flanking markers (0.0-26.0 cM). A genome-wide suggestive QTL was found for ESW 1

Discussion
The present study revealed nine QTL for 10 egg-related traits (EW 1 , ELA 1 , ESA 1 , ESS 1 , YD 1 , AW 1 , YW 1 , ESW 1 , YC-b* 2 , and AFE) in 138 F 2 female Japanese quails. To the best of our knowledge, four studies have conducted QTL analysis on egg-quality and -production traits in the Japanese quail. In these studies, fewer egg-related traits and markers were used for QTL analysis. As the Japanese quail is a model bird in Galliformes, the data obtained in quail can be applied to other bird species, especially in chickens. In addition, the Japanese quail shows close phylogenetic relatedness to chickens with similar length of genome, chromosome number, and morphology of chromosomes. Therefore, a high rate of contiguity, assembly statistics, gene content, chromosomal organization, and synteny conservation exists between the two species [42,80]. Knaga et al. [31] reported QTL for the egg number and egg production rate at positions 36 to 42 cM on chromosome 1 using 30 microsatellite markers. Minvielle et al. [52] reported QTL for egg number and age at first egg positioned at 32 and 34 cM, respectively, on chromosome 6. Recoquillay et al. [51] reported QTL for egg number on chromosomes 3 (225 cM) and 18 (3 cM) and QTL for age at first egg on chromosomes 3 (302 cM) and 19 (4 cM). The results of the current study detected a QTL for AFE located on chromosome 2. However, we did not reveal any QTL for EPR and TLE in the present study. The discrepancy in the number of laid eggs between the results of the present study and those of previous studies may be related to the different durations of the egg collection period. The results of this study are in line with those of Goraga et al. [81], who found QTL for AFE on chromosome 2 in chickens.
The present study found a genome-wide significant QTL for EW at the first stage of egg-laying on chromosome 10 positioned at 21.0 cM. A QTL responsible for EW was reported on chromosome 1 at a position of 4 cM for the Japanese quail [31]. A previous study identified four QTL for the EW trait on chromosomes 1, 3, 18, and 18 positioned at 193, 156, 61, and 62 cM, respectively [51]. At position 0 cM on chromosome 6, a QTL mapped for EW was reported in Japanese quail eggs [52]. The differences in the positions and chromosome numbers of detected QTL for EW may be related to the different periods of weight determination and suggest that other genes exert an effect during the laying period rather than at a later age [31]. In chickens, QTL for the first egg weight were determined on chromosomes 4 and 8 at 213 and 40 cM, respectively [47]. In addition, QTL underlying EW at later ages have been reported on chromosomes 8 and 10 in chickens [47,82].
A genome-wide suggestive QTL was found for ESW at the first stage of egg-laying on chromosome 12 in the current study. In turn, studies of ESW in the Japanese quail have reported QTL on chromosome 1 [31] and chromosomes 1, 5, and 20 [52]. A similar result was found in the case of chickens in which QTL associated with ESW was mapped on chromosome 12 with a confidence interval of 0-71 cM [83]. The discrepancy between the findings of QTL for ESW in the Japanese quail may be due to different times during the production period for determining ESW or the lack of a single preferred way to measure the shell weight. Therefore, it is important to identify QTL that differ over time for the traits of interest.
To date, no QTL have been detected for AW and YW in Japanese quail eggs. This study is the first to report QTL for AW and YW. However, QTL detection for AW and YW has been reported in chickens to be located on chromosomes 2, 3, 4, 7, and 8 for AW and 1-6, 8, 9, 11, 13, 15, 17, 22, 23, 26, 28, and Z chromosomes for YW [12,46,50,82,84,85]. No QTL were reported for AW and YW on chromosome 10 in chickens. In the present study, the detected QTL for YW on chromosome 6 did not overlap with the previous study in chickens [82]. This might be due to the differences in crosses designed for both studies or ages of the two species used, as genetic control for YW is age-dependent. The AW and YW traits are influenced by the age of quails. AW and YW increase with the age of the females. For eggs of the same size, old quails produce larger AW and YW than do young quails. However, in chickens, the genetic correlations among AW and YW at different age points are relatively high [46,85].
To date, there have been no studies detecting QTL for ELA, ESA, ESS, YD, and YC in Japanese quail eggs. This is the first study to detect QTL for ELA 1 , ESA 1 , ESS 1 , YD 1 , and YC-b* 2 traits in the Japanese quail. Nevertheless, there are available QTL for these traits in chickens. Sasaki et al. [86] found a QTL underlying ELA on the same chromosome as that in the current study. The presence of QTL for ESA at the early stage of egg production was reported on chromosomes 10 and Z positioned at 4 and 85 cM, respectively [47]. The positions of the two QTL for ESA are not in line with the positions of QTL on chromosomes 10 and Z for ESA 1 in the current investigation. Similar to the present study, QTL for ESS have been reported on chromosome 6 in chickens [82]. However, Goto et al. [47] found a QTL for ESS on the Z chromosome at the beginning stage of egg production, positioned at 51 cM. In the present study, a QTL for YD 1 was mapped on the Z chromosome. However, no QTL have been reported for yolk index on this chromosome in chickens. Detected QTL for the yolk index in chickens were located on chromosomes 2, 8, and 17 [46,87]. The presence of novel QTL and mapping information is very important for understanding the genetic basis of egg-related traits.
An important biophysical parameter of egg yolk is its color, which plays a very important role in the perception of food and is also a key aspect of food quality. In the current study, one QTL for the yellowness of yolk color at the second stage of egg laying was detected on chromosome 2. QTL for yolk color were also identified on chromosomes 1-4, 8, 9, and 27 in chickens [46,50,88,89]. There is a discrepancy between the findings of this study and those of chickens. It is suggested that the different methods used for the measurement of YC in the current study and the literature on chickens may account for the inconsistent results.
The findings of this investigation suggest that egg-related traits are regulated by many QTL at different growth stages associated with the importance of egg quality and production in breeding programs for poultry industries. These findings may help to develop the Japanese quail QTL database, which can be applied to other Galliformes species, especially chickens. Detected QTL affecting egg-related traits provides necessary information for further molecular studies to improve quantitative traits. It would be desirable to perform a deeper analysis of detected loci using an expanded dataset and sequence information to assess the potential of the found QTL for MAS. Japanese quail can be a valuable resource for the poultry industry sector, and the results of the current study may enable poultry breeders and industries to develop good breeding strategies based on future MAS studies. In conclusion, we determined nine main-effect QTL that affect eggrelated traits using LS and NS Japanese quail strains. This is the first study to report quail QTL information determined using RAD-seq method, including the detection of many trait-regulating QTL that have never been reported before. Thus, these findings will help us understand the genetic basis of egg-related traits and assist breeders in implementing effective selection programs based on MAS.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/genes12050735/s1, Figure S1: Plot of the phenotype against the genotypes for age at first egg QTL detected on chromosome 2, Figure S2: Plot of the phenotype against the genotypes for yolk color-yellowness 2 QTL detected on chromosome 2, Figure S3: Plot of the phenotype against the genotypes for egg long axis 1 QTL detected on chromosome 4, Figure S4: Plot of the phenotype against the genotypes for yolk weight 1 QTL detected on chromosome 6, Figure S5: Plot of the phenotype against the genotypes for egg shell strength 1 QTL detected on chromosome 6, Figure S6: Plot of the phenotype against the genotypes for egg weight 1 QTL detected on chromosome 10, Figure S7: Plot of the phenotype against the genotypes for egg short axis 1 QTL detected on chromosome 10, Figure S8: Plot of the phenotype against the genotypes for albumen weight 1 QTL detected on chromosome 10, Figure S9: Plot of the phenotype against the genotypes for yolk weight 1 QTL detected on chromosome 10, Figure S10: Plot of the phenotype against the genotypes for egg shell weight 1 QTL detected on chromosome 12, Figure S11: Plot of the phenotype against the genotypes for egg short axis 1 QTL detected on Z chromosome, Figure S12: Plot of the phenotype against the genotypes for yolk diameter 1 QTL detected on Z chromosome.