Dec 2, 2023

Abstract This study was designed to perform an association analysis and identify SNP markers associated with production traits of Japanese quail using restriction-site-associated DNA sequencing. Weekly body weight data from 805 quail were collected from hatching to 16 weeks of age. A total number of 3990 eggs obtained from 399 female quail were used to assess egg quality traits. Egg-related traits were measured at the beginning of egg production (first stage) and at 12 weeks of age (second stage). Five eggs were analyzed at each stage. Traits, such as egg weight, egg length and short axes, eggshell strength and weight, egg equator thickness, yolk weight, diameter, and colour, albumen weight, age of first egg, total number of laid eggs, and egg production rate, were assessed. A total of 383 SNPs and 1151 associations as well as 734 SNPs and 1442 associations were identified in relation to quail production traits using general linear model (GLM) and mixed linear model (MLM) approaches, respectively. The GLM-identified SNPs were located on chromosomes 1–13, 15, 17–20, 24, 26–28, and Z, underlying phenotypic traits, except for egg and albumen weight at the first stage and yolk yellowness at the second stage. The MLM-identified SNPs were positioned on defined chromosomes associated with phenotypic traits except for the egg long axis at the second stage of egg production. Finally, 35 speculated genes were identified as candidate genes for the targeted traits based on their nearest positions. Our findings provide a deeper understanding and allow a more precise genetic improvement of production traits of Galliformes, particularly in Japanese quail. Introduction Japanese quail is a model bird from the Galliformes order that is raised for production and biological research purposes. The Japanese quail holds significant importance as a bird species owing to its compact size, efficient productivity, quick generational turnover, and early sexual maturity, typically around 6 weeks of age 1 . Quail production traits are important characteristics in the poultry industry, and quantitative trait loci (QTLs) associated with these traits have been identified 2 , 3 . Molecular markers have become a key prerequisite for association mapping and are of increasing importance in animal breeding and genetics. Association analysis is an advantageous technology that identifies QTLs underlying phenotypic traits 4 and provides a link for breeders to select based on genetic information 5 . Among the molecular markers, single nucleotide polymorphism (SNP) markers are widely used 6 , 7 , 8 , 9 due to their abundance in any genome and the cost-efficient identification methods 10 . Genotyping by sequencing (GBS) is a next-generation sequencing (NGS) platform that opens new possibilities for SNP marker identification and can be used in a simple, highly multiplexed system for constructing reduced representation libraries involving inexpensive barcoding, reduced sample handling, no size fractionation, and which requires fewer polymerase chain reaction and purification steps 11 . Restriction-site associated DNA sequencing (RAD-seq) 12 , 13 , 14 is a GBS that can identify, verify, and score thousands of SNPs simultaneously, reduce complexity across genomes, deliver high-resolution population genomic data, and is convenient for non-model species at a reasonable cost 15 . Existing genotyping platforms represent efficient tools for association studies to increase the output of marker-assisted selection (MAS). Association analysis is a viable approach for linking phenotypes and genotypes in poultry genetics as well as dissecting complex traits. This type of analysis can boost the information obtained from QTL studies through MAS implementation. Quail production traits are complex traits that are controlled by QTLs and can be affected by environmental factors. The QTLs related to production traits of Japanese quail have been previously identified. Reference 2 identified 22 QTLs for body weight, egg weight, number of eggs laid, and age of the first egg using SNP markers. References 16 , 17 discovered four QTLs for body weight and nine QTLs for egg-related traits using the RAD-seq method. Moreover, microsatellite markers have been used to detect QTLs for traits associated with body weight traits 18 , 19 , 20 , 21 , 22 , 23 , growth and egg production 3 , 24 , and egg-laying curves 25 . Previous studies have also investigated the QTLs underlying features other than production traits of Japanese quail 26 , 27 , 28 , 29 , 30 . Although molecular breeding helps breeders select the production traits expressed in Japanese quail, association analysis studies in poultry have mostly focused on chickens rather than Japanese quail. In recent years, genome-wide association studies (GWASs) have been widely used in studies of different chicken breeds using a high-density genotyping platform to identify SNPs and candidate genes associated with production traits 31 , 32 , 33 , 34 , 35 , 36 , 37 . Numerous studies used GWASs to identify SNPs that control the candidate genes responsible for body weight traits of chickens 38 , 39 , 40 , 41 . Similar studies have been conducted for identifying SNPs associated with chicken egg-related traits 42 , 43 , 44 . References 45 , 46 , 47 , 48 performed GWASs for body composition and meat quality traits of chickens. Studies of association analysis in Japanese quail are limited and mainly concentrated on the association between polymorphisms of gonadotropin-releasing hormone genes and growth traits 49 , prolactin receptor genes and growth traits 50 , and gut microbial architecture of efficiency traits 51 . Existing reports on the genetic map of Japanese quail have been integrated and aligned with assembled chicken sequence data 52 . As Japanese quail shows close phylogenetic relatedness to chickens, having a similar genome length (1.2 × 109 base pairs), chromosome number (2n = 78), and homology of chromosome morphology, a high rate of synteny conservation is expected between the two species 53 . In addition to studies on the quail genome 54 , 55 , reports on the chicken genome will help breeders investigate genetic analysis and identify SNP markers associated with the targeted traits within the poultry industry. However, to date, there are no published reports on association studies of production traits of Japanese quail using GLM and MLM models with RAD-seq data. Therefore, a study pursuing the identification of SNP markers associated with production traits will provide breeders with a useful tool to assist in selecting high-production lines, specifically in Japanese quail breeding programs and poultry in general. Here, we performed an association analysis and identified SNP markers associated with production traits of Japanese quail using the RAD-seq method. We expect these results to help improve poultry breeding programs and increase production levels through molecular breeding using SNP markers with MAS. Materials and methods Experimental design Animal care, experimental protocols, and blood sample collections were approved and conducted in accordance with the Rules on Experimental Animals and Animal Experiments at Hiroshima University, Graduate School of Integrated Sciences for Life, Laboratory of Animal Breeding and Genetics (Approval No. C20-15) and the protocol described in the Guidelines for Proper Conduct of Animal Experiments, Science Council, Japan. https://www.scj.go.jp/ja/info/kohyo/pdf/kohyo-20-k16-2e.pdf . Also, we confirm that all methods have been conducted in adherence to the ARRIVE guidelines ( https://arriveguidelines.org ). Experimental birds were reared at the Research Farm of Hiroshima University, Japan. Large- and normal-sized (LS and NS) Japanese quail strains were selected as parental strains. The LS Japanese quail is known for its high body weight, reaching around 200 g in females and 170 g in males upon maturity, making it suitable for meat production purposes. Conversely, the NS strain of Japanese quail exhibits a normal body weight, approximately 130 g in females and 100 g in males at the age of maturity, and has been primarily used for egg production 1 . These two strains originate from genetically distinct backgrounds concerning the specific traits of interest, which could potentially aid in accurately identifying targeted SNPs. One hundred quails from each of the two strains, consisting of 50% males and 50% females, were reared as parents. The cross involved an LS male paired with three NS females, and an NS male paired with three LS females, resulting in 100 F1 progeny in a reciprocal cross (50 males and 50 females). Subsequently, a total of 505 F2 birds (256 males and 249 females) were generated from a reciprocal cross between six F1 males and 18 F1 females. Newly hatched chicks were instantly pedigree leg-banded and weighed before being moved to heated brooders, where they were reared until four weeks of age. Thereafter, quail were housed in individual steel wire meshed cages (15 cm deep, 18 cm wide, and 18 cm high) equipped with feeders and drinkers. Birds were fed ad libitum a standard starter diet (22% crude protein (CP); 2900 kcal metabolizable energy (ME)/kg−1) from 0 to 4 weeks of age and then a grower diet (17% CP and 2850 kcal ME/kg−1) from 4 to 16 weeks of age. Quail were reared under a 24-h light photoperiod for 4 weeks, followed by a 14 h:10 h light:dark cycle. Trait records Weekly body weight data of 805 quail (399 females and 406 males) from the parental, F1 and F2 generations were collected from hatching until 16 weeks of age (WK 0–16). A total number of 3990 eggs obtained from 399 females were used to assess egg-related traits. Egg-related traits were measured at two different egg production stages: at the beginning of the egg production stage (first stage) and at 12 weeks of age (second stage). The two stages are indicated with subscript letters ‘1’ and ‘2’ in the abbreviations of traits. The first five eggs from each female were evaluated at each developmental stage. External and internal egg-related traits, including egg weight (EW), egg long axis (ELA), egg short axis (ESA), eggshell strength (ESS), eggshell weight (ESW), egg equator thickness (EET), yolk weight (YW), yolk diameter (YD), yolk colour, lightness (L* value) (YC-L*), redness (a* value) (YC-a*), yellowness (b* value) (YC-b*), albumen weight (AW), 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 evaluated using the methodology described in our previous publication 17 . The average data from parental, F1 and F2 generations were subjected to one-way analysis of variance (ANOVA), followed by Tukey’s HSD test using JMP v. 11.0.2 (SAS Institute Inc., Tokyo, Japan). DNA collection and RAD library sequencing Blood samples from all birds were collected using the method described by 56 , and genomic DNA was extracted from each collected sample using phenol–chloroform and DNeasy Blood & Tissue Kits (Qiagen, Venlo, Netherlands) according to the manufacturer's protocol. DNA was quantified using a Qubit 3.0 assay fluorometer (Thermo Fisher Scientific Inc. Waltham, MA, USA). DNA concentrations were adjusted to 20 ng/µL to ensure consistent measurements and were subsequently used for library preparation. Birds used for matings along with 505 F2 quail were used for RAD-seq. RAD libraries were prepared according to the method described by Ref. 57 . A 50-bp single-end adapter using EcoRI and BglII enzymes was performed on an Illumina HiSeq 2500 (Illumina, San Diego, CA, USA) for RAD library sequencing. RAD-seq read data were deposited in the DDBJ Sequence Read Archive (accession no. DRA011153) https://ddbj.nig.ac.jp/search . The RAD-seq reads were trimmed using the TrimGalore 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 with default settings 58 . Further, the reads were converted to binary sequence alignment/map format (BAM) files using SAMtools 59 . Variant detection was initially performed for the F1 generation. The BAM files of F1 lines were processed using the SAMtools mpileup and varscan2 mpileup2cns 60 with default parameters and then changed to min-coverage 5. The variant call format (vcf) files were merged with bcftools 61 , and the merged file was further screened using vcftools 62 with the following parameters: minDP 5, min-meanDP 5, maxDP 100, min-alleles 2, max-alleles 2, and max-missing 1. The screened sites that were heterozygous for all F1 populations are summarised in the position list. Within the position list, polymorphisms of all samples, including birds from P1, P2, F1, and F2, were named using the SAMtools mpileup and varscan with the aforementioned parameters and merged using bcftools. After the polymorphism detection steps, only the GT fields were exported using vcftools and used for further analyses. Association analysis Association analysis was performed using general linear model (GLM) and mixed linear model (MLM) approaches 63 implemented by TASSEL v. 5 64 . GLM evaluation incorporated population structure (Principal component analysis, (PCA)) and MLM used both population structure (PCA) and kinship matrix in the association analysis. The PCA and covariance analysis were performed using the method reported by Ref. 65 to describe the population structure. A kinship matrix was produced using the Kinship tool with the Scald_IBS method to avoid spurious associations due to relatedness and population structure. In the "Scald_IBS" method, genotypes were encoded based on the count of the alleles at that locus. By utilizing this encoded data, the kinship matrix was estimated. The efficient mixed-model association (EMMA) algorithm 66 and population parameter previously determined (P3D) variance component estimation were implemented to decrease the computing time for a large dataset 67 . Statistically significant associations with p-values < Pthreshold = 0.01/N, where N is the total number of SNP markers, were identified, and a standard correction was performed by applying a Bonferroni procedure at a minimum value of 3.88806E−36 for GLM and 7.03385E−22 for MLM 68 . The coefficient of determination (R2) was subsequently calculated to quantify the extent to which genetic markers explained the phenotypic variance. The significant value and the marker effect for each SNP were exported, and Manhattan as well as quantile–quantile (QQ) plots were generated in the R Project for Statistical Analysis v. 3.6.1 69 using the qqman package 70 . The data were subjected to a correction procedure to account for the influence of fixed effects, specifically birth date, sex, and dame. This adjustment was performed using the least squares method in R-Project for Statistical Analysis v. 3.6.1 69 . Ethical approval Animal care, experimental protocols, and blood sample collections were approved and conducted in accordance with the Rules on Experimental Animals and Animal Experiments at Hiroshima University, Graduate School of Integrated Sciences for Life, Laboratory of Animal Breeding and Genetics (Approval No. C20-15) and the protocol described in the Guidelines for Proper Conduct of Animal Experiments, Science Council, Japan https://www.scj.go.jp/ja/info/kohyo/pdf/kohyo-20-k16-2e.pdf . Also, we confirm that all methods have been conducted in adherence to the ARRIVE guidelines ( https://arriveguidelines.org ). Results Filial generations obtained from reciprocal matings are presented separately. The body weights of female and male birds in the LS and NS strains and their F1 and F2 generations are shown in Tables 1 and 2 , respectively. In females, LS F1 (LS♀ ✕ NS♂) quail showed the highest body weight at hatching (8.19 and 8.54 g, respectively). The lowest body weight was detected in NS quail from hatching to 16 weeks of age. No difference was found between F2 generations in WK 0–2. From 4 to 16 weeks of age, the LS strain presented the highest body weight compared to that of the other groups. In males, the LS F1 (LS♀ ✕ NS♂) quail showed the highest body weight at birth. However, no difference was observed between the weights of the other quail groups, except for NS. Moreover, no significant differences were observed between the F1 generations in WK 1–16. The NS quail showed the lowest body weight during the entire period analysed. The LS strain exhibited the highest body weight from 4 to 16 weeks of age compared to those observed at all generations. Table 1 Body weight (mean ± standard error) of LS and NS female birds and their F1 and F2 hybrids from 0 to 16 weeks of age. Table 3 shows the means and standard errors for egg-related traits in the LS and NS strains of Japanese quail and their F1 and F2 hybrids in the first and second egg-laying stages. The LS strain EW, ESA, ESW, and AW exhibited the highest value among all groups in both egg-laying stages. NS females showed the highest values for the YC-a*1 and YC-b*2 traits across generations. F2 (LS♂ ✕ NS♀) females showed the lowest value in YC-a*1 (2.64) compared to the those observed in the other groups. No difference was observed in EET1 between the parental and filial generations. No significant differences were observed in ELA1 between NS and both filial generation females. Both F2 groups presented the lowest ESS1,2 among all groups, whereas no differences were observed between the parental and F1 quail in the first and second stage of egg production. F2 (LS♂ ✕ NS♀) birds started egg laying at a later age (49.01 days) than the other groups. F1 birds laid more eggs than the parental and F2 generations. Moreover, egg production in F2 was lower than that in the parental generation. Table 3 Means and standard errors for egg-related traits in LS and NS strains of Japanese quail and their F1 and F2 hybrids. Illumina HiSeq 2500 produced 733,233,057 RAD-seq reads for the treated samples. A total of 25,631 SNPs that passed the quality control were used for the GLM and MLM analyses. This subset of SNP genotypes was used to produce a kinship matrix for MLM analysis. The same subset of SNP markers was also used to perform the PCA in this study. General and mixed linear models produced 1,941 and 2,986 associations, respectively, after applying statistically significant loci with Bonferroni corrections. Statistically significant associations with p-values lower than the threshold were excluded, and 383 SNPs for 1,151 associations (Supplementary File S1 ) and 734 SNPs for 1,442 associations (Supplementary File S2 ) were identified to be significantly associated with production traits in the GLM and MLM approaches, respectively. GLM identified SNPs that were located on chromosomes 1–13, 15, 17–20, 24, 26–28, and Z, underlying all analysed traits except EW1, AW1, and YC-b*2 (Table 4 ). The maximum number of SNPs and associations was identified on the Z chromosome (122 and 611, respectively), followed by chromosome 1 (61 and 95, respectively) and chromosome 4 (52 and 92, respectively), while the minimum number of SNPs and associations (1 and 1, respectively) was found on chromosomes 13, 20, and 28 underlying WK 0, ESS1, and WK 3, respectively. A total of 1065 and 86 associations were detected for body weight and egg-related traits, respectively (Table 4 ). No associations for body weight traits were identified on chromosomes 12 (5,472,506 bp) and 20 (11,476,337 bp). Additionally, no associations for egg-related traits were detected on chromosomes 13, 24, and 26–28. The percentage of phenotypic variance explained by marker (R2) detected by the associations ranged between 0.003 and 0.097 using the GLM approach. Figures  1 and 2 represent the Manhattan and QQ plots for ESS2, respectively, wherein 14 associations were found on chromosomes 1, 2, 4, 5, 11, and Z. The Manhattan and QQ plots for the remaining traits from the GLM are shown in Supplementary File S4 . Table 4 Summary of chromosomes contained SNPs associations in production traits in GLM approach. Figure 1 The Manhattan plot shows the association of SNP markers with ESS2 in the GLM. Each dot represents an SNP. The figure illustrates the level of statistical significance (y-axis) as measured by the negative log of the corresponding p-value for each SNP. Each SNPs type is indicated by dots of different colours, which are arranged by chromosomal location (x-axis). The horizontal red line indicates the threshold of 5% Bonferroni genome-wide significance (p = − log10 (5e−08)), and the underlined blue line presents a genome-wide suggestive at (p = − log10 (0.3e−05)). As shown in Table 5 , MLM recognised SNPs that were positioned on all chromosomes except for chromosome 16 and associated with all traits except for ELA2. The maximum number of associations was detected on the Z chromosome (N = 516), followed by chromosomes 1 (N = 208), 3 (N = 105), and 2 (N = 96), while the minimum number of associations was observed on chromosome 25 (N = 1) for YC-L*2 and chromosome 28 (N = 1) for egg production rate. However, a maximum number of 149, 97, 76, and 69 SNPs were distributed on chromosomes 1, Z, 2, and 3, respectively. Following MLM and Kinship analyses, 739 and 703 associations were observed for body weight and egg-related traits, respectively. No SNP markers on chromosomes 11–12, 19–23, 25, and 28 were associated with body weight. The percentage of phenotypic variance detected by the associations ranged between 0.0004 and 0.1401 in the MLM analysis. The Manhattan and QQ plots for WK 1 are shown in Figs. 3 and 4 , respectively. For WK 1 trait, four associations were detected on chromosomes 1 (N = 2), 3 (N = 1), and Z (N = 1). MLM plots for each Manhattan and QQ (Supplementary File S5 ) were plotted for the associated body weight and egg-related traits. Table 5 Summary of chromosomes contained SNPs associations in production traits in MLM approach. Figure 3 The Manhattan plot shows the association of SNP markers with body weight at WK 1 in the MLM. Each dot represents an SNP. The figure illustrates the level of statistical significance (y-axis) as measured by the negative log of the corresponding p-value for each SNP. Each SNPs type is indicated by dots of different colours, which are arranged by chromosomal location (x-axis). The horizontal red line indicates the threshold of 5% Bonferroni genome-wide significance (p = − log10 (5e−08)), and the underlined blue line presents a genome-wide suggestive at (p = − log10 (0.3e−05)). The combination of association markers detected in both the GLM and MLM analyses is explained in Table 6 . A total of 191 SNPs associated with body weight and egg-related traits were shared in both models. Markers associated with all body weight traits were identified in the GLM and MLM analyses. However, no markers were found to be associated with EW1, ELA2, ESW1, EET1, YC-b*2, AW2 in both models, and EW2, ESA1,2, YD2, YC-L*1, YC-a*1,2, AW1, TLE, and EPR based on the GLM. In the GLM, 18 SNP markers were identified on chromosome 1 for ELA1 (N = 1), ESW2 (N = 4), EET2 (N = 1), YW1,2 (N = 1,1), YD1 (N = 1); chromosome 2 for ESS2 (N = 2); chromosome 4 for AFE (N = 1); chromosome 5 for ESS2 (N = 1); chromosome 9 for ESW2 (N = 2); chromosome 12 for ESS1 (N = 1); chromosome 19 for YC-L*2 (N = 1); and the Z chromosome for YC-L*1 (N = 1). Alternatively, in the MLM, 45 SNPs were observed on chromosome 1 for EW2 (N = 3), ESW2 (N = 5), and YW1,2 (N = 1,1); chromosome 2 for ESS2 (N = 2) and TLE (N = 1); chromosome 3 for TLE (N = 3) and EPR (N = 3); chromosome 5 for ELA1 (N = 1) and ESS2 (N = 1); chromosome 7 for YC-L*1 (N = 2); chromosome 8 for AW1 (N = 1); chromosome 9 for EW2 (N = 2), ESA2 (N = 1), ESW2 (N = 1), EET2 (N = 2), YW2 (N = 1), and YC-a*1 (N = 1); chromosome 12 for ESA1 (N = 1); chromosome 15 for YC-b*1 (N = 1); chromosome 19 for AFE (N = 1); chromosome 27 for TLE (N = 1) and EPR (N = 1); and Z chromosome for ELA1 (N = 3), ESA1 (N = 1), YC-a*2 (N = 1), and TLE (N = 3). Among the detected SNPs, 66 markers were associated with the same traits in both models across all chromosomes. Table 6 SNP markers found in both GLM and MLM approaches. Discussion The evaluation of 44 phenotypic traits from 567 birds in the F2 revealed 383 SNPs with 1151 associations and 734 SNPs with 1442 associations for GLM and MLM approaches, respectively. Our identified SNPs were associated with all targeted traits except EW1, AW1, and YC-b*2 in the GLM and ELA2 in the MLM analysis. To the best of our knowledge, no such association analysis of production traits in Japanese quail have been reported using GLM and MLM approaches with RAD-seq data. However, several studies on QTL analysis of production traits in Japanese quail have been reported. Reference 20 recognised QTLs underlying body weight at hatching and at four weeks of age on chromosome 1 at 12–13 cM. Here, we identified associations on the same chromosome at 11.9 cM for GLM and 12.38 cM for MLM corresponding for WK 0 and TLE traits, respectively. Reference 21 reported QTL associated with body weight from hatching to three weeks of age on chromosome 3 between 27 and 30 cM. The results of the present study revealed trait associations for SNPs located between 28.73 and 37.74 cM on chromosome 3 in both GLM and MLM analyses. Moreover, using GLM and MLM approaches, we showed that associations underlying body weight traits were mostly located on chromosomes 1–10 and Z. Associations affecting hatching weight were identified on chromosomes 1–9, 13, 17, 27, and Z for GLM, and 1, 3, 7, 9, 13, and 27 for MLM approaches. Detected SNPs located at 20.16 cM on chromosome 7 were associated with hatching weight in both GLM and MLM analyses, and these results were consistent with the previously identified location of QTLs associated with this trait on chromosome 5 (19–20 cM) 22 . The detected QTL affecting liveweight measurement at week 5 has been shown to be located on chromosome 1 in a study investigating the genetic mapping of QTLs affecting body weight 18 , which is in agreement with the SNPs associated with body weight and located on chromosome 1 in our GLM and MLM analyses. Reference 3 identified a QTL for body weight at 4 weeks of age at the centromere of chromosome 2 which is consistent with the locus detected for this trait on chromosome 2 in both association analyses performed in the present study. The study detected three QTLs associated with body weight at 1, 4, and 6 weeks of age, located in the initial region of chromosome 2 between 0 and 15 cM. Reference 2 reported putative pleiotropic loci on chromosome 3 (52.6–56.7 Mb) affecting both weight and egg traits, where a QTL for egg weight co-localised with a QTL for body weight at 65 days of age. Interestingly, our results indicated an association in the same chromosome (56.08 Mb) controlling WK 1 in both the GLM and MLM analyses. This candidate region could underlie the genetic correlation already observed in quail between body weight and egg traits 71 , 72 , 73 . QTLs for early and late growth stages (5 and 70 weeks of age) have been previously detected on chromosome 1, positioned 18–19 cM 24 . In our MLM analysis, YW associations were identified on chromosome 1 between 17.68 and 21.16 cM. It is likely that the QTLs for BW and YW located on chromosome 1 represent a single gene. In addition to QTLs for body weight traits, Ref. 24 identified QTLs for eggshell weight on chromosomes 1 (191 cM), 5 (12 cM), and 20 (21 cM); egg weight, egg number until the age of 69 weeks, and age of first egg on chromosome 6 positioned at 0, 32, and 34 cM, respectively. The results of the present study revealed associations for YC-L*2 and YC-a*2 on chromosome 1 (190.34 cM) as well as TLE and EPR on chromosome 6 (31.8 cM) in GLM. Alternatively, we identified EET2 QTLs on chromosome 5 (12.02 cM), YW2 on chromosome 6 (0 cM), and YW2 as well as YD2 on chromosome 20 (20.35 cM) in the MLM. These findings are consistent with the results presented by Ref. 24 . In turn, Ref. 3 recognised QTLs for the number of eggs laid and the egg production rate on chromosome 1 between 36 and 42 cM. The differences among the results of these studies may be explained by the different durations of egg production. In our previous report, we identified a QTL for growth-related traits on chromosome 1 using the same Japanese quail strains 16 . Similarly, in the present study, two associations were found on the same locus (150,370,977 bp) on chromosome 1 for BW 4 and 5. In addition, we detected QTLs for ESA1 and YD1 on the Z chromosome between 22,757,726 and 31,279,475 bp 17 . Associations underlying WK 9 in GLM, EET1 in MLM, and WK 7–16 in both approaches were identified within the same chromosome and at the same positions. This was expected as positive correlations have been observed between the related egg trait and body weight 73 . Only chromosomal regions were explored in the above-mentioned studies on the Japanese quails; therefore, future studies focusing on the association between SNPs and candidate genes would be of great value for improving QTL resolution in Galliformes. GWASs are widely used to distinguish SNPs associated with production traits in chickens 74 , 75 . In the first study of its kind on Japanese quail, we used association analysis with GLM and MLM, based on RAD-seq SNPs, to examine production traits like body weight and egg-related traits. Comparative genomic studies of Japanese quail and chickens based on cytogenetics 53 , orthologous genes 76 , and linkage analysis 52 have shown a high rate of synteny-conserved karyotypes and genomic structure between these species. Therefore, the GWAS results from chickens were used to identify potential candidate genes in the present study. However, due to differences in growth conditions between chickens and quail, body weight at the early growth stage, sexual maturity age, and maximum growth rate were considered in both species. The detection of associations related to body weight during juvenile age is likely due to the strong correlation between body weight traits during early growth stages. It is anticipated that numerous genes control body weight during this period because growth is a highly intricate trait influenced by multiple loci that affect appetite, feed intake, body composition, nutrient utilization, and physical activity. In the present study, six shared SNPs between GLM and MLM approaches were identified on chromosomes 1 (147.6 and 167.7 Mb), 3 (29.5 and 40.4 Mb), 13 (13.2 Mb), and 27 (1.9 Mb) associated with body weight at hatching, a trait that is strongly influenced by maternal effects 71 . Reference 77 identified a chicken QTL on chromosome 27 in the same region as the one discovered in the present study. This region on chromosome 27 contains two genes, LIM and SH3 protein 1 (LASP1) and phosphatidylinositol-5-phosphate 4-kinase, type II, beta (PIP4K2B), which regulate the decrease in mice body weight 78 , 79 . Associations between gene expression at the initial ages have also been reported in chickens. Reference 80 identified a genomic region (169.8–175.3 Mb) on chromosome 1 associated with body weight traits. In a study of genetic dissection of growth traits in a Chinese indigenous breed with a commercial broiler cross 81 , the major QTL of body weight was mapped to the end of chromosome 1 (173.7 Mb). Here, associations underlying early growth stages were located on chromosome 1 between 168.74 and 172.99 Mb as indicated by both GLM and MLM analyses. Retinoblastoma 1 (RB1) was determined to be the significant gene in this position using a haplotype approach 82 . Maximum growth rate and body weight at sexual maturity and later ages are similar traits that describe mature body weight and have a high genetic correlation of nearly 1 83 . The sexual maturity in quail is between 6 and 7 weeks of age, while that in chickens is between 16 and 24 weeks of age, depending on the breed 84 . We considered quail’s SNPs associated with body weight around the age of the birds in our study for comparison with the information obtained from chicken studies. The genetic regulation of body weight at maturity and later ages is a multifaceted process involving a network of genes that interact with various physiological systems, hormones, and environmental factors. The complexity of these interactions contributes to the expectation that numerous genes are involved in shaping body weight traits in quails. In a GWAS of carcass traits of Jinghai Yellow chickens 47 , five SNPs were detected for foot weight on chromosome 4 located between 75.54 and 75.67 Mb; thus, indicating that this region is located within 2.08–2.38 Mb away from genes, such as family with sequence similarity 184 member B (FAM184B), quinoid dihydropteridine reductase (QDPR), and LIM-domain binding factor 2 (LDB2). In addition to these genes, Ref. 40 indicated that the F box and leucine-rich repeat protein 5 (FBXL5) has a significant influence on chicken growth traits and important biological functions. They identified this gene to be located on chromosome 4 from 72.9 to 77.9 Mb and to be strongly associated with body weight for weeks 6–16. QDPR and LDB2 were associated with shank circumference 32 and body weight in Beijing You chickens 38 . Reference 85 detected two SNPs that were associated with the body weight at the time of oviposition and are located separately at 78.8 Mb and 78.9 Mb on chromosome 4 within the FAM184B gene. The non-SMC condensing 1 complex, subunit G (NCAPG) gene, which is located within 0.1 kb downstream of the FAM184B gene, is a well-known candidate gene for body frame and carcase traits in cattle, and is considered to modulate the body and carcase weight of various breeds 86 , 87 . At 73.9 and 76.1 Mb, two SNPs associated with ESS2 and ESS1 on chromosome 4, respectively, were identified by both GLM and MLM analyses. Since ESS is positively correlated with body weight 84 , the influence of these genes on body weight is of considerable importance. Moreover, Ref. 88 reported the importance of the myostatin (MSTN) gene located on chromosome 7 on chicken body weight at 112 days of age. Similarly, we detected numerous associations underlying body weight at maturity and later ages on chromosome 7 at 22.65 Mb in both GLM and MLM analyses. The detection of QTL at later ages has also been reported in chickens 89 . The growth hormone-releasing hormone receptor (GHRHR) gene was mapped at position 1.7 cM and polymorphism analysis of this gene revealed three SNPs in the promoter region, which had a significant effect on body weight at 7, 9, 11, 13, and 17 weeks of age 90 . These findings are in line with our SNP association results. In MLM analysis, we identified 10 SNPs that explain over 10% of the phenotypic variance. This underscores their genetic importance and supports our SNP association findings. However, further investigation is needed to uncover the functional roles of these genes, shedding light on their contribution to the observed phenotypic variations at 10–16 weeks of age. In addition to body weight, egg quality and production traits are the major selection criteria for poultry breeding. The findings of the association study on the TLE and AFE in this study hold significant implications for egg-type poultry breeding. The results have identified specific associations associated with variations in egg production traits, offering valuable insights into the genetic basis of egg production. However, the results of egg production may vary across different studies, and this variation may be attributed to the diverse durations of the egg collection period. In a GWAS investigating egg production and quality traits in chickens, Ref. 91 found the most significant SNP associated with egg number located within intron 12 of the growth factor receptor-bound protein 14 (GRB14) gene. GRB14 mRNA is highly expressed in the ovary, liver, kidney, and skeletal muscles of humans and mammals 92 , 93 . Similarly, our MLM suggested that an SNP associated with TLE is located on the same chromosome and position (chromosome 7; 21.46 Mb). It has been speculated that GRB14 influences egg production in layers 91 . The age of first egg is an important indicator of sexual maturation in female birds and is influenced by genetic and environmental factors 94 . Here, we detected 7 associations that might influence AFE, and they were located on chromosomes 1–4 and Z based on the GLM, and 18 associations positioned on chromosomes 1–3, 5, 8, 12, 19, and 21 based on the MLM analysis. A SNP in intron 2 of the odd oz/ten-m homolog 2 (ODZ2) gene has been previously identified to be significantly associated with chicken AFE 91 . This gene is expressed in the developing chicken brain and may affect the sexual maturity 95 . Additionally 96 , demonstrated that the gonadotropin-releasing hormone I (GnRH-I) gene is associated with chicken AFE. GnRH stimulates the synthesis and secretion of gonadotrophins, which induce steroidogenesis in the gonads, culminating in ovarian follicle growth and ovulation for egg production 97 . Reference 98 identified two novel forkhead box L2 (FOXL2) and growth differentiation factor-9 (GDF9) genes associated with egg production in the Chinese Dagu chickens. These polymorphisms play a critical role in the regulation of ovarian development in hens. This result supports our findings in the present study, in which 53 associations were detected for EPR by both the GLM and MLM analyses. Our association study focused on external egg quality traits. Eggshell quality is a major concern in poultry breeding due to its implications for both reproductive performance and human consumption. In this study, we successfully identified several associations associated with external egg traits in both egg production stages. These findings shed light on the genetic basis of egg quality and provide valuable information for targeted breeding efforts aimed at enhancing egg production and ensuring high-quality eggs for consumers. A number of QTL regions for egg weight, egg length, and eggshell characteristics have been already identified 17 . In the present study, 49 associations were identified for EW1,2 by both the GLM and MLM analyses. The associations for EW on chromosome 4 between 35.49 and 65.59 Mb are in line with the SNP for the same trait and chromosome at 49.28 Mb located in the shroom family member 3 (SHROOM3) gene 99 . Shroom3 is a cytoskeletal protein involved in regulating cell shape (arrangement and remodelling) in certain tissues 100 . Among the external egg traits, eggshell characteristics are the most important and play a significant role in the reproductive performance and human consumption. Reference 91 identified three SNPs for chicken ESW on chromosomes 2 (86.11 Mb), 3 (110.09 Mb), and 11 (9.59 Mb), and two SNPs for EST on chromosome 1 (171.22 and 179.35 Mb). Moreover, they introduced a gene for each detected SNP: polypeptide N-acetylgalactosaminyltransferase 1 (GALNT1), BLK, zinc finger protein 536 (ZNF536), ENOX1, and LOC18918. These findings support the results of our study, in which two SNPs were detected for ESS2 on chromosomes 1 (172.99 Mb) and 11 (14.96 Mb) based on the GLM, and 12 associations were found for ESS1,2, ESW1, and EET1,2 on chromosomes 1 (171.51–174.10 Mb), 2 (84.95–89.81 Mb), and 11 (5.16–14.70 Mb) based on the MLM analysis. The sodium channel (SCNN1) gene family is expressed in the active uterus during eggshell mineralisation and plays an essential role in eggshell formation 101 . Polymorphisms of eggshell organic matrix genes were considered to be related to eggshell thickness, eggshell strength, and dynamic stiffness 102 . SCNN1 has four family members, SCNN1a, SCNN1b, SCNN1g, and SCNN1d, which affect eggshell formation 103 . A study investigating eggshell quality traits identified SCNN1a, SCNN1b, SCNN1d, and SCNN1g genes on chromosomes 1 (80.03–80.04 Mb), 14 (7.002–7.01 Mb), 21 (2.43–2.44 Mb), and 14 (7.01–7.02 Mb), respectively, and associated them with egg weight, eggshell weight, eggshell percentage, eggshell strength, and eggshell thickness 101 . Our MLM analysis suggested associations for ESS2 (chromosome 1, 78.06 Mb), ESW1 (chromosome 14, 8.36 Mb), and EW1 (chromosome 21, 2.58 Mb), located near the same chromosome, as stated by Ref. 101 . In addition, an association study identified three candidate genes, phosphatidylinositol-4-phosphate 3-kinase catalytic subunit type 2 gamma (PIK3C2G), inositol 1,4,5-trisphosphate receptor type 2 (ITPR2), and non-SMC condensin I complex subunit G (NCAPG) to be implicated in the dynamic eggshell quality and located on chromosomes 1 and 4 from 57.3 to 71.4 Mb 104 . PIK3C2G possesses the C2 domain and acts as lipid binding motif, ITPR2 has been shown to be important to the process of eggshell calcification, and NCAPG gene was discovered to be associated with eggshell weight for young hens in a genome-wide association analysis 104 . Moreover, one association analysis determined that another lipid-related gene, low-density lipoprotein receptor-related protein 8 (LRP8), a new member of the egg shell matrix protein family, was significantly associated with eggshell traits 105 . Our MLM findings are in support of the aforementioned studies, in which 6 associations were identified on chromosome 1 between 60.20 and 67.94 Mb for EW2 and 5 associations were detected for EW1, ESS2, and ESW1 on chromosome 4 between 58.94 and 69.85 Mb. In our egg quality association study, we extended our investigation to internal aspects, such as yolk weight, diameter, color, and albumen weight. Understanding the genetic associations influencing internal egg quality is crucial due to the shift towards egg products in consumption. Through genetic marker analysis, we identified significant associations in internal egg characteristics, providing valuable insights for targeted breeding to enhance yolk weight and meet consumer preferences in the egg industry. Reference 91 identified the ataxia telangiectasia mutated (ATM) gene located on chromosome 1 (184.63 Mb) to be associated with the YW trait. In the present study, based on MLM analysis, we introduced 17 associations for YW1,2 on chromosome 1 between 67.94 and 172.99 Mb. Cathepsin D (CTSD) is another candidate gene that might affect egg characteristics 106 . Yolk formation involves cholesterol uptake and transport mediated by the very-low-density lipoprotein receptor on the membrane, and CTSD is the key enzyme regulating this process 107 . In this study, the detection of a wide range of associations (354) on all chromosomes (except for 22, 24, and 28) suggests that CTSD is a candidate gene associated with egg yolk quality traits. The albumen makes up approximately two-thirds of an egg’s weight. Numerous QTLs have been reported for egg white characteristics according to the chicken QTL database ( https://www.animalgenome.org/cgi-bin/QTLdb/GG/index ). Two promising genes, dopamine receptor D1 (DRD1) and msh homeobox 2 (MSX2) were found to be associated with albumen characteristics 108 . DRD1 encodes the D1 subtype of the dopamine receptor 109 . In birds, dopamine has been manifested to be involved in both stimulating and inhibiting prolactin (PRL) secretion, which has been illustrated to play an important role in the onset and maintenance of incubation behaviour 110 . Moreover, it has been shown to be associated with egg production traits 109 . MSX2 is a member of the msh homeobox family and is present in various embryonic tissues 108 . In the developing chick, MSX2 is expressed in the apical ectodermal ridge and the ectoderm of the genital tubercle; thus, playing a crucial role in the growth and patterning of the limb mesoderm 111 . Since MSX2 and DRD1 are involved in embryonic and ovarian development 108 , they could be treated as candidate genes associated with egg quality. Here, 23 associations were discovered for AW1,2 located on chromosomes 1–3, 8–9, 14, 21, 23, and Z in the GLM and MLM analyses. Yolk colours are important internal egg traits that are influenced by genetic and environmental factors, most often through the components which are included in the quail’s diet 112 . A total of 231 associations were found in the present study for SNPs that might influence yolk colours (L*, a*, and b* values) in both the GLM and MLM analyses. These associations were located on all chromosomes except for chromosomes 13, 22, 24, and 28. Genes underlying yolk colours have been previously reported in Galliformes. Reference 113 introduced three candidate genes, glucagon-like peptide 1 receptor (GLP1R), calcitonin receptor (CALCR), and tissue factor pathway inhibitor 2 (TFPI2), and suggested they are significantly associated with yolk colour. GLP1R is a member of the glucagon receptor family of G protein-coupled receptors found on beta cells in the pancreas and on neurons, and is involved in controlling blood sugar levels 114 . CALCR encodes a calcitonin receptor whose activity is mediated by G proteins which activate adenylyl-cyclase and is involved in calcium homeostasis, bone formation and metabolism, as well as lipid metabolism 115 . TFPI2 is a Kunitz-type serine protease inhibitor that inhibits matrix metalloproteinase activation and extracellular matrix degradation 116 . These genes are located on chromosome 5 within the 21.06–21.36 Mb genomic region 113 . In the present study, based on the MLM analysis, the nearest located SNP associations were in the 18.53–25.60 Mb range on chromosome 5 and responsible for YC-L*1 and YC-b*2 (Supplementary File S3 ). This similarity in positional candidate genes may indicate that our associations control yolk colour characteristics. In this study, we applied the Bonferroni-corrected GLM and MLM to investigate the genetic basis of targeted traits. Our R2 values conform to established statistical practices, considering multiple comparisons. The variability in R2 values can be attributed to trait complexity and genetic diversity. Comparative analysis with existing research supports the significance of our findings, contributing to the understanding of targeted traits. By the assessment of 44 phenotypic traits using GLM and MLM approaches, we showed the existence of SNP associations with all targeted traits, with the exception of EW1, AW1, and YC-b*2 in the GLM and ELA2 in the MLM analysis. The observed differences in GLM and MLM results may be influenced by factors like population structure, handling of missing data, genetic architecture, and assumptions made by each model. In agreement with previous studies, we hypothesize that multiple positional candidate genes are associated with body weight and egg-related traits. Future studies focusing on genotyping the SNPs within these candidate genes or linked markers would be of great use to improve the associations at the gene level. To the best of our knowledge, this is the first study to perform an association analysis of the production traits of Japanese quail. A total of 383 SNPs and 1151 associations were obtained following GLM, and 734 SNPs and 1442 associations were obtained following MLM analyses. The identified associations were significantly correlated with body weight and egg production. Moreover, 35 annotated genes were introduced as candidate genes for the targeted traits based on their nearest positions. Identifying associations and candidate genes could contribute to better understanding the genetic factors controlling body weight, egg production, and quality traits in quails particularly, and chickens generally, and may accelerate the genetic progress in breeding strategies. Data availability References Tsudzuki, M. Mutations of Japanese quail (Coturnix japonica) and recent advances of molecular genetics for this species. J. Poult. Sci. 45, 159–179 (2008). Ori, R. J., Esmailizadeh, A. K., Charati, H., Mohammadabadi, M. R. & Sohrabi, S. S. Identification of QTL for live weight and growth rate using DNA markers on chromosome 3 in an F2 population of Japanese quail. Mol. Biol. Rep. 41, 1049–1057 (2014). Charati, H., Koshkoiyeh, A. E., Ori, R. J., Moradian, H. & Mehrgardi, A. A. Detection of quantitative trait loci affecting carcass traits and internal organs on chromosome 3 in an F2 intercross of Japanese quail. Anim. Sci. Pap. Rep. 32, 369–383 (2014). Sun, Y. et al. Genome-wide linkage analysis and association study identifies loci for polydactyly in chickens. G3: Genes, Genomes, Genet. 4, 1167–1172 (2014). Sun, Y. et al. Genomewide association study of immune traits in chicken F2 resource population. J. Anim. Breed. Genet. 133, 197–206 (2016). Vollmar, S. et al. The gut microbial architecture of efficiency traits in the domestic poultry model species Japanese quail (coturnix japonica) assessed by mixed linear models. G3: Genes, Genomes, Genet. 10, 2553–2562 (2020). Kayang, B. B. et al. Integrated maps in quali (Coturnix japonica) confirm the high degree of synteny conservation with chicken (Gallus gallus) despite 35 million years of divergence. BMC Genom. 7, 101 (2006). Rights and permissions Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ .