Functional Mapping of Genes Modulating Plant Shade Avoidance Using Leaf Traits

Shade avoidance syndrome (SAS) refers to a set of plant responses that increases light capture in dense stands. This process is crucial for plants in natural and agricultural environments as they compete for resources and avoid suboptimal conditions. Although the molecular, biochemical, and physiological mechanisms underlying the SAS response have been extensively studied, the genetic basis of developmental variation in leaves in regard to leaf area, petiole length, and leaf length (i.e., their allometric relationships) remains unresolved. In this study, with the recombinant inbred line (RIL) population, the developmental traits of leaves of Arabidopsis were investigated under two growth density conditions (high- and low-density plantings). The observed changes were then reconstructed digitally, and their allometric relationships were modelled. Taking the genome-wide association analysis, the SNP genotype and the dynamic phenotype of the leaf from both densities were combined to explore the allometry QTLs. Under different densities, leaf change phenotype was analyzed from two core ecological scenarios: (i) the allometric change of the leaf area with leaf length, and (ii) the change of the leaf length with petiole length. QTLs modulating these two scenarios were characterized as ‘leaf shape QTLs’ and ‘leaf position QTLs’. With functional mapping, results showed a total of 30 and 24 significant SNPs for shapeQTLs and positionQTLs, respectively. By annotation, immune pathway genes, photosensory receptor genes, and phytohormone genes were identified to be involved in the SAS response. Interestingly, genes modulating the immune pathway and salt tolerance, i.e., systemic acquired resistance (SAR) regulatory proteins (MININ-1-related) and salt tolerance homologs (STH), were reported to mediate the SAS response. By dissecting and comparing QTL effects from low- and high-density conditions, our results elucidate the genetic control of leaf formation in the context of the SAS response. The mechanism with leaf development × density interaction can further aid the development of density-tolerant crop varieties for agricultural practices.


Introduction
Over the course of evolution, plants have acquired an avoidance strategy to adjust their growth when exposed to suboptimal conditions, similar to an escape strategy in animals. Under highly shaded conditions, most plants tailor their photosynthetic strategy to decreased light levels while at the same time expressing shade-avoidance syndrome (SAS). Plants exhibit SAS by positioning leaves into well-lit microsites; accelerating the elongation of hypocotyls, internodes, and petioles; and elongating shoots away from the shade of neighboring plants to reach light. This plays out both in nature and in agricultural According to the advantage of functional mapping, we design two density conditions and incorporate the dynamic petiole and leaf traits from two density sets into a functional mapping framework. From a methodological viewpoint, this strongly support the detection of SAS QTL, but through a QTL × density interaction approach. Taking an Arabidopsis recombinant inbred line (RIL) as the mapping population, we propose the use of bivariate functional mapping to study the genetic architecture of SAS in this study. In particular, how genes modulate the dynamic covariation of petiole with leaf length, how leaf length covaries with leaf area, and how genes mediate the interaction of leaf development and density environment were all comprehensively analyzed in this study.

Shape Reconstruction
We used the shape reconstruction method described by Fu et al [16]. On principle, with the known size of a black/white square, the coordinate system can be build using the corner point of the checkerboard grid. Figure 1 shows the results for two types of Arabidopsis leaf, one with a very short petiole ( Figure 1A) and the other with a tenuous petiole ( Figure 1B). Using this shape reconstruction method, a series of coordinates along the leaf boundary were determined for both the leaves ( Figure 1C,D), which match the photographic images. Figure 1C,D fully possessed the feature of the petiole base from Figure 1A,B. Additionally, the shape with the leaf tip and other parts are reserved, suggesting a strong power with a such method. After the shape reconstruction, the traits of leaf length, leaf area, and petiole length can further be quantified.
these trait relationships and, thus, deepen our understanding of developmental regulation of SAS in leaves [13,14].
According to the advantage of functional mapping, we design two density conditions and incorporate the dynamic petiole and leaf traits from two density sets into a functional mapping framework. From a methodological viewpoint, this strongly support the detection of SAS QTL, but through a QTL × density interaction approach. Taking an Arabidopsis recombinant inbred line (RIL) as the mapping population, we propose the use of bivariate functional mapping to study the genetic architecture of SAS in this study. In particular, how genes modulate the dynamic covariation of petiole with leaf length, how leaf length covaries with leaf area, and how genes mediate the interaction of leaf development and density environment were all comprehensively analyzed in this study.

Shape Reconstruction
We used the shape reconstruction method described by Fu et al [16]. On principle, with the known size of a black/white square, the coordinate system can be build using the corner point of the checkerboard grid. Figure 1 shows the results for two types of Arabidopsis leaf, one with a very short petiole ( Figure 1A) and the other with a tenuous petiole ( Figure 1B). Using this shape reconstruction method, a series of coordinates along the leaf boundary were determined for both the leaves ( Figure 1C, 1D), which match the photographic images. Figure 1C and 1D fully possessed the feature of the petiole base from Figure 1A and 1B. Additionally, the shape with the leaf tip and other parts are reserved, suggesting a strong power with a such method. After the shape reconstruction, the traits of leaf length, leaf area, and petiole length can further be quantified.

Longitudinal Phenotypes of Leaf Traits under Different Density Conditions
Next, using 1680 seedlings, we depicted the phenotypic changes in A, L, and P at the progeny level ( Figure 2). Taking all eight timepoints for low-density and seven timepoints for high-density conditions, all three traits showed S-shaped curves during leaf development. In addition, curve-fitting at the population level and between the two conditions, indicated a high variance for all three traits. Overall, as compared with the Plants 2023, 12, 608 4 of 17 early stage, the high-density condition led to higher growth rates and larger traits in later stages, supporting the general understanding that the elevation of plant density will cause the alteration of light in quality and quantity, to affect the plant growth or elongation synthetically. indicated a high variance for all three traits. Overall, as compared with the early stage, the high-density condition led to higher growth rates and larger traits in later stages, supporting the general understanding that the elevation of plant density will cause the alteration of light in quality and quantity, to affect the plant growth or elongation synthetically.
Low-density conditions led to more variance among phenotypes at nearly all timepoints, suggesting an activated SAS response. As shown in Figure 2, by comparing the mean curve from the two-density condition, their crossover points were located at 4.05, 3.67, and 3.23 weeks for the leaf area (Figure 2A), leaf length ( Figure 2B), and petiole length ( Figure 2C), respectively. The earliest crossover for the petiole trait might imply the sensitivity of petiole growth regulation when facing the changing density condition, suggesting that petiole growth regulation is sensitive to density conditions, and this trait may be the first to trigger the process of SAS. Figure 2. Dynamic growth pattern of leaf traits under high-density and low-density conditions. (A)-(C), respectively, present the ontogenic growth curve of the leaf area, leaf length, and petiole length. Yellow curves denote trait patterns from a low-density environment while blue curves denote that from a high-density environment. The two thick lines represent the fitted mean growth curve, while thin lines represents the observed growth phenotype with each progeny in the population. The vertical line indicates the timepoint when the two mean curves for the two densities crossed. Figure 2. Dynamic growth pattern of leaf traits under high-density and low-density conditions. (A-C), respectively, present the ontogenic growth curve of the leaf area, leaf length, and petiole length. Yellow curves denote trait patterns from a low-density environment while blue curves denote that from a high-density environment. The two thick lines represent the fitted mean growth curve, while thin lines represents the observed growth phenotype with each progeny in the population. The vertical line indicates the timepoint when the two mean curves for the two densities crossed.

Mathematical Modeling of Leaf Allometry
Low-density conditions led to more variance among phenotypes at nearly all timepoints, suggesting an activated SAS response. As shown in Figure 2, by comparing the mean curve from the two-density condition, their crossover points were located at 4.05, 3.67, and 3.23 weeks for the leaf area (Figure 2A), leaf length ( Figure 2B), and petiole length ( Figure 2C), respectively. The earliest crossover for the petiole trait might imply the sensitivity of petiole growth regulation when facing the changing density condition, suggesting that petiole growth regulation is sensitive to density conditions, and this trait may be the first to trigger the process of SAS.

Mathematical Modeling of Leaf Allometry
To verify the fitness between the allometric growth equation and the measured phenotypes, the R-squared value for plants grown under each condition was evaluated using the least squares method. Under high-density conditions, more than 71% of progenies had R-squared values greater than 0.73 for the P-A relationship. For the L-A and P-L relationships, more than 93% and 92% of progenies had R-squared values greater than 0.96 and 0.92, respectively. Under low-density conditions, 75%, 92%, and 90% of progenies had goodness-of-fit values of 0.75, 0.94, and 0.95, respectively. These results indicate that the allometry between P and A was not strong. Therefore, we subsequently focused on analyzing allometry for the latter two pairs of traits, i.e., L-A and P-L. Figure 3 shows the allometric patterns for L-A and P-L during development. L-A has a concave-shaped curve ( Figure 3A,B), whereas that for P-L is convex ( Figure 3D,E). Although the allometric pattern for each progeny varied, those associated with each density were remarkably similar. The data in Figure 3C indicate that, for a given lamina area, high-density conditions led to slightly longer leaves, indicating that leaf elongation occurs to ensure light acquisition in such environments. Nevertheless, petiole lengths were nearly equal between the two densities ( Figure 3F). This suggests that, under high density, petiole elongation does not necessarily account for the increase in leaf length. This coincides with the nonsignificant goodness-of-fit of the allometry data for P-A. Considering the substantial role of the petiole in changing the leaf position and the considerable variation in P-L under both densities, the allometric relationships between P and L can provide important insights into the genetic dissection of the SAS mechanism. By estimating the power coefficients from Equation (1), how one trait scales with the other can be illustrated, and the scaling pattern between the two densities can be compared.

QTLs for SAS
By integrating phenotypes from the high-and low-density groups into the QTL mapping framework, a bivariate functional mapping model incorporating equation (1) was implemented. We tested how QTL regulates the allometric scaling of L-A and P-L during

QTLs for SAS
By integrating phenotypes from the high-and low-density groups into the QTL mapping framework, a bivariate functional mapping model incorporating Equation (1) was implemented. We tested how QTL regulates the allometric scaling of L-A and P-L during leaf ontogeny. We argue that L-A and P-L allometry has different developmental significance for leaves. For L-A, at any given length, leaves with a smaller A are narrower, and this might lead to better photosynthetic ability than larger leaves, increasing competitiveness under crowded and shaded conditions.
To help the interpretation of QTLs, we performed a gene enrichment analysis using the combined genes of two set of analysis. The enrichment result based on Gene Ontology (GO) showed six classes of functions with p-values less than 0.05: red, far-red light phototransduction (GO:0009585, 2.56 × 10 −7 ), gibberellin 20-oxidase activity (GO:0045544, 3.28 × 10 −5 ), shade avoidance (GO:0009641, 0.0014), response to light stimulus (GO:0009416, 0.038), leaf shaping (GO:0010358, 0.021), regulation of defense response (GO:0031347, 0.015). In comparison with the background gene number, 8/32, 4/8, 4/26, 2/12, 2/28, 5/83 genes were involved in these functions, respectively. In particular, for leaf development, the detailed annotation of 'leaf pavement cell development', 'leaf shaping', and 'leaf vascular tissue pattern formation' were found. For class 6, the 'defense response to bacterium', 'regulation of systemic acquired resistance', 'defense response to fungus', 'response to UV-B', and 'defense response to bacterium, incompatible interaction' were found. All these have pointed out the importance of leaf development, immune system, far-red light regulation to the SAS, and have also promoted the screening of key genes for the further dissection of SAS.

ShapeQTLs: How Leaf Area Scales with Leaf Length
Supporting the general understanding that the elevation of plant density will cause the alteration of light in quality and quantity, this will affect the plant growth or elongation synthetically. Of these, 280, 82, 86, 90, and 7 were distributed between chromosomes 1 to 5 ( Figure 4A). All were annotated using the TAIR website (www.arabidopsis.org, accessed on 1 October 2022) according to their genomic position, resulting in 30 annotated functional proteins. In terms of protein function, there were fourteen photosensory receptor genes, four phytohormone genes, six immune pathway genes, and one leaf shape-related gene (Table 1), which implies that SAS is regulated through a multi-step mechanism. The photosensory receptor pathway had the greatest number of genes, suggesting its paramount role in regulating the SAS response. For these key genes, their phenotype variation explained (PVE) was calculated, ranging from 0.99% to 1.91%. and 'defense response to bacterium, incompatible interaction' were found. All these have pointed out the importance of leaf development, immune system, far-red light regulation to the SAS, and have also promoted the screening of key genes for the further dissection of SAS.

ShapeQTLs: How Leaf Area Scales with Leaf Length
Supporting the general understanding that the elevation of plant density will cause the alteration of light in quality and quantity, this will affect the plant growth or elongation synthetically. Of these, 280, 82, 86, 90, and 7 were distributed between chromosomes 1 to 5 ( Figure 4A). All were annotated using the TAIR website (www.arabidopsis.org, accessed on 1 October 2022) according to their genomic position, resulting in 30 annotated functional proteins. In terms of protein function, there were fourteen photosensory receptor genes, four phytohormone genes, six immune pathway genes, and one leaf shape-related gene (Table 1), which implies that SAS is regulated through a multi-step mechanism. The photosensory receptor pathway had the greatest number of genes, suggesting its paramount role in regulating the SAS response. For these key genes, their phenotype variation explained (PVE) was calculated, ranging from 0.99% to 1.91%. Next, we focused on an SNP located at AT4G01895.1, which encodes systemic acquired resistance (SAR) regulatory proteins (NIMIN-1-related) belonging to the immune regulation pathway. As shown in Figure 5A, each genotype at this shapeQTL showed a Next, we focused on an SNP located at AT4G01895.1, which encodes systemic acquired resistance (SAR) regulatory proteins (NIMIN-1-related) belonging to the immune regulation pathway. As shown in Figure 5A, each genotype at this shapeQTL showed a difference in the pattern of how leaf area and leaf length change with time. With the two-dimensional expansion of the leaf, A showed a period of linear growth from 2 to 4 weeks and increased substantially, compared to the increase in L. Indeed, the curves for A showed greater differences between density conditions than between genotypes. Over time, the allometric difference between the two densities narrowed ( Figure 5B). Regarding genotypes, CC and GG showed significant differences between the two densities. In addition, all 10 shapeQTLs showed a weakening impact on the L-A relationship during the middle-late stages of growth after a slight increase in the early stage ( Figure 6A). Once a leaf reaches 4 cm in length, the additive effects decrease. In comparing effects between the two densities, the inversive patterns of 10 additive effect curves were disclosed ( Figure 6B). showed greater differences between density conditions than between genotypes. Over time, the allometric difference between the two densities narrowed ( Figure 5B). Regarding genotypes, CC and GG showed significant differences between the two densities. In addition, all 10 shapeQTLs showed a weakening impact on the L-A relationship during the middle-late stages of growth after a slight increase in the early stage ( Figure 6A). Once a leaf reaches 4 cm in length, the additive effects decrease. In comparing effects between the two densities, the inversive patterns of 10 additive effect curves were disclosed ( Figure  6B).

