MATERIALS AND METHODS

Ethical approval

All procedures involving animals were reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of Mahasarakham University, Mahasarakham, Thailand (Approval No. IACUC-MSU-026-016/2025). The study was conducted in accordance with institutional guidelines for the ethical care and use of animals in research and complied with the principles of the 3Rs (Replacement, Reduction, and Refinement). All procedures were designed to minimize animal discomfort and distress, and sample collection and animal handling were performed by trained personnel in accordance with accepted veterinary and animal welfare practices.

Study period and location

The study was conducted from January 2018 to November 2020 using Pradu Hang Dum Thai native chickens maintained at the Chiang Mai Livestock Research and Breeding Center, Chiang Mai, Thailand. The production system represented a tropical semi-intensive environment typical of sustainability-oriented native chicken breeding programs. Ambient temperatures ranged from 21.3°C to 32.5°C, with an average relative humidity of approximately 52%.

Study design

A longitudinal, population-based genetic study was conducted on 563 female Pradu Hang Dum chickens from two consecutive generations. Phenotypic, pedigree, and genotypic information was integrated to investigate relationships among growth performance, sexual maturity, egg production, and egg-laying persistency. Multivariate analyses, genetic parameter estimation, and candidate–gene association analyses were performed to characterize growth–maturity–egg production trade-offs and to evaluate marker effects on EBVs.

Animals and management

Animals used in this study were obtained from the Chiang Mai Livestock Research and Breeding Center (Chiang Mai, Thailand). A total of 563 female Pradu Hang Dum chickens from two generations (2018–2020) were included. This population structure under tropical semi-intensive management was considered representative of sustainability-oriented native chicken breeding systems.

The breeding population was maintained as a closed nucleus population with within-breed selection, in which approximately 30 cocks and 150 hens were selected in each generation (effective population size ≈100). Each cock was mated with five hens through artificial insemination, and mating plans were designed to minimize inbreeding. Pedigree records comprised 5,676 individuals. Approximately 2,500 chicks were produced per generation, from which approximately 40 cocks and 270 hens were retained for performance testing. Hens exhibiting poor laying performance or failure to lay were removed before selecting the top 30 cocks and 150 hens for the subsequent generation.

Birds were individually identified using wing bands and reared in open-sided houses with natural ventilation. After a 21-day incubation period, chicks were stocked at a density of 8 birds/m² until 18 weeks of age, then transferred to individual battery cages.

Birds were fed according to age. From hatch to 4 weeks, chicks received a starter diet containing 21% crude protein (CP) and 3,000 kcal/kg metabolizable energy (ME). From 4 to 12 weeks, a grower diet containing 19% CP and 2,900 kcal/kg ME was provided. Feed and water were supplied ad libitum throughout the experimental period.

Continuous lighting (24 h light:0 h dark) using 100-W heating lamps was applied from hatch to 4 weeks, followed by natural daylight from 4 to 12 weeks of age. After 18 weeks of age, birds were maintained under ambient tropical conditions. At 20 weeks of age, hens were transferred to individual cages (20 × 45 × 40 cm; floor area ≈0.09 m²) and remained individually housed throughout the laying period. All birds received the same commercial layer diet (19% CP and 2,900 kcal ME/kg) at 110 g/day with free access to water. A consistent 12 h light/day schedule was maintained during the laying period.

All birds were managed in accordance with Good Agricultural Practices, and vaccination schedules were implemented by the farm veterinarian in accordance with Thai Agricultural Standards. Under semi-intensive management, Thai native chickens typically attain market weights of 1.2-1.4 kg at 12-14 weeks of age, with feed conversion ratios of 3.0-3.5, survival rates exceeding 90%, and carcass yields ranging from 68% to 72%.

Data collection and measurements

Growth performance traits: Growth performance was evaluated using standardized procedures. Individual body weights were recorded using calibrated scales at hatch (BW0) and at 4, 8, 12, and 16 weeks of age (BW4, BW8, BW12, and BW16, respectively). All measurements were recorded individually to ensure accuracy and reliability.

