Genome-Wide Association Study on Root System Architecture and Identification of Candidate Genes in Wheat (Triticum aestivum L.)

The root tissues play important roles in water and nutrient acquisition, environmental adaptation, and plant development. In this study, a diversity panel of 388 wheat accessions was collected to investigate nine root system architecture (RSA) traits at the three-leaf stage under two growing environments: outdoor pot culture (OPC) and indoor pot culture (IPC). Phenotypic analysis revealed that root development was faster under OPC than that under IPC and a significant correlation was observed between the nine RSA traits. The 660K single-nucleotide polymorphism (SNP) chip was used for a genome-wide association study (GWAS). Significant SNPs with a threshold of −log10 (p-value) ≥ 4 were considered. Thus, 36 quantitative trait loci (QTLs), including 13 QTL clusters that were associated with more than one trait, were detected, and 31 QTLs were first identified. The QTL clusters on chromosomes 3D and 5B were associated with four and five RSA traits, respectively. Two candidate genes, TraesCS2A01G516200 and TraesCS7B01G036900, were found to be associated with more than one RSA trait using haplotype analysis, and preferentially expressed in the root tissues. These favourable alleles for RSA traits identified in this study may be useful to optimise the root system in wheat.


Introduction
Bread wheat (Triticum aestivum L.) is one of the most important staple crops worldwide. The cultivated area of wheat is approximately 217 million hectares, annually producing approximately 765 million tonnes [1]. However, the demand for wheat continues to increase owing to a rapid increase in population. According to the speculation by the Food and Agriculture Organisation of the United Nations, approximately 840 million tonnes of wheat will be required annually by 2050 [2]. However, wheat production is facing serious challenges with the increasing global climate change in recent years [3,4]. Wheat breeding for high yield and stress tolerance will remain the main strategy for a long time. Roots, the underground parts of plants, absorb nearly all of the nutrients required for plant development and are primarily exposed to osmotic stress [5]. Optimal RSA improves the wheat yield and stress tolerance [6,7]. Therefore, it is important to study the genetic mechanisms of crop RSA traits to optimise the root system for further breeding.

Population Structure and Linkage Disequilibrium (LD) Decay Analysis
The population structure was determined using STRUCTURE V2.3.4 according to the ∆K method of Bayesian clustering. We determined the slope breakpoint at K = 8. Thus, the panel used in this study was divided into eight subpopulations (Sp1-Sp8, Figure 1a, Table S1). Sp1 comprised 11 exotic varieties, 6 landraces, and 49 local cultivars. Sp2 consisted of 7 exotic varieties, 4 landraces, and 29 local cultivars. Sp3 comprised 5 exotic varieties, 1 landrace, and 58 local cultivars. Sp4 contained 13 exotic varieties, 15 landraces, and 19 local cultivars. Sp5 was composed of 10 exotic varieties, 3 landraces, and 28 local cultivars. Sp7 consisted of 1 exotic variety, 1 landrace, and 28 local cultivars. Sp8 contained 19 exotic varieties, 7 landraces, and 40 local cultivars (Table 1, Figure 1b). quency correlation (r 2 ). A plot of the LD estimates (r 2 ) as a function of physical distance (Mb) indicated clear LD decay with increasing physical distance (Figure 1c). A comparison of LD between the sub-genomes and chromosomes showed variances in the LD decay. Overall, the average LD decay distance in the whole genome was approximately 3.6 Mb. LD decayed faster in the D genome (1.6 Mb) than in the A (3 Mb) and B (5.6 Mb) genomes ( Figure 1c). The faster LD decay in the D genome is compatible with the evolutionary history of wheat.  The different-coloured curves represent LD decay fits for sub-genomes A (purple), B (blue), D (green), and the whole genome (red). The thick horizontal black line represents the population-specific critical r 2 value (0.1) above which LD may be due to linkage.  Sp1  66  11  6  49  Sp2  41  7  5  29  Sp3  64  5  1  58  Sp4  47  13  15  19  Sp5  41  10  3  28  Sp6  33  0  0  33  Sp7  30  1  1  28  Sp8  66  19  7  40 We estimated the degree of LD using 161,801, 184,658, and 58,479 SNP markers in the A, B, and D sub-genomes, respectively. LD was analysed based on 717,701,068 pairwise comparisons of 404,938 SNPs. Pair-wise LD was estimated using the squared allelefrequency correlation (r 2 ). A plot of the LD estimates (r 2 ) as a function of physical distance (Mb) indicated clear LD decay with increasing physical distance (Figure 1c). A comparison of LD between the sub-genomes and chromosomes showed variances in the LD decay. Overall, the average LD decay distance in the whole genome was approximately 3.6 Mb. LD decayed faster in the D genome (1.6 Mb) than in the A (3 Mb) and B (5.6 Mb) genomes ( Figure 1c). The faster LD decay in the D genome is compatible with the evolutionary history of wheat.