PositionQTLs: How Leaf Length Scales with Petiole Length
There were 418 significant SNPs identified for P-L allometry, involving 280, 30, 17, 59, and 32 SNPs on chromosomes 1 to 5 ( Figure 4B). Functional annotation identified 24 protein genes classified into the categories listed in Table 2. The pathways and gene categories were the same as those for shapeQTLs, clearly indicating their shared function in regulating SAS. Therefore, positionQTLs not only regulate leaf position by changing L or P, but also affect traits that are regulated by shapeQTLs. These key genes owned PVEs ranging from 1.01% to 1.88%. Next, we focused on a QTL located at AT3G22170.1, which encodes far-red elongated hypocotyls 3 protein. Under both densities, P changed with time; although its growth rate was slower than that of L, both P and L changed more at the density level than at the genotype level ( Figure 5C). High-density conditions led to a larger P than low-density conditions, which suggests that this QTL plays a role in developing a longer petiole to move leaves towards well-lit microsites. For L, the curves diverged between genotypes after 1.5 weeks under low density. Using the estimated parameters from Equation (1), we analyzed the allometry pattern at the genotype level. The positionQTLs fit the allometry scaling relationship for P-L well ( Figure 5D). The allometric patterns among QTL genotypes under high density were indistinguishable. This suggests a G × E interaction property with this QTL. For genotypes AA and TT, the differences in allometry scaling between the two densities were significant. Calculations of the additive effects under low density indicated that QTLs 4, 5, and 7 showed a downward parabola-shaped additive effect. QTLs 3, 6, 9, and 10 showed a sharp decrease with an increase in P ( Figure 6C). Under high density, QTLs 1, 2, 3, 6, 9, and 10 showed a decreasing role in controlling allometry in relation to L. QTLs 5,7,8,and 9 showed an increasing role with P ( Figure 6D).

Pleiotropism among QTLs
Interestingly, 89 significant SNPs were identified for shapeQTLs that were also significant SNPs for positionQTLs, which suggests that these SNPs simultaneously govern allometric scaling of both L-A and P-L. There were 55, 6, 5, 21, and 2 shared SNPs located on chromosomes 1 to 5. As the two allometric relationships have different biological meanings, these QTLs appear to be pleiotropic. AT1G10240.1, AT1G69010.1, and AT5G38860.1 were directly associated with the SAS response. These SNPs encode FAR1-related sequence 11, BES1-interacting Myc-like protein 2, and BES1-interacting Myc-like protein, respec-tively, and were all involved in photosensory pathways. The QTL at AT1G10240.1 had a significant role in regulating L-A ( Figure 7A,B) and P-L allometry under both densities ( Figure 7C,D). Such a QTL may behave more sensitively to alter the phenotype of the leaf when density changes, so that plants can produce a flexible mechanism to better adapt to shade.
icant SNPs for positionQTLs, which suggests that these SNPs simultaneously govern allometric scaling of both L-A and P-L. There were 55, 6, 5, 21, and 2 shared SNPs located on chromosomes 1 to 5. As the two allometric relationships have different biological meanings, these QTLs appear to be pleiotropic. AT1G10240.1, AT1G69010.1, and AT5G38860.1 were directly associated with the SAS response. These SNPs encode FAR1related sequence 11, BES1-interacting Myc-like protein 2, and BES1-interacting Myc-like protein, respectively, and were all involved in photosensory pathways. The QTL at AT1G10240.1 had a significant role in regulating L-A ( Figure 7A, 7B) and P-L allometry under both densities ( Figure 7C, 7D). Such a QTL may behave more sensitively to alter the phenotype of the leaf when density changes, so that plants can produce a flexible mechanism to better adapt to shade.