Egg production traits: Seven egg production-related traits were recorded, including body weight at first egg (AFE_WT), AFE, egg weight at first egg (AFE_EW), egg weight at 270 days (EW270), egg weight at 360 days (EW360), cumulative egg number at 270 days (EN270), and cumulative egg number at 360 days (EN360). Measurements were obtained over a complete 365-day laying cycle. Unequal sample sizes across traits resulted from incomplete records and were addressed using mixed-model analyses.

Egg-laying persistency traits: Because repeated longitudinal egg production records over the entire laying cycle were unavailable, simplified persistency indicators based on late-stage egg production were used instead of model-based persistency measures.

The persistency ratio (PR) was defined as

\documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} \[ PR=\frac{{R}_{late}}{{R}_{early}} \] \end{document}

Where

\documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} ({R}{early}=\frac{E{N}{early}}{\left(270,AFE\right)}) \end{document}and \documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} ({R}{late}=\frac{E{N}{late}}{90}) \end{document}

Persistency late share (PLS) was calculated as

\documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} \[ PLS=\frac{E{N}_{late}}{E{N}_{360}}=\frac{E{N}_{late}}{E{N}_{early}+E{N}_{late}} \] \end{document}

where \documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} (E{N}{early}) \end{document}represents the number of eggs produced from AFE to 270 days of age and \documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} (E{N}{late}) \end{document}represents the number of eggs produced between 271 and 360 days of age. Therefore, the early laying period was standardized according to individual AFE. These simplified indices were considered appropriate for indigenous chickens, which characteristically exhibit clutch-based laying patterns with prolonged but gradually declining production under tropical conditions.

DNA extraction and candidate gene genotyping

Blood samples (1 mL) were collected from the wing vein into 1.5-mL microcentrifuge tubes containing 0.5 M EDTA. Genomic DNA was extracted from whole blood using the guanidine hydrochloride method [30].

Polymorphisms in the candidate genes DRD2, VIP, MTNR1C, and NPY were analyzed using polymerase chain reaction–restriction fragment length polymorphism (PCR-RFLP). Gene-specific primers were used to amplify DRD2 (chromosome 24; GenBank accession number 428252), VIP (chromosome 3; GenBank accession number 396323), MTNR1C (chromosome 4; GenBank accession number JQ249896), and NPY (chromosome 2; GenBank accession number M87298).

PCR amplification was performed using annealing temperatures of 60°C for DRD2 and NPY and 58°C for VIP and MTNR1C. Amplified fragments were digested with _Bse_GI, _Vsp_I, _Mbo_I, and _Dra_I, respectively, and separated on 2.5% agarose gels. Genotypes were assigned according to restriction fragment patterns [22, 23] (Figure 1). Samples with ambiguous banding patterns were excluded. To confirm genotyping accuracy, a subset of samples was randomly re-genotyped. All polymorphisms investigated in this study have previously been reported as synonymous or non-coding variants [22, 23].

Figure 1: Genotyping patterns of candidate gene polymorphisms associated with growth and reproductive traits in Thai native chickens using polymerase chain reaction–restriction fragment length polymorphism analysis. (A) DRD2 digested with BseGI; M = 100 bp DNA ladder; TT genotype (196 and 52 bp), TC genotype (248, 196, and 52 bp), and CC genotype (248 bp). (B) VIP digested with VspI; II genotype (306 bp), ID genotype (306, 154, and 152 bp), and DD genotype (154 and 152 bp). (C) MTNR1C digested with MboI; AA genotype (372 bp), AG genotype (372 and 333 bp), and GG genotype (333 bp). (D) NPY digested with DraI; LL genotype (240 bp), MM genotype (240, 161, and 79 bp), and SS genotype (161 and 79 bp).

Statistical analysis

PCA: Trait clustering was evaluated using a PCA-derived path diagram to visualize relationships among principal components and phenotypic performance traits at the individual bird level. The analysis was performed using JASP statistical software version 0.95.4 (JASP Team). All phenotypic variables were standardized before analysis, and rotated component structures were obtained using the oblique via cluster method to improve interpretability. Complete rotated component loadings and communalities are provided in Supplementary Table-S1. This approach facilitated the identification of biological and physiological relationships among trait groups.

