Genome-Wide Association Study and Selective Sweep Analysis Reveal the Genetic Architecture of Body Weights in a Chicken F2 Resource Population.

Rapid growth is one of the most important economic traits in broiler breeding programs. Identifying markers and genes for growth traits may not only benefit marker-assisted selection (MAS)/genomic selection (GS) but also provide important information for understanding the genetic architecture of growth traits in broilers. In the present study, an F2 resource population derived from a cross between the broiler and Baier yellow chicken (a Chinese local breed) was used and body weights from 1 to 12 weeks of age [body weight (BW) 1-BW12)] were measured. A total of 519 F2 birds were genome re-sequenced, and a combination of genome-wide association study (GWAS) and selective sweep analysis was carried out to characterize the genetic architecture affecting chicken body weight comprehensively. As a result, 1,539 SNPs with significant effects on body weights at different weeks of age were identified using a genome-wide efficient mixed-model association (GEMMA) package. These SNPs were distributed on chromosomes 1 and 4. Besides, windows under selection identified for BW1-BW12 varied from 1,581 to 2,265. A total of 42 genes were also identified with significant effects on BW1-BW12 based on both GWAS and selective sweep analysis. Among these genes, diacylglycerol kinase eta (DGKH), deleted in lymphocytic leukemia (DLEU7), forkhead box O17 (FOXO1), karyopherin subunit alpha 3 (KPNA3), calcium binding protein 39 like (CAB39L), potassium voltage-gated channel interacting protein 4 (KCNIP4), and slit guidance ligand 2 (SLIT2) were considered as important genes for broiler growth based on their basic functions. The results of this study may supply important information for understanding the genetic architecture of growth traits in broilers.

Introduction

Chicken is one of the most economically important food production animals supplying meat and eggs to human beings. After intensive selection for the past 60 years, substantial advances have been made in improving body weight (BW) in modern commercial meat-type broilers (1). The selection for rapid growth will continue to be one of the most important economic traits in broiler breeding programs. Identifying genetic markers, especially single nucleotide polymorphisms (SNPs), and causal genes affecting BW can provide vital information for marker-assisted selection (MAS) and genomic selection (GS). Additionally, chicken is also considered an essential model for animal genomic studies (2). Therefore, the identification of genomic regions and potential candidate markers/genes can not only help understand the molecular mechanisms involved in the regulation of performance traits in the chicken, but also provide important information for the study in other species.

To date, 4,776 quantitative trait loci (QTLs) for chicken growth traits, including average daily gain and BW at different days of age, are hosted in the Chicken QTL database (release 45) (3). However, many of these QTLs, especially QTLs detected in previous studies, are coarsely mapped, which means that the confidence intervals of these QTLs are large and contain too many genes. The F2 design population is beneficial to QTL mapping of traits due to creation of larger genetic variation and trait segregation through the DNA recombination (4). In our previous study, several QTLs for growth and carcass traits were identified using microsatellite markers in the F2 resource population from a broiler × a Chinese local breed cross and these QTLs spanned large regions of the genome (5–9).

In the past decade, a genome-wide association study (GWAS) has been used to identify loci significantly associated with traits of interest of domestic animals using high-density SNP panels and genome resequencing technique. In chicken, a host of markers or genes important for growth, meat quality, fertility and so on, were identified in different populations using GWASs (10–16). A GWAS together with a selection signature analysis is widely used to identify SNPs and candidate genes associated with quantitative traits at the genome wide level.

In the present study, birds from an F2 resource population that was constructed by crossing broiler cocks derived from Arbor Acres with high abdominal fat content and Baier yellow chicken dams (a Chinese native breed) were genome re-sequenced, and a GWAS together with a selection signature analysis was carried out to comprehensively characterize the genetic architecture affecting BW in chickens. The results of this study can provide important information for understanding the genetic background of growth traits in chickens.

Materials and Methods

Animals and Phenotypic Measurements

All animal experiments were conducted according to the guidelines for the Care and Use of Experimental Animals established by the Ministry of Science and Technology of the People's Republic of China (approval number: 2006–398) and approved by the Laboratory Animal Management Committee of Northeast Agricultural University. This study used the F2 chicken resource population, which was described previously by Liu et al. (5). To establish this population, we crossed broiler sires derived from a high abdominal fat line divergently selected for abdominal fat with Baier yellow dams (a Chinese native breed). The F1 birds were intercrossed to produce an F2 population. All F2 birds had free access to feed and water. Commercial corn- and soybean-based diets that met all (17) requirements were provided in the study. From hatch to 3 weeks of age, the birds received a starter feed (3,000 kcal of ME/kg and 210 g/kg of CP) and from 3 to 12 weeks of age, the birds were fed a grower diet (3,100 kcal of ME/kg and 190 g/kg of CP) (18). The body weights (BWs) of a total of 519 F2 individuals (male and female birds) were measured at hatch and weekly up to 12 weeks of age. Then, quality control of BW was performed. Normality test was conducted to check the distribution of BW at every week using the Shapiro-Wilk test with JMP statistical software version 11.0 (SAS Institute Inc., Cary, NC, USA). If the traits were skewed from the normal test, outlier values were stepwise removed until the traits follow or roughly follow normal distribution. Then, phenotypic data was used for descriptive statistical analysis and GWAS.

DNA Library Preparation and Sequencing

Total genomic DNA was extracted from the chicken using the reagent test kit. For each bird, a single individual was used for genome sequencing on the Illumina HiSeq PE150 platform with an average depth of 3 ×. Library construction and sample indexing were done as described.

Population SNP Detection

Paired-end reads were mapped to the GCF_000002315.6_GRCg6a reference genome with Burrows-Wheeler Aligner (Version: 0.7.8) (19). The command line was “BWA mem -t 4 -k 32 –M.” After sorting, the “rmdup” command was used to remove potential PCR duplicates: only the pair with the highest mapping quality was retained, while multiple read pairs had identical external coordinates. After alignment, we performed SNP calling on a population scale with the package SAMtools (20). We then calculated genotype likelihoods from reads for each individual at each genomic location, and the allele frequencies in the sample were determined with a Bayesian approach. The “mpileup” command was used to identify SNPs with the parameters as “-q 1 -C 50 -S -D -m 2 -F 0.002 –u.” Then, to exclude SNP calling errors caused by incorrect mapping, only high-quality SNPs [coverage depth ≥2, root mean square (RMS) mapping quality ≥20, minor allele frequency (MAF) ≥0.05, and miss ≤0.3] were kept for subsequent analysis. After filtering from 15,868,916 raw SNPs, 10,889,955 SNPs remained. The missing genotypes of F2 individuals were imputed using 26 sequencing F0 individuals with 10-fold. Imputation was performed using BEAGLE 4.0 with default parameter settings (21). The genotypes for each individual were assumed to be unphased, and no relationships between individuals were used. Then, further quality control was conducted (filtered by MAF ≥0.05, missing rate ≤ 0.1, depth ≥2, and LD <0.6). Imputation accuracy (r) was calculated per SNP using the correlation between the observed and imputed genotypes. A total of 7,895,409 SNPs were left after the imputed 10,889,955 SNPs were filtered for the 519 individuals.

Functional Annotation of Genetic Variants

SNP annotation was performed according to the GCF_000002315.6_GRCg6a reference genome using the package ANNOVAR (Version: 2013-05-20) (22). Based on the genome annotation, SNPs were classified into exonic regions (overlapping with a coding exon), intronic regions (overlapping with an intron), splicing sites (within 2 bp of a splicing junction), upstream and downstream regions (within a 1 kb region upstream or downstream from the transcription start site), and intergenic regions. SNPs in coding exons were further grouped into synonymous SNPs (did not cause amino acid changes) or non-synonymous SNPs (caused amino acid changes). Candidate genes were screened in the 40 kb region upstream and downstream of each top SNP. Besides, mutations causing stop gain and stop loss were also classified into this group. Only the high-quality SNPs were annotated.

SNP-GWAS

In our association panel containing 519 samples, a total of 7,895,409 SNPs (left and filtered by MAF ≥0.05, missing rate ≤0.1, depth ≥2, and LD <0.6) were used in our GWAS for BW at 1–12 weeks of age (BW1–BW12). Association analysis was conducted using the genome-wide efficient mixed-model association (GEMMA) software package (23). For the mixed linear model analysis, we used the equation:

where y represents the phenotype; S is the incidence matrix of fixed effects and β is the vector of corresponding coefficients including the intercept; sex was included as a fixed effect to build up the S matrix. X represents the vector of SNP genotype and α is the corresponding effect of the marker; K is incidence matrix for μ, μ is the vector of random additive genetic effects following the multinormal distribution N (0, ), in which G is the genomic relationship matrix based on identity by state (IBS) [Genomic kinship f between individual i and j based on IBS is calculated using the following formula: . Where g(g) is the genotype of the i-th(j-th) bird at the k-th SNP. The frequency p is for the major allele and n is the number of SNPs (24)], and is the polygenetic additive variance. e represents random residual with a distribution of N (0, ). We performed principal component analysis (PCA) and tested the significance of top 10 PCAs using the EIGENSTRAT software, and the results showed there are no significant PCs in this population, indicating that there is no striking stratification. Therefore, PCs were not eventually included in the mixed model. The significant level was set as 0.05/N (P-value = 6.33 × 10−9) to control the genome-wide type 1 error rate and N is the number of informative SNPs.

Genome-Wide Selective Sweep Analysis

The BWs at each week in 519 birds were ranked. Selection signature analysis of BW was conducted between the two groups (15 birds per group) divided based on the highest and lowest BWs. Using VCFtools (25), we calculated the genome-wide distribution of fixation index (F) values and θπ ratios for the defined group pairs (high– and low– BW groups, 40-kb windows sliding in 10-kb steps) to characterize genome-wide selective sweeps related to the selection for growth rate. The θπ ratios were log2-transformed. Subsequently, the empirical percentiles of F and log2 (θπ ratio) in each window were estimated and ranked. The windows with the top 5% F and log2 (θπ ratio) values simultaneously were considered as candidate outliers under strong selective sweeps. All outlier windows were assigned to corresponding SNPs and genes. In other ways, the analysis of the allele frequency differences between the two groups was realized using VCFtools.

Results

Descriptive Statistics

The number of animals, means, and standard errors of BW1–BW12 are given in Table 1. Standard deviation (SD) and coefficient of variation (CV) of BWs vary from 10.09 to 392.62 grams and from 13.40 to 18.66%, respectively, indicating that there is large variability in BWs.

Table 1

Number of animals (N), mean (M), standard deviation (SD), minimum (MIN), maximum (MAX), and coefficient of variation (CV) of body weight at 1–12 weeks of age (BW1-BW12, in grams) of F2 chickens.

Trait	N	M	SD	MIN	MAX	CV (%)
BW1	492	75.30	10.09	50.5	102.4	13.40
BW2	490	169.12	22.17	116.1	235.9	13.11
BW3	499	308.70	43.66	170.0	425.0	14.14
BW4	499	479.87	70.48	250.0	660.0	14.69
BW5	489	643.71	94.06	365.0	885.0	14.61
BW6	499	844.44	128.03	480.0	1,160.0	15.16
BW7	504	1,091.51	172.27	660.0	1,535.0	15.78
BW8	500	1,290.74	213.07	790.0	1,845.0	16.51
BW9	493	1,542.36	264.30	945.0	2,185.0	17.14
BW10	504	1,730.23	301.01	1045.0	2,550.0	17.40
BW11	515	1,927.90	350.55	1170.0	2,880.0	18.18
BW12	519	2,104.00	392.62	1240.0	3,185.0	18.66

GWAS for BW at Different Weeks of Age

The GWAS for BW1–BW12 was carried out using the mixed-model statistical software package GEMMA (23) in 519 individuals of the F2 resource population (Figure 1; Supplementary Figure 1). A total of 1,539 SNPs with significant effects (P < 6.33 × 10−9) on BW1–BW12 were detected (Table 2; Supplementary Table 1). These SNPs with significant effects on BW1–BW12 were distributed on chromosomes 1 and 4. Two lead SNPs responsible for BW12 were detected at the 171,411,019bp (rs316877904) on chromosome 1 (P = 8.49e-19), and at the 74,526,009bp (rs13774694) on chromosome 4 (P = 3.54e-10), respectively. Linkage disequilibrium (LD) analyses showed that these two lead SNPs were in high LD with some near SNPs (Supplementary Figure 2). Additionally, an interesting phenomenon is that the number of significant SNPs and genes within the region identified on chromosome 1 is consecutively increasing accompanied with weeks of age (from 5 to 12 weeks of age). The genes in 40-kb regions of the SNPs with significant effects on BW1–BW12 were extracted, and 265 genes were detected (Table 2).