Discussion
SAS triggers phenotypic changes in plants, as they attempt to avoid shade and seek enhanced light interception. For plants, shade varies both in the vertical and horizontal directions and also temporally, and shade avoidance can involve any one or more of these axes. All these hold important implications for agriculture practice. For example, in the context of decreasingly available arable land and the increasing population, rational close cultivation is one of the effective strategies to increase crop yield per unit area. To date, although studies have identified dominant components of SAS regulation, the dissection on genetic control of SAS plasticity had largely not factored the co-varying relationships among distinct phenotypes of leaf.
In this study, we focused on the allometric relationships among leaf traits. Specific traits were investigated longitudinally and nondestructively by modeling the changes observed in laboratory experiments. We also modeled the allometric relationships among traits and determined how QTLs may control allometry. The association analysis for each SNP have been used to determine the genetic architecture of complex traits [17,18]. Functional mapping, a QTL identification framework that utilizes mathematical curves to model biological growth, was conceptualized two decades ago [13]. Although only one mapping population used in this study, we hope to construct a second mapping population to validate the QTL results to improve the reliability. Compared to GWAS using static traits, the use of dynamic phenotypes can provide deeper biological insights by correlating genotypes with phenotypes. Thus, since this method was introduced, it has been extensively used, for example, to dissect gravitropism and phototropism based on dynamic phenotypes of seedlings [19], to study leaf allometry between leaf area and leaf mass in the common bean [20], and to model phenotypic plasticity using growth trajectories of rice [21].
In this study, L-A and P-L relationships were used as a predictor of the capacity of Arabidopsis to respond to shade. Through functional mapping, the genetic architecture underlying allometry was dissected. Focusing on traits of leaf development and incorporating phenotypes from plants grown under two different densities, we were able to map shapeQTLs and positionQTLs. We also identified pleiotropic genes. Among genes modulating the change in A as a function of L, we identified 545 functional genes in the shapeQTLs group, 30 of which encode a molecular function: these included photosensory receptor genes, phytohormone genes, and genes regulating leaf development and leaf shape. Most photosensory receptors were phytochromes, which respond to low R:FR ratios caused by dense vegetation [22]. AT3G63300.1, encoding FORKED 1 (FKD1), plays a role in patterning leaf vascular and can coordinate leaf size with vein density [23]. Interestingly, AT1G75540.1, which encodes salt tolerance homolog2 (STH2), was also associated with shade avoidance. It interacts with COP1 to control de-etiolation and positively regulates photosynthesis [24]. Interestingly, high-density growth conditions lead to more severe plant diseases [25] which may be due to the change in light environment (e.g., a decrease in UV-B radiation typically decreases plant resistance to biotic stressors).
Genes modulating changes in L as a function of P (i.e., positionQTLs) fell into the same pathway categories as those for shapeQTLs but the gene sets were different. This implies that changes in leaf position and leaf shape are linked to the SAS response. The overlap of gene categories between L-A and P-L may reflect their combined role in driving changes in leaf position and shape. On the other hand, this overlapping further illustrated the need for QTL mapping of SAS using distinct allometric trait pairs; indeed, doing so identified another phenotype regulating SAS under high-density conditions, namely, the A trait. This indicates that there are QTL-density interactions. Remarkably, a previous study reported the role of AT3G22170 in mediating the elongation of hypocotyls in response to FR light in concert with FAR1 [26,27]. However, in our study, this QTL was related with P-L, revealing another role during later development.
The allometry underlying complex traits is pervasive. For example, the metabolic theory of ecology explaining metabolic rate scales with the three-fourths power of body mass [28]. Another example is the quarter-power scaling relationship between separate organs and overall body size [29]. As a practical case, in forms of height-diameter allometry [15] identified genes modulating the tree stem growth. Using Reeve and Huxley's allometric equation [30], the increase in the dry weight of soybean seed was examined using a functional mapping method [31]. Therefore, the allometry principle deserves a widespread application to the dissection of more complex traits. In this study, we used functional mapping to map allometry-related QTLs, with an emphasis on SAS. We utilized an RIL population, which can give sufficient replicate samples, as each progeny within the RIL population can propagate many replicates through self-reproduction. This can reduce spurious associations when mapping QTLs. In addition, we reduce other factors by growing plants in chambers with controllable conditions. All these steps improved the precision of QTL detection.
For rosette species such as Arabidopsis, vegetative internode elongation is highly suppressed and we, therefore, used leaf traits to analyze the genetic architecture of SAS. However, for other non-rosette plants, accelerated elongation of internodes and the resulting increase in stem length have been reported as the most striking SAS component [32]. This might also be accompanied by enhanced apical dominance and stimulated vertical growth of the main stem. In such species, it may be better to adjust the model to consider the growth relationship between internode length and stem height. The further modeling of such potential allometry could reveal molecular mechanisms underlying dwarf or semidwarf phenotypes, which could aid the development of modern crop species that can grow at high densities.