Phenotypic Variation of the RSA Traits
A total of 4656 (388 wheat accessions × 6 replications × 2 environments) seedling roots were collected and observed for the nine RSA traits namely, total root length (TRL), total root surface area (TRA), average root diameter (ARD), total root volume (TRV), number of root tips (NRT), root length (RL), root surface area (RA), root volume (RV), and root tips (RT) ( Table S1). The frequency distribution of the investigated traits was normal, as indicated by the kurtosis (bk) and skewness (bs) (Figure 2). This means that the traits were quantitative and suitable for GWAS. Phenotypic variation analysis showed that the nine RSA traits exhibited greater variation under OPC than under IPC, indicated by the coefficient of variation (Table 2). In addition, the nine RSA traits were significantly higher under OPC than under IPC (p ≤ 0.001). In other words, the root development was faster under OPC than that under IPC. The broad-sense heritability (H 2 ) of the nine RSA traits ranged from 66.78% to 77.83%, indicating that the variations observed in these RSA traits were mainly regulated by genotype.

Correlations between RSA Traits
The relationships of the nine RSA traits between the two environments were a lysed. The results showed that each RSA trait showed a significant correlation (p ≤ 0.0

Correlations between RSA Traits
The relationships of the nine RSA traits between the two environments were analysed. The results showed that each RSA trait showed a significant correlation (p ≤ 0.001) between the two environments. This means that these traits were mainly regulated by genotype. The relationships between the nine traits were calculated ( Figure 3). The results showed that significant positive correlations were observed between eight RSA traits, except for ARD. ARD showed a significant negative correlation with NRT, RL, RA, RV, and RT in the two environments. The correlation coefficients between RL, RA, and RV were higher than 0.94 in the two environments, indicating that these three traits were strongly correlated. 43 6 of 20

GWAS of RSA Traits
After filtering low-quality SNPs (minor allele frequency < 0.05 and missing data > 0.1), a total of 411,605 polymorphic SNP markers with effective chromosome information

GWAS of RSA Traits
After filtering low-quality SNPs (minor allele frequency < 0.05 and missing data > 0.1), a total of 411,605 polymorphic SNP markers with effective chromosome information were available for GWAS using a univariate mixed linear model of GEMMA (Table S2). Of these, 404,938 SNPs were mapped to 21 wheat chromosomes, with 161,801, 184,658, and 58,479 SNPs in the A, B, and D sub-genomes, respectively. The marker density varied among different chromosomes ( Figure S1, Table 3). We discovered that the minimum marker density was 8.1 markers per Mb on chromosome 4D, whereas the maximum was 55.09 markers per Mb on chromosome 3B. Compared with the A and B sub-genomes, the D sub-genome had fewer SNP markers and effective markers, especially chromosome 4D, which contained 1702 SNPs. The markers of known chromosome positions were used to analyse the genetic diversity. We found that there was little difference in the genetic diversity among the three sub-genomes. The mean values of the expected heterozygosity (He) and polymorphism information content (PIC) were 0.66-0.74 and 0.27-0.30 among the three sub-genomes, respectively.
We performed a GWAS using the best linear unbiased predictions (BLUPs) of the nine RSA traits. Significant SNPs for each trait were selected with a threshold of -log 10 (p-value) ≥ 3.5 (Table S3, Figure S2). After analysing the data, 133 QTLs evenly distributed on 21 chromosomes were identified for the nine RSA traits, among which 38 were associated with multiple RSA traits. Eight QTLs were overlapped with previously identified QTLs (Table S3).
We focused on significant SNPs with a threshold of -log 10 (p-value) ≥ 4. For convenience, we selected tagged SNPs for each QTL exhibiting the strongest association with RSA traits (Table 4), yielding 36 QTLs on 17 chromosomes (Figure 4), of which 13 were associated with more than one RSA trait. Among these QTLs, five QTLs, mapped on chromosomes 5B, 6B, 6D, 7A, and 7B, were co-localised with those identified in previous studies ( Table 4), indicating that our results were reliable. Table 3. Summary of the genetic diversity in three sub-genomes and chromosomes of this wheat panel and evaluation of the effective number of independent SNPs, including the suggested pvalue thresholds.

Chromosome
No  QTLs (Table S3). We focused on significant SNPs with a threshold of -log10 (p-value) ≥ 4. For convenience, we selected tagged SNPs for each QTL exhibiting the strongest association with RSA traits (Table 4), yielding 36 QTLs on 17 chromosomes (Figure 4), of which 13 were associated with more than one RSA trait. Among these QTLs, five QTLs, mapped on chromosomes 5B, 6B, 6D, 7A, and 7B, were co-localised with those identified in previous studies (

QTL Clusters for Root System Architecture Traits
QTL clusters refer to the QTLs that were associated with more than one RSA trait. We identified 13 QTL clusters located on chromosomes 2A, 2B, 3A, 3B, 3D, 4A, 5B, 5D, 6B, and 7B. Chromosomes 2A, 3B, and 7B contained two QTL clusters, and chromosomes 2B, 3A, 3D, 4A, 5B, 5D, and 6B harboured a single QTL cluster. One QTL cluster on chromosome 5B (571.23 Mb) was associated with five RSA traits: RL, RA, RV, RT, and NRT. One QTL cluster on chromosome 3D (549.74 Mb) was associated with four RSA traits: RL, RA, RV, and TRL. Two QTL clusters on chromosome 2A (61.28 Mb and 740.01 Mb) and one QTL cluster on chromosome 2B (694.26 Mb) were associated with TRA and TRV. Two QTL clusters on chromosome 7B were associated with TRL and TRA. One QTL cluster on chromosome 3A, two on chromosome 3B, one on chromosome 4A, one on chromosome 5D, and one on chromosome 6B were associated with RT and NRT. Therefore, these QTL clusters showed a pleiotropic effect.

Candidate Genes Associated with RSA
A total of 649 genes were identified in the 13 QTL clusters (Table S4). The genes encoding the bHLH transcription factor, WUSCHEL-RELATED HOMEOBOX 4 (WOX4), ROOT HAIR DEFECTIVE 6 (RHD6), GLABRA 3 (GL3), and ROOTHAIRLESS-LIKE 1 (LRL1), which are closely related to RSA traits, were found in the QTL clusters. We analysed the significant SNPs in thirteen QTL clusters, and identified nine genes with nonsynonymous polymorphisms across the wheat panel (Table S5). The phenotypic variation related to each haplotype of the above nine genes was characterised, and the phenotypic differences between the haplotypes of two genes, TraesCS2A01G516200 with the SNP AX-110564036 (GG/AA) and TraesCS7B01G036900 with the SNP AX-112287343 (CC/TT), were significant (p ≤ 0.05) in more than one RSA trait.
TraesCS2A01G516200 (AX-110564036), mapped to chromosome 2A between 739.8 Mb and 740.8 Mb, encoded a short-chain dehydrogenase/reductase (Figure 5a). AX-110564036 SNP caused an amino acid change from isoleucine (Ile) to valine (Val) at 167 bp in this gene (Figure 5b). The 388 accessions were divided into two haplotypes (GG or AA) according to their genotype at this SNP. We determined that 258 and 93 accessions harboured the GG and AA haplotypes, respectively, whereas the remaining 37 accessions were heterozygous or lacked genotype information. Under the two environments, the accessions carrying the AA allele showed significantly higher values of the RSA traits TRL, TRA, and NRT than those carrying the GG allele (p ≤ 0.05), indicating that the AA haplotype had a positive effect on RSA (Figure 5c). RT-qPCR analysis showed that this gene was preferentially expressed in the root tissues at the jointing and heading stages (Figure 5d). These results showed that TraesCS2A01G516200 plays an important role in the RSA traits of wheat.
TraesCS7B01G036900 (AX-112287343), mapped to chromosome 7B between 35 Mb and 36 Mb, encoded a development and cell death domain protein (Figure 6a). AX-112287343 SNP caused an amino acid change from valine (Val) to alanine (Ala) at 396 bp in this gene (Figure 6b). The 388 accessions were divided into two haplotypes (CC or TT) according to their genotype at this SNP. We determined that 348 and 36 accessions harboured the CC and TT haplotypes, respectively, while the remaining four accessions were heterozygous or lacked genotype information. Under the two environments, the accessions with the CC haplotype exhibited markedly higher RSA trait values for TRL, TRA, TRV, RL, RA, and RV than those with the TT haplotype (p ≤ 0.05), indicating that the CC haplotype had a positive effect on RSA (Figure 6c). The RT-qPCR analysis indicated that this gene was preferentially expressed in the root tissues at the jointing and heading stages (Figure 6d). These results showed that TraesCS7B01G036900 plays an important role in the RSA traits of wheat.

Soil Culture and 660K SNp Chip Are Ideal Methods for GWAS of RSA Traits in Wheat
Bread wheat is a staple crop worldwide, whose RSA traits strongly affect crop development and yield [49]. The spatial distribution of roots affects the absorption of water and nutrients and plant development. Zhu et al. [50] established that the length and number of seminal roots are particularly vital for the uptake of immobile nutrients. Devaiah et al. [51] found that the differences between RSA traits lead to changes in nutrient absorption. In addition, a deeper root system has higher water absorption capacity owing to its better access to soil water [52]. More importantly, White et al. [53] indicated that shoot and root architectures are inherited independently, and could be improved to increase the absorption and utilisation of water and nutrients. Therefore, it is important to analyse the RSA traits and explore the QTLs for these characteristics for further screening of key genes.
Some studies have identified several QTLs for RSA traits in wheat. To conveniently obtain whole roots, plants were usually cultivated in Hoagland's nutrient solution, sand, or agarose gel [45,54]. Hydroponics is the most common approach used to obtain the whole roots in wheat [22,31,32,42]. Although the QTLs for RSA traits that were identified through hydroponic cultivation coincided with those for yield and nutrient absorption in wheat [43,44,55], we were not sure whether the phenotype of RSA traits in hydroponics was similar to that in soil cultivation. Therefore, soil cultivation is reliable to investigate RSA traits, although it is time-consuming and difficult. In this study, 388 wheat accessions were cultivated in soil under OPC and IPC, and the whole roots were sampled and cleaned up at the three-leaf stage to obtain the RSA traits for further analysis. The nine RSA traits showed high phenotypic variation at the seeding stage, indicating that the wheat panel used in this study exhibited high genetic diversity. The nine RSA traits under the two environments were significantly correlated, and the correlation between these RSA traits was significant, indicating that some QTLs may contribute to more than one RSA trait. The mean values of these RSA traits under OPC were significantly higher than those under IPC. We found that the soil under the two environments was comparable, and the mean temperature under IPC was higher than that under OPC. However, the root development was faster under OPC than that under IPC. It may be due to the climatic conditions from 15 October 2019 to 9 November 2019 that were optimal for wheat seedlings in Xinxiang (35.29 • N, 113.83 • E), Henan, China. The broad-sense heritability of the nine RSA traits ranged from 66.78-77.83%, indicating that these traits were suitable for QTL exploration.
Because of the vastness of genomic data, it is expensive to explore QTLs through genome re-sequencing in wheat. A series of high-density SNP chips, including 820K [56], 50K [57], 55K [58], 90K [33], and 660K chips, were designed and utilised for markerassisted breeding in bread wheat. The 660K SNP chip contains 229,266 SNPs in the gene or promoter regions of 66,834 genes, representing 63.52% of the annotated genes in wheat. It is an optimal chip for marker-assisted breeding and has been widely used for QTL exploration for quality traits, agronomic traits, and disease resistance in bread wheat [48]. In the present study, the wheat panel was genotyped for GWAS using the 660K SNP chip.
It is reliable to screen for significant QTLs associated with RSA traits.

QTL Analysis for RSA Traits
A substantial number of studies have performed QTL mapping for RSA traits in different plants, of which 26 studies focused on RSA traits in wheat. They identified 410 QTLs for RSA traits (Table 5). In this study, 36 QTLs with a threshold of −log 10 (p-value) ≥ 4 were identified for nine RSA traits. Among these, five QTLs overlapped with previously reported QTLs, and 31 new QTLs were identified. This indicates that these genomic regions for the RSA traits are reliable and provide new information for further studies. NRT, TRA, TRL, TRV, MRL [22,29,30,32,34,37,38,40,45] This study identified thirteen QTL clusters associated with more than one RSA trait that were distributed on ten chromosomes, indicating that there were co-localisations among these RSA traits. From the relationship among the nine RSA traits, we found that eight RSA traits showed significant positive correlations, and ARD showed a significant negative correlation with NRT, RL RA, RV, and RT in the two environments. From these results, we found that root traits are closely associated with each other, and some traits could be controlled by the same QTL. This might be because of the pleiotropic effects of the same gene or some RSA traits that might have a common genetic basis.

QTL Analysis Reveals the Key Genes for RSA
Key genes regulate RSA and root development. Therefore, the screening of key genes is important for RSA improvement. Zhang et al. [59] found that WOX4, encoded by AT1G46480, plays a central role in cambium development in Arabidopsis roots. Here, we found TraesCS2A01G514000, a homologue of AT1G46480, in the QTL on chromosome 2A, which was associated with TRA and TRV. Menand et al. [60] found that RHD6, encoded by AT1G66470, was a positive regulator for root hair development. TraesCS3B01G007700 was found to be homologous with AT1G66470 and mapped to the QTL of chromosome 3B that was associated with NRT and RT. Bernhardt et al. [61] found that a GL3 protein, encoded by AT5G41315 in Arabidopsis, partially regulates root epidermal cell fates. Here, we found that TraesCS2D01G575600, a homologue of AT5G41315, was located in the QTL of chromosome 2D that was associated with ARD, which might influence root epidermal cell fate. Karas et al. [62] found that LRL1, encoded by AT2G24260, is a partial regulator of root hair development in Arabidopsis. Here, we found that TraesCS7D01G158300, a homologue of AT2G24260, was located in the QTL of chromosome 7D that was associated with ARD. Therefore, these genes were considered candidates for the respective RSA traits.
In addition, we manually screened the significant SNPs for haplotype analysis, and found that nine genes showed non-synonymous polymorphisms across the wheat panel. Two genes, TraesCS2A01G516200 (AX-110564036) and TraesCS7B01G036900 (AX-112287343), have attracted our attention. The different combinations of haplotypes for the two genes showed significant differences for RSA traits. Importantly, the two genes were highly expressed in the root tissues. Functional annotation analysis revealed that TraesCS2A01G516200 encoded a short-chain dehydrogenase/reductase and TraesCS7B01G036900 encoded a development and cell death domain protein. These results indicate that these two genes may be vital candidate genes for root architecture.

Plant Materials and Experimental Treatment
Based on pre-diversity assessment of more than 5000 wheat accessions using KASP markers or 660 K SNP chip, a diversity panel of 388 wheat accessions was collected, including 38 local landraces, 66 exotic cultivars, and 284 cross-derived cultivars in China. The name and origin of each wheat accession are listed in Table S1. Wheat seeds of uniform size were washed thrice with distilled water and germinated on wet filter papers in Petri dishes for germination at 20°C. After 2 days, nine germinating seeds of each accession were evenly planted in one pot, and each pot was thinned to six seedlings after 10 days of sowing. Each pot (height 15 cm; diameter 20 cm) was filled with 4.5 kg of soil (nutritional soil: sieved tillage soil, 1:1) to easily obtain the whole root system.
For OPC, the pots were randomly buried in a field on 15 October 2019 at Xinxiang (35.29 • N, 113.83 • E), Henan, China. To facilitate management, 10 blocks were set in the field, each containing approximately 40 pots. The plants were watered to 70-80% of the field capacity. The seedlings developed to the three-leaf stage after 25 days of sowing with a mean temperature of 14.85 • C and 0 mm of rainfall. For IPC, the pots were randomly placed in a phytotron. The seedlings were grown under 80% relative humidity, a 16 h light (20 • C) / 8 h darkness (18 • C) photoperiod, and 8000 Lux light intensity. The seedlings developed to the three-leaf stage after 30 days of sowing. The roots were cleaned with water. For each wheat accession under the two environments, six biological replications were sampled for further RSA trait measurement.

RSA Trait Measurement
The root samples were placed in plexiglass trays (200 mm × 250 mm) containing 4-6 mm of distilled water. The roots were manipulated to minimise tangling and overlap. Nine RSA traits (Table 6), TRL, TRA, ARD, TRV, NRT, RL, RA, RV, and RT, were measured using a recording scanner (Epson 1680, Suwa, Japan). The images were analysed using Win RHIZO Pro Vision 5.0a (Regent Instruments Inc., Quebec, QC, Canada), which was operationally semi-automated.
where σ 2 G is the genetic variance, σ 2 GE is the genotype × environment, σ 2 e is the error variance, n represents the number of environments, and r represents the number of replications. The mean values were used to analyse the differences using a t-test. The CV was calculated for all the RSA traits [63].

Population Structure
Principal component analysis of the population was performed using the GCTA software [64]. For population structure analysis, a Bayesian model-based (Markov Chain Monte Carlo) clustering approach was used in the STRUCTURE software v2.3.4 with unlinked markers (r 2 = 0) [65]. The K-value, representing the number of subgroups, was set from 2 to 20, and five independent runs for each K were performed with the burn-in length set to 20,000 and iterations set to 10,000. The most likely number of subpopulations was determined using the ∆K method based on the rate of change in LnP (D) between successive K-values [66].

LD Analysis
Genome-wide LD analysis for the A, B, and D genomes was performed using the PLINK software. The squared allele-frequency correlations (r 2 ) were used to estimate the LD by applying pair-wise comparisons among the filtered SNP markers. The values for genomes A, B, and D were plotted against the genetic distance to determine LD decay. The parameters for calculating r 2 were set to -r2 -ld-window-kb 30000 -ld-window 1000 -ld-window-r2 0. A locally weighted polynomial regression (LOESS) curve was drawn to fit the data using second-order local weighted scatter plot smoothing in the R program. The critical distance spanning the QTL was defined based on the intersection point of the fitted LOESS curve with the LD (r 2 = 0.1).

Genotyping Data Using 660K SNp Chip
The wheat panel was genotyped using the Affymetrix Wheat 660K SNP Chip [48], which was provided by the Beijing CapitalBio Technology Company (http://www.capitalbiotech.com, accessed on 8 January 2022). SNP genotyping calling and allele clustering were performed using the polyploid version of the Affymetrix Genotyping Console software. To ensure genotyping accuracy, quality control, and filtering of raw data with missing data > 10%, a minor allele frequency < 5%, or Hardy-Weinberg Equilibrium > 0.01 were performed. The remaining high-quality SNP markers were used for GWAS. The sequence of each polymorphic SNP was identified using Basic Local Alignment Search Tool (BLAST)n 2.2.29+, and the physical positions of all SNPs were determined by blasting them with the wheat reference genome IWGSC RefSeq v1.0. The parameters of the blast result were set to: ratio (AL/QL) ≥ 0.95, identity ≥ 0.95, and E-value < 1 × e −10 . Two genetic diversity parameters, PIC and He, were calculated for each SNP marker and each chromosome using the formulas described by Botsteinet et al. [67] and Nei [68], respectively.

GWAS
GWAS was conducted using a univariate linear mixed model in the GEMMA software [69]. The suggestive threshold for p-value (p = 1/Ne) was calculated using a modified Bonferroni correction (genetic type 1 Error calculator, version 0.2), where Ne represents the effective number of SNP markers [70]. A threshold was established at −log 10 p ≥ 4 to determine the marker-trait association significance. Moreover, the p-value of −log 10 p ≥ 3.5 was used to establish a second, less restrictive threshold. Significant markers from GWAS were visualised in Manhattan plots, and important p-value distributions were visualised using quantile-quantile (Q-Q) plots drawn with the qqman package in R 4.1.1. The phenotypic variance (R 2 ) explained by significant SNPs was evaluated using the GCTA software [64].

Candidate Gene Identification, RNA Extraction, and RT-qPCR Verification
Candidate genes were selected using functional annotation and haplotype analysis. For haplotype analysis, the significant SNPs in QTL clusters were used to identify nonsynonymous polymorphisms. The phenotypic differences between the haplotypes of the nine RSA traits were analysed to identify the candidate genes.
The root, stem, and leaf tissues were sampled from wheat plants of cv. Chinese Spring grown in the field at the jointing and heading stages, respectively. The total RNA was extracted using the TRIzol reagent (TaKaRa, Ohtsu, Japan) according to the manufacturer's instructions. One microgram of the total RNA was used for first-strand cDNA synthesis using the HiScript III First-Strand cDNA Synthesis Kit (Vazyme, Nanjing, China). qPCR was performed on an ABI 7500 real-time PCR system (Applied Biosystems, Forest City, CA, USA) using SYBR Premix Ex Taq Kit (TaKaRa, Japan). Three biological replications were used for each experiment. TaActin3 was used as an internal control [71], and the relative expression value was calculated using the 2 −∆∆Ct method. All primers used for RT-qPCR were listed in Table S6.

Conflicts of Interest:
The authors declare no conflict of interest.