GWAS and Linkage Disequilibrium (LD) Analysis

GWAS was performed with EMMAX program (Kang et al. 2010). Briefly, we first used PLINK to get the genotype file, and the ETN phenotype data were collected from two locations, Jiangxi and Hainan. The genotype and phenotype file, as well as the Balding-Nichols kinship matrix were used to fit the mixed linear model in EMMAX. The Manhattan and quantile-quantile (Q-Q) plots were generated using the R package “qqman” (Turner. 2018). The pairwise linkage disequilibrium (LD) (r²) was calculated using PopLDdecay software (Zhang et al. 2019). The Pairwise r² was calculated for the SNPs in 500 kb and averaged in 1 kb across the whole genome. The LD decay rate was measured as the chromosomal distance at which the average pairwise correlation coefficient (r²) decreased to half of its maximum value (Huang et al. 2010). The regional Manhattan plot and LD heatmap were obtained using LDBlockShow software (Dong et al. 2020).