Figure 1

Results of genome-wide association studies for body weight at 12 weeks of age (BW12) using the GEMMA package. The results are presented as the Manhattan plot in the left panel and the Q–Q plot in the right panel. The solid line indicates the threshold to control the genome-wide type I error of 5% (P < 6.33 × 10−9). The Q-Q plot was used to estimate the difference between observed and expected chi-square statistic values of quantitative traits, indicating that the potential candidate loci related to the traits were not caused by population stratification and the statistical model was reasonable.

Table 2

Number of SNPs with significant effects on body weight at different weeks of age and number of selected windows by selective sweep analysis.

Traits	Number of significant SNPs	Number of genes detected by GWAS	Number of selected windows	Number of genes detected by selective sweep analysis
BW1	0	0	1,581	507
BW2	0	0	1,924	468
BW3	0	0	1,837	499
BW4	0	0	1,842	479
BW5	12	9	1,867	502
BW6	20	15	1,713	429
BW7	79	23	1,973	459
BW8	292	49	1,794	498
BW9	122	30	1,945	523
BW10	203	34	2,020	501
BW11	384	50	2,137	476
BW12	427	55	2,265	538

Selective Sweep Detected by FST Combined With Pi Methods

We ranked 519 birds according to BW values at every week. The two groups (15 birds per group) were divided in the light of the highest and lowest BWs (Supplementary Table 2). The genome-wide selection signatures were detected using both fixation index (F) values and θπ ratios. Windows with the top 5% F and outliers of log2 (θπ ratio) values simultaneously were considered as essential regions under strong selective sweeps. A total of 1,581–2,265 windows under selection were identified for BW1–BW12 (Table 2; Figure 2, Supplementary Figure 3). After deleting the overlaps, 8,554 selected windows were left (Supplementary Table 3). These selection signatures were distributed on nearly all chromosomes with several peaks on chromosomes 1, 4, 5, and Z. All windows under selective sweeps were assigned to corresponding SNPs and genes, and about 429–538 annotated genes were identified. After deleting the overlaps, 2,812 genes were left.

Figure 2

Selective sweep analysis for body weight at 12 weeks of age (BW12) detected by FST and Pi methods. Blue and green colors indicate windows with the top 5% FST and log2 (θπ ratio) values simultaneously, which were considered as the selective sweeps.

Candidate Genes for BW at Different Weeks of Age

After compared the GWAS and selection signature gene lists, we found 42 genes overlap (Table 3). Some genes, including diacylglycerol kinase eta (DGKH), deleted in lymphocytic leukemia (DLEU7), forkhead box O17 (FOXO1), karyopherin subunit alpha 3 (KPNA3), calcium binding protein 39 like (CAB39L), potassium voltage-gated channel interacting protein 4 (KCNIP4), and slit guidance ligand 2 (SLIT2), were reportedly thought to be necessary for chicken growth based on their basic function study.

Table 3

Overlap genes detected by GWAS and selective sweep analysis.