Subsequently, PCA was performed using combined genotype data from DRD2, VIP, MTNR1C, and NPY to characterize the multivariate genetic structure underlying variation in individual performance. Unlike previous applications of PCA to egg production stages in commercial layers or egg quality traits, the present study incorporated growth, sexual maturity, and two distinct persistency indicators (PR and PLS) to resolve biological trade-offs in a slow-growing indigenous population.

Genotype and allele frequencies, polymorphism information content, and expected heterozygosity were calculated as described by Falconer and Mackay [31]. Conformity with Hardy–Weinberg equilibrium was evaluated using the chi-square (χ²) test.

Genetic parameter estimation: Contemporary group, defined by hatch and generation, was included as a fixed effect for growth traits and AFE in Model 1. For egg production traits, including AFE_WT, AFE_EW, egg weights at 270 and 360 days (EW270 and EW360), egg numbers at 270 and 360 days (EN270 and EN360), and persistency traits (PR and PLS), AFE was included as a covariate in Model 2.

Variance components and heritability estimates were obtained using single-trait animal models, whereas genetic correlations were estimated using bivariate linear animal models. Variance components were estimated using the average information restricted maximum likelihood algorithm implemented in the BLUPF90 program [32]. Heritability estimates were obtained using single-trait animal models, and genetic correlations were estimated using bivariate analyses. Convergence was assumed when changes in log-likelihood values between successive iterations were <10⁻¹².

EBVs derived from single-trait models were subsequently used in candidate gene association analyses.

The bivariate animal model is expressed as:

\documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} \[ \begin{aligned} {y}_{1}={X}_{1}{\beta }_{1}+{Z}_{1}{a}_{1}+{\varepsilon }_{1} \\ {y}_{2}={X}_{2}{\beta }_{2}+{Z}_{2}{a}_{2}+{\varepsilon }_{2} \end{aligned} \] \end{document}

where y₁ and y₂ are vectors of observations for traits 1 and 2, respectively; β represents the vector of fixed effects; a denotes additive genetic effects; and ε represents residual effects. X and Z are incidence matrices relating observations to fixed and random effects, respectively. Additive genetic effects were assumed to follow:

\documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} \[ a∼N(0,A{\sigma }_{a}^{2}) \] \end{document}

where A is the additive relationship matrix.

Heritability (h²) was calculated as:

\documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} \[ {h}^{2}=\frac{{\sigma }_{a}^{2}}{{\sigma }_{a}^{2}+{\sigma }_{e}^{2}} \] \end{document}

where σ²a is additive genetic variance and σ²e is residual variance.

Genetic (rg) and residual (re) correlations between traits were estimated as:

\documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} \[ \begin{aligned} {r}_{g}=\frac{{\sigma }_{a12}}{\sqrt{{\sigma }_{a1}^{2}{\sigma }_{a2}^{2}}} \\ {r}_{e}=\frac{{\sigma }_{e12}}{\sqrt{{\sigma }_{e1}^{2}{\sigma }_{e2}^{2}}} \end{aligned} \] \end{document}

where σa12 and σe12 represent additive genetic and residual covariances between traits, respectively.

Candidate gene association analysis: Associations between candidate gene polymorphisms and EBVs were evaluated using a linear mixed model. Genotype was fitted as a fixed effect for growth traits (BW0, BW4, BW8, BW12, and BW16) and reproductive traits (AFE, AFE_WT, AFE_EW, EW270, EW360, EN270, and EN360). Analyses were performed using SAS software version 9.4 [33].

The model was expressed as:

\documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} \[ {y}_{ijk}=\mu +{Gene}_{j}+{\varepsilon }_{ijk} \] \end{document}

where \documentclass{article} \usepackage{amsmath} \usepackage{amssymb} \begin{document} ({y}_{ijk}) \end{document}represents the EBV for trait i, μ is the overall mean, Genej is the fixed effect of genotype j, and εijk is the random residual error.

No additional covariates were included because EBVs had already been adjusted for pedigree structure and fixed environmental effects through animal-model analyses. Consequently, using EBVs rather than raw phenotypic values enabled candidate gene associations to be evaluated while accounting for pedigree relationships and environmental influences.