Plant Materials
Arabidopsis ecotype Ler (Landsbergerecta) and Sha (Shakdara) were used as parental material. F1 seeds were obtained by manual crossing and subsequent self-fertilization for 10 generations to develop a RIL population. In all, 84 progenies were used, with each progeny, 20 replicates were taken for both the high-and low-density conditions. The genome DNA for the whole RIL population was extracted with a TIANGEN DP305 DNA extraction kit, using 100 mg of young leaf. DNA quality was inspected by agarose electrophoresis and imaged with an ultraviolet imager. DNA samples without degradation were sequenced subsequently. The RIL population was genotyped at Total Genomics Solution (TGS) Institute, China. Using the Illumina Hiseq 2000 platform, the 125 bp paired-end reads were produced at a depth of 9.91×. After filtering the low-quality sequence, 1.10 GB data per sample were generated on average. The percentage of Q30 and Q20 bases were more than 90% and 95%, respectively, with a normal GC distribution. Taking the reference genome from the TAIR website (https://www.arabidopsis.org/download/index-auto.jsp? dir=%2Fdownload_files%2FSequences%2FAssemblies, accessed on 1 October 2022), SNP detection was performed with GATK, a widely used mutation detection software. The SNP screening criteria were customized as follows: the sequencing depth ≥ 4, the SNP base quality ≥ 20, the variation detection quality ≥ 50. SNPs with the allele frequency ≥ 5%, or with the missing rate ≥ 50% were filtered out. With these criteria, a total of 417,495 high-quality SNPs were determined for QTL mapping [33].

Experimental Design
From the 84 progenies, 40 seeds were sowed, from which 20 seedlings were used for high-and low-density planting (with 3 × 3 cm and 5 × 5 cm row spacing for the high and low conditions, respectively). For each density condition, a randomized block design of 20 replicates with each progeny, 5 blocks each with 4 replicates, was used. Plants were grown in an artificial climate chamber at Beijing Forestry University, China. Light and temperature conditions in the chamber were 16 h light and 8 h dark per day, with a light intensity of 150 µmol·s -1 ·m -2 and a constant room temperature of 22 • C with 80% relative humidity. A mixture of turfy soil, perlite, and vermiculite at a 1:1:1 ratio was used for cultivation.

Phenotype Collection
One week after transplanting, the third true leaf of each plant was chosen for phenotypic investigation. Selected leaves were photographed weekly using a digital camera (NIKON COOLPIX A10). To assess phenotypes, an image calibration card (a printed checkerboard, each square having a size of 10 × 10 mm) was used as illustrated in Figure 1. At the time of taking the photograph, the leaf was flattened against the blank space of the card, making the checkerboard and leaf integrity show on photo ( Figure 1A,B). The information was imported into the CameraCalibrator module of MATLAB for further analysis. Then, following a previously described method [16], all leaf shapes were reconstructed. Furthermore, the three traits under study (P, L, and A) were assessed via MATLAB.

Statistical Models
Mean phenotypes were calculated for each progeny in the RIL population, and each of the three traits were calculated. Equation (1) was used to model allometric growth: where x i is the longitudinal phenotype of one trait for progeny i, µ j is the mean vector of another trait for SNP genotype j, and α j and β j are parameters for the allometric law.
To incorporate the dynamic phenotype from the two densities into one model, taking the scaling relationship of L-A as an example, let x i = (x ih ; x il ) = (x ih (t 1 ),(x ih (t 2 ),..., (x ih (t Th ); x il (t 1 ),(x il (t 2 ),..., (x il (t Tl )) denote the joint phenotype of L for progeny i, where h represents the high-density environment and l represents the low-density environment. Similarly, y i = (y ih ; y il ) = (y ih (t 1 ),(y ih (t 2 ),..., (y ih (t Th ); y il (t 1 ),(y il (t 2 ),..., (y il (t Tl )) expresses the joint phenotype of A. For the two-density condition, we allowed the starting measurement time to vary depending on the growth after seedling transplantation.
Next, assuming two genotypes for each SNP, n progenies can be grouped into two sets. The likelihood for the dynamic trait from the two-density condition is formulated as follows: According to this equation, at a specific SNP, the given progeny i carrying QTL genotype j, the function f j|i (y ih ; y il ) µ j , Σ i is a multivariate normal distribution, where the expected mean vector µ j is expressed as µ j = µ jh ; µ jl ) = µ jh (t 1 ) , µ jh (t 2 ), .., µ jh (t Th ); µ jl (t 1 ), µ jl (t 2 ) , µ jl (t Tl For (co-)variance matrix ∑ i , because x i in Equation (1) denotes one trait within the allometric trait pair, x i has progeny-specific values. To synthesize the two density conditions, ∑ i can be expressed as where ∑ i(hh) and ∑ i(ll) represent a (T h × T h ) and (T l × T l ) residual covariance matrix of growth over the y trait, respectively, and ∑ i(hl) = (∑ i(lh) )' is a (T h × T l ) covariance matrix between the two density conditions. Based on a previous derivation of multi-trait functional mapping [34], combined with a one-order structured ante-dependence model (SAD (1)) [35] we further rewrite (4) as where I · · · 0 . . . . . . . . . . . . . . .
A maximum likelihood estimation was used to estimate curve parameters (α j ,β j ) = (α jh ,α jl ,β jh ,β jl ) for the mean vector, and (co-)variance parameters (ϕ 1 , ϕ 2 , σ 2 ). Under likelihood (2), whether a QTL for SAS regulation exists can be tested by formulating the null and alternative hypotheses as follows: H 0 : α j = α, β j = β H 1 : At least one of the above expression does not hold (7) For this hypothesis, we used the likelihood ratio test to determine QTLs, the threshold for which can be determined through 1000 replicates of a permutation test, using a significance level of 0.01. LR = -2(logL 0 -logL 1 ) We can also test how QTLs may change allometric scaling relationships under different densities to impact the SAS response. This can give deeper insights into the genetic control of leaf traits to optimize photosynthetic properties. The equation can be expressed as follows: H 0 : α jh = α jl , β jh = β jl H 1 : At least one of the above expressions does not hold (9) We performed the same calculations for P and L to determine their allometric relationship. Then, considering the three correlated leaf traits (P, A, and L), the two allometric scaling relationships, and the two (high/low) density conditions, the QTLs responsible for the SAS response were comprehensively assessed.

Functional Annotation
Using the models above, all SNPs with a high LOD value larger than the threshold were taken as QTLs. All SNPs were aligned against the Arabidopsis genome; thus, the annotation result with all SNPs can be determined, including their genomic position, and their reported function if SNPs were located at functional proteins.

Conclusions
Using the developmental trait of leaf, this study performed the QTL mapping of plant SAS with the Arabidopsis RIL population. With the aid of bivariate functional mapping, leaf traits investigated under a high-and low-density planting environment were mathematically modeled from the viewpoint of allometry. For shapeQTL and positionQTL results using the L-A and P-L allometry, a total of 30 and 24 functional SNPs were determined to mediate the plant SAS, respectively. Through functional annotation, immune pathway genes, photosensory receptor genes, and phytohormone genes were identified to be involved in the SAS response, among which, systemic acquired resistance (SAR) regulatory proteins (NIMIN-1-related) and salt tolerance homolog2 (STH2) were identified. By dissecting and comparing QTL effects, the results elucidate the genetic control of leaf formation in the context of the SAS response.