No	Gene	Description	Location
1	DGKH	diacylglycerol kinase eta	chr1:167536310-167702128
2	LOC112531567		chr1:170650308-170655717
3	LOC112531569		chr1:171083270-171100945
4	DLEU7	deleted in lymphocytic leukemia, 7	chr1:171143861-171152064
5	SERPINE3	serpin family E member 3	chr1:171398844-171416703
6	WDFY2	WD repeat and FYVE domain containing 2	chr1:171476437-171542083
7	FOXO1	forkhead box O1	chr1:171900263-171963540
8	GTF2F2	general transcription factor IIF subunit 2	chr1:168988717-169082382
9	GPALPP1	GPALPP motifs containing 1	chr1:168950859-168965646
10	FNDC3A	fibronectin type III domain containing 3A	chr1:170318524-170431952
11	KPNA3	karyopherin subunit alpha 3	chr1:170597160-170650244
12	FAM124A	family with sequence similarity 124 member A	chr1:171336721-171377902
13	SETDB2	SET domain bifurcated 2	chr1:170526801-170564892
14	LOC112531568		chr1:171008615-171039942
15	MRPS31	mitochondrial ribosomal protein S31	chr1:171849779-171873661
16	LHFP	lipoma HMGIC fusion partner-like 1	chr1:172287762-172427229
17	NHLRC3	NHL repeat containing 3	chr1:172554893-172565969
18	LOC100859822		chr1:169095074-169099300
19	SLC25A30	solute carrier family 25 member 30	chr1:169101373-169110759
20	COG3	component of oligomeric golgi complex 3	chr1:169110893-169143086
21	CDADC1	cytidine and dCMP deaminase domain containing 1	chr1:170447991-170463866
22	CAB39L	calcium binding protein 39 like	chr1:170465092-170526727
23	ARL11	ADP ribosylation factor like GTPase 11	chr1:170586604-170597668
24	LOC107051704		chr1:170658231-170671446
25	RNASEH2B	ribonuclease H2 subunit B	chr1:171220939-171282281
26	TPTE2	transmembrane phosphatase with tensin homology	chr1:171808785-171827969
27	LOC101750153		chr1:172061063-172091470
28	COG6	component of oligomeric golgi complex 6	chr1:172218064-172269404
29	LOC107052027		chr1:172532497-172555759
30	PROSER1	proline and serine rich 1	chr1:172566162-172586926
31	LOC107052035		chr1: 172584113-172597742
32	C4A	complement C4A (Rodgers blood group)	chr1:172599919-172652724
33	FREM2	FRAS1 related extracellular matrix protein 2	chr1:172657847-172776158
34	LOC112532426		chr1:74546421-74555772
35	RUBCNL	rubicon like autophagy enhance	chr1:169433550-169453972
36	LRCH1	leucine rich repeats and calponin homology domain containing 1	chr1:169520469-169644616
37	DCLK1	doublecortin like kinase 1	chr1:174079577-174317535
38	SMIM20	small integral membrane protein 20	chr4:73372750-73376776
39	SLC34A2	solute carrier family 34 member 2	chr4:73421043-73440963
40	ANAPC4	anaphase promoting complex subunit 4	chr4:73476279-73507544
41	SLIT2	slit guidance ligand 2	chr4:74981753-75225786
42	KCNIP4	potassium voltage-gated channel interacting protein 4	chr4:74568444-74948433

Discussion

Deciphering of the genetic underpinning of chicken growth traits is conducive to further genetic improvement in breeding program. Owing to the genetic hitch-hiking effect and relatively large QTL confidence intervals of F2 population (26), we leveraged an integrative strategy that couples GWAS with selection signatures analysis to dissect the genetic determinants of chicken growth traits at the genome-wide level.

A GWAS has been widely used to identify SNPs and candidate genes associated with important quantitative traits at the genome-wide level. Some important regions associated with production, reproduction, and disease resistance traits in chickens have been identified using GWASs (11, 27–32). Growth trait, especially body weight, is one of the most important economic traits in the poultry industry. Therefore, in this study, we carried out a GWAS for BW from 1 to 12 weeks of age using an F2 resource population established by crossing broiler sires with Baier yellow dams. A total of 1,539 SNPs with significant effects on BW1-BW12 were distributed on chromosomes 1 and 4, indicating that BWs are complex traits controlled by multiple genetic determinants. Intriguingly, it could be observed from Supplementary Figure 1 and Table 2 that more significant SNPs were detected over increased ages and BW (from BW5 to BW12). We speculate that there are two possibilities: one is that more genes are likely to be involved in chicken growth and development at the late developmental stages, another is that the involved genes probably play growingly important roles over increased ages. Additionally, according to the numbers of significant SNPs in the region identified on the chromosome 1 for BW5–BW12 (12, 20, 79, 292, 122, 202, 374 and 414, respectively), we found that the top 10 significant SNPs did not change over ages. For instance, 3 SNPs (Chromosome1: 171411019, 171710979, 171925771) are consecutively significant from 5 to 12 weeks; 2 SNPs (Chromosome1: 170926098, 170930631) from 6 to 12 weeks; 5 SNPs (Chromosome1: 170713988, 170714148, 170715097, 170926202, 170926223) from 7 to 12 weeks. These findings reveal that there are shared genetic determinants for BW5–BW12, and these SNP are pleiotropic variants.

