Genome-Wide Association Study Reveals Candidate Genes for Control of Plant Height, Branch Initiation Height and Branch Number in Rapeseed (Brassica napus L.)

Plant architecture is crucial for rapeseed yield and is determined by plant height (PH), branch initiation height (BIH), branch number (BN) and leaf and inflorescence morphology. In this study, we measured three major factors (PH, BIH, and BN) in a panel of 333 rapeseed accessions across 4 years. A genome-wide association study (GWAS) was performed via Q + K model and the panel was genotyped using the 60 k Brassica Infinium SNP array. We identified seven loci for PH, four for BIH, and five for BN. Subsequently, by determining linkage disequilibrium (LD) decay associated with 38 significant SNPs, we gained 31, 15, and 17 candidate genes for these traits, respectively. We also showed that PH is significantly correlated with BIH, while no other correlation was revealed. Notably, a GA signaling gene (BnRGA) and a flowering gene (BnFT) located on chromosome A02 were identified as the most likely candidate genes associated with PH regulation. Furthermore, a meristem initiation gene (BnLOF2) and a NAC domain transcriptional factor (BnCUC3) that may be associated with BN were identified on the chromosome A07. This study reveals novel insight into the genetic control of plant architecture and may facilitate marker-based breeding for rapeseed.


INTRODUCTION
Rapeseed (Brassica napus L., AACC, 2n = 38) is a globally significant oilseed crop used in the production of vegetable oil and oil meals, and ∼95% of its total cultivation is in Asia, Europe, and North America. According to a new report from the United States Department of Agriculture (USDA), rapeseed production in 2015/16 is estimated at 66.9 million tons, down 4.7% from the previous year. Harvested area is estimated at 32.9 million hectares, down 5.5% from the previous year (http://www.pecad.fas.usda.gov/). To meet growing demands for both edible products and biofuel, it is critical to increase rapeseed productivity through efficient breeding.
Rapeseed productivity is determined by three components: number of siliques per plant (NSP), number of seeds per silique (NSS), and seed weight per plant (Clarke and Simpson, 1978).
Simultaneously, rapeseed plant architecture is a crucial determinant of its yield and is itself determined by plant height (PH) and branch number, as well as inflorescence morphology, which indirectly influences its yield by affecting NSP (Qiu et al., 2006;Chen et al., 2014). Plant height was first selected to improve crop yields: taller plants were more easily affected by lodging and therefore showed decreased yield, while dwarf plants showed resistance to lodging and their yield could be further increased by the use of nitrogen fertilizers. In the late 1960s and early 1970s, the dwarf gene was rapidly adopted in wheat and rice to cultivate new varieties, which greatly increased their yield and precipitated the "Green Revolution" (Peng et al., 1999;Evenson and Gollin, 2003). In addition, most domesticated crops also contain the plant domestication gene TB1 (Doebley et al., 1997;Takeda et al., 2003;Aguilar-Martinez et al., 2007;Lewis et al., 2008). Recently, ideal plant architecture (IPA) has been proposed as a plant type with increased yield potential (Jiao et al., 2010;Miura et al., 2010).
In Arabidopsis and rice, the mechanisms regulating PH are well-known, and many phytohormones, including Gibberellins (GAs), Brassinosteroids (BRs), and Strigolactones (SLs) participate in this process. Genes related to phytohormone biosynthesis and signal transduction contain mutations that affect internode elongation, which in turn regulates PH (Ikeda et al., 2001;Ishikawa et al., 2005;Sun et al., 2010;Clouse, 2011;Jiang et al., 2013;Zhou et al., 2013). During the floral transition, axillary meristems are transformed into branch meristems to form branches (Teo et al., 2014). Phytohormone, including Auxin, Cytokinins (CKs), and SLs biosynthesis, transport, and signaling genes also play important roles in branch formation (Ferguson and Beveridge, 2009;Janssen et al., 2014), as do floral meristem identity genes (i.e., MYB family RAX proteins, ALOG family proteins, other meristem identity genes, flowering time genes, etc.; Liljegren et al., 1999;Hiraoka et al., 2013;Liu et al., 2013). In Arabidopsis, most integrated regulatory networks controlling the formation of branching are regulated by the gene SHOOT MERISTEMLESS (STM) (Long et al., 1996;Lenhard et al., 2002). Despite significant advances in other plants, our understanding of the molecular mechanisms of plant architecture regulation in rapeseed remains limited.
Rapeseed (B. napus) and the model plant Arabidopsis are members of the Brassicaceae. Unlike Arabidopsis thaliana, B. napus is polyphyletic, and is related to B. rapa (AA) and B. oleracea (CC) (Allender and King, 2010). So far, genomes have been sequenced for B. rapa (AA), B. oleracea (CC), and B. napus (AACC) (Wang et al., 2011;Chalhoub et al., 2014;Liu et al., 2014). B. napus is a young allopolyploid and most orthologous genes are duplicated compared to the respective progenitor genomes (Chalhoub et al., 2014). Consequently, obtaining a detailed characterization of each copy is difficult.
Genetic mapping studies have employed bi-parental mapping populations to identify many quantitative trait locus (QTLs) for yield traits in rapeseed. However, most of these QTLs are localized to a large interval (10-20 cM) on chromosomes (Butruille et al., 1999;Quijada et al., 2006;Udall et al., 2006;Chen et al., 2007;Mei et al., 2009;Shi et al., 2009). The genomewide association study (GWAS) has become a powerful tool to identify multiple related candidate genes regulating important traits in crops (Huang et al., 2010;Brachi et al., 2011;Li H. et al., 2013). This method is performed by scanning a genome using abundant single-nucleotide polymorphisms (SNPs) representing broad genetic variability. Presently, the GWAS in rapeseed has been used to identify loci and candidate genes related to oil content, plant architecture, flowering time, and other yield traits (Schiessl et al., 2015;Liu et al., 2016;Lu et al., 2016;Sun et al., 2016a,b;Sun F. M. et al., 2016;Wang et al., 2016;Xu et al., 2016).
In this study, we report a GWAS for plant architecture traits (PH, BIH, and BN) in rapeseed using the 60 K Brassica Infinium SNP array on a panel with 333 accessions that cover a broad range of genetic diversity. We further identify several possible candidate genes for the three traits underlying these QTLs by LD decays harboring significant SNPs. Although some candidate QTLs are consistent with those found in previous studies, we have also identified novel QTLs for these traits. Moreover, the SNPs and candidate genes in our study may facilitate markerbased breeding to improve plant architecture in order to increase rapeseed yield.

Plant Materials and Trait Measurement
The association population used in this study consisted of 333 diverse rapeseed accessions (20 are winter type, 308 are semiwinter type and 5 are unknown) collected as part of a recently published study (Sun F. M. et al., 2016). The rapeseed accessions were grown in two replicates over the course of 4 years (2012)(2013)(2014)(2015) in Wuhan, China (E 114.32 • , N 30.52 • ). Individuals from germplasm populations were genotyped using the 60 K Brassica Infinium SNP array. Accessions were grown in plots of 2 rows, with 10-15 plants in each row. Between eight and ten plants were selected from the center of each row to measure the three traits at maturity: PH, BIH, and BN.
The length of the plant from the base of the stem to the tip of the main inflorescence was noted as PH. The length from the base of the stem to the first primary branch base was recorded as BIH. We also measured the number of primary branches arising from the main shoot (BN). Correlation analyses of PH, BIH, and BN for the association panel were performed in R3.3.0 (Ihaka and Gentleman, 1996).

In silico Mapping and Linkage Disequilibrium (LD) Analysis
The probe sequences of 52,157 SNPs were used to perform a BLASTN search through the B. napus genome database (Chalhoub et al., 2014). Only top blast-hits with an e-value threshold of e −15 were used. Some SNPs were excluded due to a match in BLAST with multiple loci with identical scores. SNPs that showed a minor allele frequency (MAF) of <5% or that showed a maximal missing frequency of >2% were excluded (Brachi et al., 2011), remaining 32,297 SNPs identified, covering all 19 chromosomes of the B. napus genome and providing approximately one SNP per 25 kb. An estimate of the linkage disequilibrium (LD) was calculated using squared allele frequency correlations (r 2 ) between all pairs of SNPs using TASSEL 5.2.28 (Bradbury et al., 2007).

Genome-Wide Association Study
All the SNPs were used to calculate the principal component analysis (PCA) matrix using the GCTA tool, and a subset of 3,571 SNPs (MAF ≥ 0.2) evenly distributed across the whole genome (every 150 k) were selected to perform population structure (Q) and relative kinship analysis (K). STRUCTURE v2.3.4 was used to calculate a Q matrix, the putative number of genetic groups (k-value) setting from 1 to 10 with five independent runs (Pritchard et al., 2000). The length of the burn-in period and the number of Markov Chain Monte Carlo (MCMC) replications after burn-in were set to 50,000 and 100,000 iterations, respectively. The relative kinship K matrix was performed by the SPAGeDi software package (Hardy and Vekemans, 2002;Yang et al., 2011). The k method was used to determine the most likely number of groups or subpopulations, as described by Evanno et al. (2005).
Association analysis was performed using TASSEL 5.2.28 (Bradbury et al., 2007), using a mixed linear model (MLM) to calculate the association in all analyses, incorporating Q matrix/PCA and kinship data (K) (Zhao et al., 2011). The MLM was applied using default settings (P3D for variance component analysis and compression set to optimum level). For MLM (Q + K), the significance threshold for significantly associated markers was set to p ≤ 4.06 × 10 −4 [−log 10 (p-value) = 3.39].

Candidate Gene Mining
To identify candidate genes, local LD decay was calculated using TASSEL5.2.28 to capture flanking regions of up to 450 kb on either side of significant SNPs, with a cut-off value of r 2 = 0.2. Gene sequences that correspond to putative orthologs in A. thaliana were based on GO annotations (http://www.arabidopsis. org/index.jsp). The genes with GO terms for auxin, GA, IAA, SL, CK, and flowering time were highlighted, and the closest one of these genes was considered as the most likely candidate.

Phenotypic Variation among Accessions for PH, BIH, and BN
Extensive phenotypic variations for PH, BIH, and BN were observed in the 333 accessions that were grown over the course of 4 years (2012)(2013)(2014)(2015). PH varied from 86.2 to 206.0 cm, with 1.6 to 2.4-fold variations across the 4 years ( Figure 1A, Table 1). BIH varied from 6.9 to 122.0 cm, with 5.9 to 17.3-fold variations across the 4 years ( Figure 1B, Table 1). BN varied from 1.7 to 15.0, with 3.5 to 8.6-fold variations across the 4 years ( Figure 1C, Table 1).
In addition, we assessed these three traits for any significant correlation. PH and BIH were significantly correlated across the 4 years with a Pearson's correlation coefficient of 0.3 to 0.8 (Table S1). This suggests that these two traits may show genetic linkage or that some genes have pleiotropic roles in controlling these phenotypes. However, no other significant correlations were found in our accession panel, between the PH and BN, or between BN and BIH (Table S1).

SNP Quality Control, Performance, and in silico Mapping
SNP genotyping was performed using the Brassica 60 K SNP array. We blasted the SNP probe sequences to the B. napus genome database (http://www.Genoscope.cns.fr/brassicanapus/) and a total of 34,292 SNP markers (65.7%) were selected to genotype the panel of 333 accessions ( Table 2). The C04 linkage group had the most SNPs (3,182), with a marker density of one per 20.89 kb, whereas the C05 linkage group had the least SNPs (1,085), with a marker density of one per 39.8 kb ( Table 2). Altogether, the overall mapping results demonstrate the high quality of the genotyping in this study.

Linkage Disequilibrium
We calculated LD in 333 accessions using the parameter r 2 with the 34,292 SNP markers. In this analysis, three chromosomes  (A08, C02, and C07) showed strong LD, with distances ranging from 1,300 to 1,710 kb; while others exhibited modest LD, with distances ranging from 100 to 825 kb when r 2 = 0.2 ( Table 2).
The average LD decay of chromosomes A and C were of 190 and 700 kb respectively, while the average LD decay of the whole genome (A + C chromosomes) was of 450 kb when r 2 = 0.2 (Sun F. M. et al., 2016; Figure 2A).

Genome-Wide Association Study
We performed GWAS with mixed linear modeling (MLM) across 4 years. For the groups-populations, we created a Q matrix (k = 2) for all 333 accessions, in order to obtain the highest k-value ( Figure 2B). In addition, minor differences were observed among the PCA, Q, K, Q + K, and PCA + K models, and the Q + K model is an appropriate model (Figure 3). Thus, the Q + K model was used for GWAS. As in other GWAS studies, we set a threshold of −log 10 P = 3.39 as a significant association ( Figure S1; Benjamini and Hochberg, 1995;Wang et al., 2016). Furthermore, QQ plots of PH in 2012 and 2013 were not linear. As a result, the SNPs showing the 50 highest −log 10 P-values in each year were compiled for the PH trait for these 2 years to identify potentially significant SNP clusters (Figure 4; Ueda et al., 2015). Overall, for all traits analyzed, a total of 158 SNPs associated with PH, BIH, and BN were identified (Table S2). Furthermore, LD blocks (450 kb) harboring the repeating SNPs across the 4 years (and detected in at least two of those years) were identified as regions containing putative candidate loci. This led to the selection of 38 SNPs (Table 3), which in turn could be merged into seven, four, and five loci for PH, BIH, and BN, respectively (Figures 4-6, Table 3, Figures S2-S4).
For PH, we detected nine GWAS peak SNPs in the seven loci located on chromosomes A02, A06, C03, and A07, and five of which were detected in at least 2 years (Figure 4, Table 3). These seven peak SNPs explained 35.56% of the total phenotypic variance, with the largest contribution (R 2 ∼9.79%) from SNP Bn-A02-p9505646 (Table 3). For BIH, four peak SNPs were detected on chromosomes A02, A07, A08, and A09, and three of these SNPs were detected in at least 2 years (Figure 5,  Table 3). These three SNPs explained 27.96% of the phenotypic variance (Table 3). SNP Bn-A02-p9505646 was detected in both PH and BIH traits and it accounted for 9.36% of the total phenotypic variance in BIH (Table 3). For BN, five peak SNPs were detected on chromosomes A01, A03, A07, C04, and C09 (Figure 6), two of which were detected in at least 3 years, with the largest contribution (R 2 ∼11.54%) from SNP Bn-A07-p21413042 (Table 3).

Allele-Specific SNP Markers Correlated with Physical Traits (PH, BIH, and BN)
The GWAS identified several significant SNP markers in the regulation of three traits (Table 3). To find useful SNP markers for marker-based breeding, we estimated the allelic effect of the peak SNPs across these three traits. The marker Bn-A02-p9505646 (G/A) of the GG allele showed the largest contribution (R 2 ∼9.79%) to PH regulation, measuring 21.5 cm more than those of the AA allele (P < 0.05; Figure 7A). Apart from the Bn-A02-p9505646 marker, individuals with the CC allele marker Bn-A09-p9669847 showed the largest contribution to BIH regulation (R 2 ∼9.87%) and measured 11.4 cm more than those of the AA allele (P < 0.05; Figure 7B). The average BN of individuals with the AA allele with marker Bn-A03-p7326202 showed the higher contribution (R 2 ∼8.61%) and was significantly greater than GG alleles on chromosome A3 ( Figure 7C). These SNP markers may permit marker-based breeding for rapeseed plant architecture.

Identification of Candidate Genes
We searched for candidate genes from the genomic regions with significant SNPs according to GO annotations from the Arabidopsis database (http://www.arabidopsis.org/index. jsp). Notably, 31, 15, and 17 candidate genes for PH, BIH, and BN were located in each region, respectively (Table S3).
For BIH, the C2H2 zinc finger protein, TCP 12, and SOB five-like 2 (BnaA02g13870D, BnaA02g14010D, and BnaA02g14060D) were also identified in the candidate region around the peak SNP Bn-A02-p9505646 [−log 10 (p) = 4.70] on chromosome A02 (Table 4). On chromosome A07, a CLAVATA3 protein and a LOF2 transcription factor (BnaA07g24210D and BnaA07g24240D) were identified downstream from the peak SNP Bn-A07-p15758978) [−log 10 (p) = 4.88] ( Table 4). Their homologs in Arabidopsis are also involved in meristem initiation and maintenance (Lee et al., 2009). In addition, on chromosome A09, we found a peak SNP Bn-A09-p9669847 with the highest probability of marker-trait association, with −log 10 (p) = 5.04, and with the largest contribution to BIH variation, with R 2 = 9.87%. Notably, a MATE efflux family protein (BnaA09g14730D) involved in determining the rate of organ initiation was identified close to the peak SNP (Table 4).

DISCUSSION
In our study, a MLM was used to calculate the association in the GWAS analysis, while incorporating a Q matrix and kinship data to control the false discovery rate, this model was an improvement relative to the naïve GLM (Figure 3; Larsson et al., 2013). Although it has been suggested that PCA or PCA + K models may reduce false positives in GWAS analyses Zhu et al., 2008;, many studies have used the Q + K model because it permits control of the false discovery rate (Cai et al., 2014;Sun F. M. et al., 2016;Wang et al., 2016). In our assays, even the QQ plots of PH deviated from the expected distribution, while the Manhattan plots showed clearly defined peaks, indicating that MLM might be the best approach to detect the main QTLs for our traits (Figures 4-6, Figure S1). According to previous studies, we defined the SNPs which less than the p-value of 4.06 × 10 −4 [−log 10 (p-value) = 3.39, Q + K] or in the top 50 SNPs in PH control as the candidate SNPs (Ueda et al., 2015;Wang et al., 2016). The LD decay is the key factor for mining candidate genes in a GWAS study (Yu and Buckler, 2006). We observed that LD decay in the A-and C-subgenomes is different in our panel. The LD decay in the A-subgenome was generally shorter than those in the C-subgenome, while the LD decays on A08, C02, and C07 were beyond 1.3 M. Our results for LD decay rates are similar to those found in recent GWAS studies on rapeseed (Qian et al., 2014;Liu et al., 2016), and therefore, we hold that the LD decay in our study is reliable to identify candidate genes.

Novel Genetic Control of Three Plant Architecture Traits (PH, BIH, and BN) in Rapeseed
GWAS was used to detect QTLs of PH, BIH, and BN, which are main factors regulating plant architecture in rapeseed. According to the marker-trait association SNPs, we identified seven QTLs for PH, four for BIH, and five for BN (Table 3). In the PH candidate loci, these QTLs were distributed on six chromosomes and some QTL regions were similar to those identified in previous studies. Among them, the QTLs at 1.2 and 6.4 M on chromosome A02 overlapped regions detected by Udall et al. (2006) ;Li X. N. et al. (2013), and Sun et al. (2016b). The QTL at 18.0 M on chromosome A07 was close to the QTLs detected by Udall et al. (2006) and . The QTL near the top of chromosome C03 was also detected by Ding et al. (2012) and Shi et al. (2009). Additionally, the QTL at 6.2 M on chromosome A03 was found by Mei et al. (2009) and Li X. N. et al. (2013), and also corresponds closely to the findings of . The QTL at 45.8 M on chromosome C02 corresponds closely to the QTL on the bottom of C02 detected by Udall et al. (2006). These results indicate that these QTLs were stable as they were detected by different map populations and in different environments using distinct methods of analysis, and they may also have been selected in rapeseed breeding. Notably, the QTL at 19.7 M of A06, which was detected in 2 years has not been found in previous studies, and may be a novel QTL.
For the BN candidate loci, only one QTL was located at 6.6 M on chromosome A03, as shown previously in other studies (Shi et al., 2011;. Interestingly, we found that all the other four QTLs were detected in a rapeseed accession that has a stable BN but increased branching. In our study, most of the accessions are varieties cultivated from different    Previously, only one study identified QTLs associated with the height of the lowest primary effective branch (HPB) using two populations in rapeseed, and it identified 10 QTLs distributed on chromosomes A02, A07, A08, C04, C06, and C07 (Chen et al., 2007). In our study, BIH was measured from the ground level up to the base of the lowest primary effective branch, and four QTLs were detected in our natural population. One of these QTLs was located at 18.0 M on chromosome A08, which corresponds to the QTL detected on the bottom of the same chromosome by Chen et al. (2007). Two more QTLs were found on corresponding chromosomes A02 and A07. However, their specific regions did not overlap. To verify the reliability of these QTLs, further investigations of this trait must be undertaken.
In our study, the QTLs at 6.3 M on chromosome A02 and at 18.0 M on chromosome A07 influenced PH and BIH, and a significant correlation was found between these two traits during our correlation assay. In the A02 QTL regions, PH and BIH traits shared three SNPs and the two traits were also co-localized, which suggests pleiotropic regulation of a single gene or the existence of closely linked genes. In contrast, in the A07 QTLs, the SNP associated with the PH QTL was separated from the related SNP associated with the BIH QTL by about 700 kb, while the LD (r 2 ) between the two SNPs was 0.27, indicating that they are separate QTLs. In addition, QTLs related to PH on chromosome A03 and A07 were found close to the QTLs related to BN on the same chromosomes, at a distance of 500 kb and 1 M, respectively. The QTL at 17.7 M on chromosome A07 related to BIH regulation was close to the QTL related to BN, at ∼1.3 M. These distances are also beyond the average LD decay and the LD (r 2 ) between the two peak SNPs for each group were relatively low, suggesting that they are separate QTLs.

Candidate Genes For PH and BN Regulation
In Arabidopsis, rice, and other plant species, BR, IAA, GA, and SL biosynthesis and signaling pathways are known to regulate PH, while mutations in most of these genes cause dwarf phenotypes (Wang and Li, 2008;Sun et al., 2010;Clouse, 2011). Additionally, genes involved in flowering time also play a role in PH regulation, such as DTH8 and Ghd7 in rice (Xue et al., 2008;Wei et al., 2010), and Flt-2L in wheat (Chen et al., 2009). In rapeseed, Mei et al. (2009) found an overlap between a flowering time locus and PH QTLs. Moreover, in some genes regulating cell wall formation, like cellulose synthase genes, mutations affecting cell elongation lead to dwarf phenotypes . Therefore, we defined PH regulation candidate genes in rapeseed based on similarity to orthologs of the QTL regions. Using the average LD decay, 5 genes (16.7%) are related to BR biosynthesis or signaling pathways, 6 genes (20%) are related to GA, 2 genes (6.7%) are associated with Auxin, 4 genes (13.3%) are involved in cell wall formation, 4 genes (13.3%) are related to flowering time regulation, and 2 gene (6.7%) are related to trichome branching development in Arabidopsis (Bischoff et al., 2010; Table S3). Notably, the RGA gene, encoding a DELLA protein, displays a signal close to Bn-A02-p9505646, which itself is close to the Bn-A02-p9610453 marker identified in Sun et al. (2016b). Furthermore, we found that some flowering time QTLs (BnFT) are also located in the same region of chromosome A02 . In Arabidopsis and rice, most studies showed that FLOWERING LOCUS T (FT) and the rice FT homolog HEADINGDATE3a (Hd3a) gene affected flowering and PH (Kardailsky et al., 1999;Kobayashi et al., 1999;Tamaki et al., 2007). In addition, overexpressed the FT in tobacco also exhibited the dwarf phenotype (Lewis and Kernodle, 2009). Thus, we deduce that BnRGA (BnaA02g12260D) and BnFT (BnaA02g12130D) are main candidate genes on chromosome A02 that are involved in PH regulation.
Several transcription factors (TF) and hormones (Auxin, CKs, and SLs) also affect shoot outgrowth in plant architecture (Yang and Jiao, 2016). We compared these orthologs to our main candidate genes for BN. Four genes (23.5%) are related to CK biosynthesis or signaling pathways, 6 genes (35.3%) are related to IAA, 3 genes (17.6%) are associated with flowering or flower identity formation, and one gene is involved in the SL signaling pathway. Many recent studies report that SL is functionally important in BN and PH regulation; however, only one candidate gene was identified in our GWAS. This may be due to a low density of SNP markers, or due to the candidate genes being farther away from the QTL locus than the average LD decay, or due to the confounding of the population structure. In Arabidopsis, LATERAL ORGAN FUSION2 (LOF2) encodes a MYB-domain transcription factor that functions in both lateral organ separation and axillary meristem formation, in part through interaction with CUC2, CUC3, and STM (Raman et al., 2008;Lee et al., 2009). In novel loci on chromosome A07, genes BnLOF2 (BnaA07g26170D) and BnCUC3 (BnaA07g28210D) were identified as the most likely candidate genes and may play roles in both organ separation and axillary meristem formation.
Although several candidate genes were identified following the GWAS, the function of these genes remains ambiguous. In the future, it is necessary to isolate target genes by cloning these putative genes or by developing mapping populations for the QTLs, and to illustrate their functional role in relative trait regulation by transformation experiments.

AUTHOR CONTRIBUTIONS
MZ, WH, and HW conceived and designed the research; HL performed GWAS; XW, MT, HY, XL, JL, and XS characterized the agronomic traits; MZ, CP, HL, and JX analyzed the data; and MZ wrote the manuscript. All authors read and approved the final manuscript.