To date, an array of QTLs for chicken BW have been identified, and some of them were consistent with the results of the present study. Xu et al. (27) reported that chromosomes 1 and 4 were the two critical chromosomes influencing growth traits, particularly BW in chickens. In the present study, chromosomes 1 and 4 were detected by both GWAS and selective sweep analysis; they were found to harbor important genes for chicken BW. Podisi et al. (33) also reported two significant QTLs for BW at 12 weeks of age on chromosome 1 in broiler cross-bred female chickens. Mebratie et al. (34) carried out a GWAS for BW and found that SNPs with significant effects on BW were located on chicken chromosomes 1, 6, 8, 12, 14, 23, and (35) also identified some SNPs with significant effects on chicken BW located on chromosomes 1, 2, 3, 4, 5, 6, 7, 8, 10, 14, and 21. Chromosome 1 was also identified as harboring important genes for chicken BW in our previous report using the same population by the marker-QTL linkage analysis (6, 8).

As an importantly economic trait, body weight in chickens has undergone long-term artificial selection in the past decades, which is expected to left selective signatures on chicken genomics. The detection of selection signatures can expedite the identification of genes responsible for important economic traits and better understanding the biological mechanisms affected by strong ongoing natural or artificial selection in livestock populations (36, 37). Accordingly, in this study selective signature analysis was utilized to confirm overlapping genomic regions detected by GWAS to screen important candidate genes for chicken BW.

Forty-two annotated genes of chicken were identified by both GWAS and selective sweep analysis. The basic functions of these 42 genes were extracted from the previous reports. Some genes, including DGKH, DLEU7, FOXO1, KPNA3, CAB39L, KCNIP4, SLIT2, which were found to be associated with growth traits in farm animals, were considered as important candidate genes for growth traits in broilers. DGKH was identified as a candidate gene affecting divergent growth in cattle (38), and this gene could regulate the growth of cattle by regulating the secretion of growth-related hormones (39). KPNA3 was found to be associated with chicken growth traits in a previous GWAS (10, 40). Zhang et al. (41) identified the CAB39L could be a candidate gene for growth and carcass traits by GWAS and pathway enrichment analysis in a Gushi-Anka F2 chicken population. The SNPs of DLEU7 gene were associated with growth traits in Jinghai yellow chickens (40). The replication of DLEU7 was associated with height in African-derived populations (42). FOXO1 could influence food intake and then regulate growth (43, 44). KCNIP4 was identified to be associated with BW of chicken using a GWAS (28). SLIT2 was also found to be associated with BW at 35 and 41 days of age in chickens (45). Apart from above-mentioned 7 genes reportedly associated with growth traits in farm animals, others 35 genes can be considered as novel candidate genes for chicken growth and development.

Observations at multiple time points for the same individual are called longitudinal traits, which can better describe the growth and production of farm animals than single data records (46). Chicken BWs at different weeks of age are classic longitudinal traits. In the present study GWAS was independently performed for every time point to dissect the genetic basis of BW. A better alternative strategy is to fit the growth curve and then conduct the association analysis using the fitted parameters, which could better mirror growth trajectory and provide novel insight into genetic underpinning of BW in the chicken.

Conclusions

In summary, in this study, both GWAS and selective sweep analysis were carried out to identify important SNPs and genes for chicken BW. Finally, 42 genes were detected, and some genes, including DGKH, DLEU7, FOXO1, KPNA3, CAB39L, KCNIP4, SLIT2 were identified as important candidate genes for rapid growth in chickens.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: [NCBI SRA AND PRJNA861112].

Ethics Statement

The animal study was reviewed and approved by the Laboratory Animal Management Committee of Northeast Agricultural University. Written informed consent was obtained from the owners for the participation of their animals in this study.

Author Contributions

SW and YW carried out the experiments, performed the statistical analyses, and prepared the manuscript. YL contributed to the design of the experiments and writing the manuscript. FX, HuG, HaG, and NW contributed to writing the manuscript. HZ and HL conceived and designed the study, participated in data interpretation, and contributed to writing the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the National Key R&D Program of China (Grant No. 2021YFD1300100), the Joint Guidance Project of Heilongjiang Natural Science Foundation (No. LH2021C036), the National Natural Science Foundation of China (Nos. 31572394 and 31972549), China Agriculture Research System of MOF and MARA (No. CARS-41), and Project of the Ministry of Agriculture and Rural Affairs in China (No. 19190526).

Conflict of Interest

FX, HsG, and HhG are/were employed by Fujian Sunnzer Biotechnology Development Co., Ltd, China. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.