Construction and fitting of the HDEE
The selection of the optimal expression of autogenic regulation terms \({\phi }_{i}\left(l\right)\) is based on the numerical experimental results of phenotype fitting, which include standard error, goodness of fit, and the comparative implementation of three information criteria (S Table 1, Additional file 2). In this study, the dataset of 98 samples was divided into 6 groups using Admixture software, as depicted in Fig. 2A, and the difference curves of photosynthetic phenotypes among the 6 groups with varying light intensities are presented in Fig. 2B. The trajectories of phenotypes with varying light intensities reveal a consistent average dynamic trend among different groups. The differences are primarily observed in the variance of the group phenotypes. Figures 2C-E show the distribution of group disturbance parameter values, providing a more intuitive observation that despite variations in the overall distribution range of group disturbance parameters, the average parameters among different groups are highly similar. The results of variance analysis for the parameters of the three phenotypes across different groups indicated nonsignificant population differences, as evidenced by p values much higher than 0.05 in S Table 2 [see Additional file 2]. Furthermore, we conducted an analysis of the genetic relationships among all samples (Fig. 2F) and determined that the correlation coefficient of the genetic relationships ranged from 0 to 0.48. Moreover, it was observed that more than 98.46% of the genetic relationships had correlation coefficients below 0.3. Based on these findings, it can be inferred that the impact of kinship within the natural population studied is relatively weak. In this study, the impact of population structure differences and genetic relationships among samples on phenotypic data was found to be negligible. Therefore, the disturbance term associated with population structure was omitted from the model. Combined with the results of the numerical analysis in S Table 1, HDEE can be expressed as:
$$\left\{\begin{array}{c}{y}_{1}=\frac{{K}_{1}}{1+{{a}_{1}e}^{-{b}_{1}l}}+{c}_{1}{l}^{{d}_{1}}\\ {y}_{2}={K}_{2}\left(1-{a}_{2}{e}^{-{b}_{2}l}\right)+{c}_{2}{l}^{{d}_{2}}\\ {y}_{3}={K}_{3}{e}^{-\frac{{a}_{3}}{{l}^{{b}_{3}}}}+{c}_{3}{l}^{{d}_{3}}\end{array}\right.$$
1
S Table 1 illustrates that the equations incorporating the environmental disturbance term provide a more precise characterization of the changing patterns observed in the growth and development of the three photosynthetic phenotypes compared to the classical growth equations. Figures 3A-C depict scatter plots illustrating the variation in mean values for the three phenotypes in response to changes in gradient light intensity. The overall fitted curve, obtained using the HDEE model, closely aligns with the observed mean values of the actual phenotype, indicating a strong correspondence between the model and the real data. The independent growth curves and the environmental disturbance curves elucidate the dynamics of photosynthetic phenotypes and the variation attributed to environmental disturbance during the gradient light intensity change process. As light intensity increased, the phenotypic values of ETR and qN exhibited a gradual increase, whereas the value of qP decreased. The corresponding independent growth curves consistently followed this trend. However, the trend of the environmental disturbance curves did not entirely align with that of the overall phenotype. While the environmental disturbance showed a promoting effect on the increase in phenotypic values, the environmental disturbance of ETR increased with light intensity, the environmental disturbance of qP decreased with light intensity, and the environmental disturbance of qN tended to approach zero.
In addition, Figs. 3D-F show the proportion of independent growth curves and environmental disturbance curves in the overall phenotype. It can be observed that while the independent growth of ETR exhibits an upward trend with increasing light intensity, its relative proportion in the overall dynamics gradually decreases, while the proportion of environmental disturbance dynamics gradually increases. This indicates that the influence of light intensity gradually enhances the role of "environmental disturbance" in the photosynthetic phenotype. For qP, its environmental disturbance curve shows a downward trend, but its proportion in the overall phenotype is gradually increasing, indicating that the impact of environmental interference on qP is also gradually increasing. Although the disturbance curve of qN is approximately 0, the scale diagram reveals that at low light intensities (approximately 2), the proportion of environmental interaction components is 100%, but the proportion gradually decreases to 0 when the value reaches approximately 3. This observation indicates that qN is weakly influenced by the environment under conditions of high light intensity. The intensity of environmental disturbance is clearly influenced by light intensity, and the dynamic changes in phenotypic autogenic regulation in response to gradient light intensity may be influenced by genetic regulation.
Identification of QTLs for High-dimensional Photosynthetic System
In the QTL mapping analysis, we initially conducted a genome-wide association analysis on 4,996,309 high-quality SNP data using the functional mapping approach with several classical equations (specifically focusing on the autogenic regulators in Eq. 1). This analysis led to the identification of 10,000 SNP loci as potential key QTL sets, which were considered to have the highest likelihood of influencing the photosynthetic phenotype under investigation. For the loci with high probabilities of influencing the dynamic process of photosynthetic phenotypes in response to light intensity, we employed static methods to identify sites that significantly regulate photosynthetic phenotypes across 11 light intensity groups, as depicted in S Fig. 1 [see Additional file 1]. The identified gene loci have a significant impact on the variation in photosynthetic phenotypes under varying light intensities, and the gradual increase in light intensity results in alterations in the number and pathways of associated genetic loci involved in phenotypic plasticity. To investigate how significant QTLs regulate dynamic changes in phenotypes under gradient light intensity, we performed a high-dimensional phenotype-genotype association analysis based on HDEE. This analysis allowed us to identify the top 3% of significant quantitative trait loci associated with the three phenotypes, as illustrated in Fig. 4. Based on the autoregulatory growth model, which is represented by the classical growth equation, we found that in the association analysis of ETR, qP, and qN, 127 (42.33%), 96 (32.00%), and 76 (25.33%) of the significant sites, respectively, were annotated to genes with known functions in Populus trichocarpa. By utilizing HDEE for association analysis, we identified 129 (43.00%), 101 (33.67%), and 71 (23.67%) SNPs within genes of known function, respectively.
In the comparative analysis of significant loci identified using the autoregulatory growth model and HDEE (Fig. 5), we specifically emphasize loci that exhibit significant effects in both models. This is important because it allows us to mitigate the potential presence of false positives in the QTL mapping results. The loci that exhibit significance in both models simultaneously are considered to have higher reliability. These shared loci play a crucial genetic role in the genetic variation of the corresponding photosynthetic phenotype, exerting a strong and significant effect on the associated phenotype. Furthermore, these loci were identified through screening with HDEE, incorporating environmental regulators. This suggests that the regulation of gene loci in phenotypic expression may be influenced, to some extent, by environmental perturbations in addition to autologous growth regulation. For the ETR, qP, and qN phenotypes, we identified 209, 216, and 268 significant loci, respectively, that regulate photosynthetic phenotypes in response to changes in gradient light intensity. S Fig. 2–4 displays the results of the gene ontology (GO) function enrichment analysis for the genes associated with these identified loci [see Additional file 1]. Specifically, for the ETR phenotype in response to changes in light intensity, the top three enriched functions of the regulating genes include biological processes such as nucleic acid binding and protein binding, along with a molecular function related to catalytic activity. Significant sites that regulate the dynamics of qP phenotypes exhibit important functions such as ATP binding and protein phosphorylation in biological processes, along with the molecular function of protein kinase activity. Similarly, significant gene functions at prominent sites associated with qN phenotypes also involve ATP binding, protein kinase activity, and protein phosphorylation.
Moreover, within the intersection of these pairs of significant loci, several noteworthy QTLs were found to regulate the growth and development of two phenotypes simultaneously (highlighted in the red box in Fig. 5B). Specifically, SNPs 2309065 and 3844928 influence the changes in both qP and qN in response to light intensity. Additionally, SNP 2443013 exerts an impact on both ETR and qP dynamics under varying light intensities. The pleiotropic effects of these loci are believed to be significant in the context of high-dimensional photosynthetic phenotypic systems.
QTL Network Construction for Photosynthetic Phenotypes
The overlapping QTLs between the significant loci identified using the autogenous regulatory growth model and those identified using HDEE demonstrate significant regulatory effects on the associated photosynthetic phenotypes. Furthermore, there exist interactive relationships among the effects of these QTLs on phenotypes during adaptation to varying light intensities across gradients. Utilizing the dynamic curve of genetic effects of QTLs governing phenotypic changes, we established genetic regulatory networks among QTLs based on the principles of differential evolutionary game theory and modeled how the genetic regulation of a significant QTL on phenotype is influenced by its inherent control and the strategies employed by other QTLs (Figs. 6A-C). In the context of a network node A, an outgoing link refers to a link originating from node A and pointing toward other nodes. Conversely, an incoming link of node A denotes a link directed toward node A from another node. Based on the statistical data depicted in Figs. 6D-F, it is evident that despite the potentially large number of nodes within the network structure, each node is primarily influenced by only a limited number of other nodes in the network. Specifically, the number of incoming links between nodes tends to concentrate within the range of 1–13. This observation provides empirical evidence confirming the network sparsity principle within the genetic regulatory network that governs the photosynthetic phenotype. Nevertheless, there is substantial variation in the number of outgoing links among different nodes, ranging from 0 to 107, with only a few nodes exhibiting a notably high number of outgoing links. In the global network structure, nodes characterized by a significant number of outgoing links are referred to as network hubs. In this study, the designation of hubs is attributed to nodes with outgoing link counts exceeding 20% of the total number of nodes in the network. Given the robust regulatory influence exerted by hubs on other nodes within the network, they assume a crucial role in the dynamic manifestation of photosynthetic phenotypes.
Within the ETR QTL network structure, a mere 865 links (1.99% of all possible links among 209 nodes) were identified. Among these links, 58.15% were determined to be positive, indicating their promotion of regulatory interactions, while the remaining 41.85% were classified as inhibitory links. The relationship observed among genetic loci governing ETR in response to varying light intensity primarily exhibited a cooperative nature. The network is distinguished by the presence of hubs, namely, SNP 4258266, SNP 4519045, SNP 2161006, SNP 4536279, SNP 4317096, SNP 3650470, SNP 4283036, SNP 4855429, SNP 3719550, and SNP 2691420. Based on gene function annotation results, it was determined that SNP 4519045 is positioned within the gene Potri.018G018800, which exhibits the function of 3-hydroxyisobutyryl CoA hydrolase activity (GO:0003860) in Populus trichocarpa. This specific site corresponds to the protein family Enoyl-Coa hydratase/isomerase (PF16113). Furthermore, SNP 4536279 is situated within the Potri.018 G036600, SNP 4317096 is located within the Potri.017 G017000, and SNP 4855429 is positioned within the Potri.019 G098500. Additionally, SNP 3719550 is located within the Potri.014 G118100, and SNP 2691420 is positioned within the Potri.009 G023800. For detailed information regarding the specific functions of genes in Populus trichocarpa and the homologous genes in Arabidopsis thaliana, please refer to S Table 3 [see Additional file 2].
Individuals harboring different genotypes can exhibit varying responses to environmental fluctuations, resulting in diverse phenotypic landscapes. To investigate the specific genetic structure underlying the regulatory changes in phenotypic traits, we focused on an intercross locus, namely, SNP 4536279, which was selected from the hub node within the ETR network. Based on the analysis of genotypic difference curves (Fig. 7A) for SNP 4536279, it was observed that samples with the TT genotype exhibited higher phenotypic values (asymptotic growth value) of ETR with increasing light intensity, followed by samples with the GG genotype and then the GT genotype. The independent curves showed consistent differences in line with these observations. The sample corresponding to the TT genotype exhibited a greater environmental disturbance term, while the perturbation term associated with the GT genotype was larger than that associated with the GG genotype. These findings suggest that although the sample with the GT genotype had a lower ETR value, it demonstrated a higher proportion of regulation in response to the light environment.
The disparity in phenotypic values among different genotypes represents the manifestation of the genetic effect of the QTL on the phenotype. From Fig. 7B, it is evident that the overall genetic effect value of SNP 4536279 exhibits an upward trend as the light intensity increases. This indicates that the genetic regulation of genes on the ETR response to changes in light intensity intensifies with higher light intensity. The curves depicting independent genetic effects displayed similar trends as the light intensity changed. However, when the light intensity reached 6, the value of independent genetic effects was notably lower than that of the overall genetic effects. The higher overall genetic effect of SNP 4536279 can be attributed to the promotion exerted by SNP 4258266 and SNP 4519045, specifically within the Potri.018 G018800. Interestingly, based on the observed trend in the effect curve, it can be noted that the promoting effect of SNP 4519045 gradually diminished after reaching a light intensity of 7, eventually even transitioning into a negative value, indicating an inhibitory role. Conversely, the promoting effect of SNP 4258266 demonstrated a gradual increase after the light intensity reached 7. Consequently, it is plausible to suggest that the effects of these two sites on SNP 4536279 may exhibit mutual exclusivity.
As a hub node that assumes a dominant regulatory role within the network, our primary focus lies in understanding how this particular site regulates other sites. Based on the findings depicted in Fig. 7C, SNP 4536279 exerts a regulatory effect on 65 nodes. Among these nodes, the promotion effect accounts for 55.38%, while the inhibition effect accounts for 44.62%. The intensity of regulation remained close to 0 within the light intensity range of 2–6. However, as the light intensity continued to increase, the regulatory effect significantly intensified. Notably, the most pronounced regulation occurred when the promotion effect on SNP 142656 reached 0.45 and the inhibition effect on SNP 150702 reached 0.41, both observed at a light intensity of 8.
Within the genetic regulatory network of qP, a total of 908 links were identified, which corresponds to 1.96% of all potential pairwise regulations. Among these links, approximately 66.96% were found to have a positive regulatory role, while approximately 33.04% exhibited a negative regulatory role. Cooperation emerged as the primary interaction strategy among the genetic loci involved in regulating qP changes in response to varying light intensity gradients. The hub nodes in the qP network, characterized by having the number of outgoing links surpassing 20% of the total number of nodes, encompass SNP 4363461, SNP 2637121, SNP 3075616, SNP 4889531, SNP 4389471, SNP 4289989, SNP 667824, SNP 2802797, and SNP 4535545. SNP 4289989 is situated within the candidate gene Potri.016 G130600, which exhibits GTPase activity (GO:0003924) and GTP binding (GO:0005525). Additionally, according to the Pfam annotation, it belongs to the Ras family (PF00071), a significant regulator involved in vesicle formation, motility, and fusion. We conducted an analysis of the genotypic curves (Fig. 7D) associated with the hub locus responsible for regulating the photosynthetic phenotype qP in response to varying light intensity gradients. Our findings revealed that environmental interference under different genotypes consistently exhibited a promotion effect on the qP phenotype as the light intensity gradient changed. However, it was observed that the promoting effect diminished with increasing light intensity. Additionally, it was found that the samples with genotype GA exhibited a higher degree of environmental disturbance.
Based on the analysis of the genetic effect of the hub SNP 4289989 (Fig. 7E), the effect curve exhibited an initial increase within the range of low light intensity (2–3) and subsequently decreased with the changing gradient light intensity. Moreover, the trend and variation range of the independent genetic effect curves were essentially similar to those of the overall curves. SNP 4535545 is promoted by SNP 768114 (Potri.002G095700) and inhibited by SNP 4851971. Additionally, the effect of SNP 4713423 (Potri.019G016100) exhibits a promotion trend within the range of 2–5, followed by inhibition within the range of 5–7 and subsequent promotion within the range of 7–8. In the overall network structure, SNP 4535545 exhibited regulatory effects on 35 SNPs, with the promoting effect accounting for 62.86% and the inhibiting effect accounting for 37.14% (Fig. 7F). The more significant regulatory effects of SNP 4535545 included the inhibition of SNP 2097359, which exhibited an increasing trend (reaching 0.51) within the light intensity range of 2–6, followed by a decrease. Additionally, there was promotion observed on SNP 4993953 and SNP 4851971, which initially increased and then decreased, with maximum values of 0.47 and 0.27, respectively.
In Fig. 4C, the genetic regulatory network of qN consists of 1250 links, which represents 1.75% of all possible relationships. Among these links, 770 indicate facilitation, accounting for 61.6% of the total, while 480 links exhibit inhibition, accounting for 38.4% of all existing relationships. Among them, SNP 1461140, SNP 4782328, and SNP 4505330 serve as hub nodes in the network. In this study, we conducted an analysis to investigate the genetic control of SNP 1461140 on the phenotypic trait qN. From the analysis of Fig. 7G, it is evident that environmental perturbations exhibit a mild promoting effect on the variation in the qN phenotype value under gradient light intensity. Notably, the samples with genotype GA demonstrate a stronger influence of environmental perturbations compared to the samples with genotype GG. The overall genetic effect curve exhibited an initial increase followed by a decrease, whereas the independent genetic effect curve generally showed an upward trend, with a slight decline observed between light intensities 4–6 (Fig. 7H). In the light intensity range from 2 to 6, the overall genetic effect surpassed the independent genetic effect, primarily due to the promotion effect of SNPs 2776614, 4667554, and 4594627 on the genetic control of phenotypic response to light intensity at this locus. Although SNP 1461140 plays a pivotal role in the genetic regulatory network of qP, its genetic control of the phenotype primarily arises from the promotion effect exerted by other SNPs. However, when the light intensity reached 7, the overall genetic effects became inferior to the independent genetic effects, mainly due to the suppression exerted by SNP 2145968 and SNP 4675672. We observed that although SNP 1461140 is not located within known genes, it is influenced by several loci within genes, including SNP 2776614 (Potri.010G001600), SNP 4594627 (Potri.018G060900), SNP 2145968 (Potri.006G199000), and SNP 4675672 (Potri.018G112200). The influence of these SNPs on SNP 1461140 may explain its significance as a central node in the network structure of qN. Additionally, we observed that the genetic effects of 63 other SNPs were transitively promoted, while the genetic effects of 22 SNPs were inhibited in the network (Fig. 7I). Notably, the promotion effect on SNP 212941 increased to 0.33 with the progressive rise in light intensity, which was considerably higher compared to other SNPs.
Analysis of the pleiotropic effects of QTLs on multiple phenotypes
Based on the comparative analysis results presented in Fig. 5, it can be observed that certain loci within the high-dimensional photosynthetic phenotype system exhibit simultaneous regulation of multiple photosynthetic phenotypes in response to changes in gradient light intensity. This phenomenon indicates the presence of pleiotropic control, where gene regulation influences the expression of multiple phenotypic traits. We selected a specific locus, SNP 2309065, situated on chromosome 7, which exhibits an impact on the response of both qP and qN values to changes in light intensity. The objective of our analysis was to investigate the genetic control exerted by this locus on the two photosynthetic phenotypes. As depicted in Figs. 8A and 8C, it is evident that the overall qP curves, under the influence of SNP 2309065 genotypes, exhibit a higher magnitude compared to the independent dynamic curves. This observation suggests that the perturbation of the light environment positively influences the phenotype qP, with a more pronounced effect observed in samples with the CG genotype. Regarding the photosynthetic phenotype qN, the CG sample exhibited lower phenotypic values than the CC sample, but their values tended to converge under high light conditions. Specifically, the overall dynamic curve of the CC samples closely aligned with the independent curve, while the independent growth curve of the CG samples significantly surpassed the overall dynamic curve when exposed to high-intensity light perturbation. Notably, the magnitude of the perturbation term was positively correlated with the light intensity. Furthermore, as anticipated, the impact of light intensity perturbation on the overall photochemical quenching exhibited a progressive increase with rising light levels. In contrast, the proportion of the effect on nonphotochemical quenching gradually diminished. Among them, the proportion of light environment disturbance in qP was found to be higher in the CC genotype than in the CG genotype, accounting for 22.34% versus 12.11% under low light conditions and 69.45% versus 43.07% under high light conditions. In contrast, the proportion of the disturbance term in phenotypic qN was lower, at 1.58% versus 19.74%. These results indicate that Populus simonii individuals with the CC genotype at SNP 2309065 exhibit greater potential for photosynthetic activity. On the other hand, individuals with the CG genotype demonstrate a stronger ability to implement self-photoprotection through the dissipation of excess excitability via heat dissipation.
In the network structure (Fig. 8B-D), the independent genetic effects of SNP 2309065 on the regulation of qP and qN responses to light were found to be smaller than the overall genetic effects. However, as the light intensity increased, the independent effects on regulating qP gradually decreased, and the disparity between the independent effects and the overall effects increased. Under high light intensity, SNPs 4363461, 4745041, SNP 1027933, and 3825266 exhibit a shared positive promotion effect on the genetic regulation of the qP phenotype. However, the genetic control of the locus on the phenotype is reversely inhibited by SNP 3825266. In the qN network, the independent effect gradually increased, causing the overall effect curve to converge with the increase in light intensity despite the presence of regulatory effects at other sites. The promoting effects of SNP 2445348 (Potri.008G032901), SNP 4803797, SNP 3651672, and SNP 1461140, as well as the inhibiting effects of SNP 1956825 (Potri.006G097600) and SNP 3278883 (Potri.012G044701), offset each other.
Based on the quantile analysis of growth parameters of the phenotypic curves under different genotypes (Fig. 9), SNP 2309065 exhibited genetic control over the five growth parameters of qP and qN. However, significant differences were primarily observed in the regulation of autogenic regulation parameters. In the regulation of phenotypic qP, the median of the asymptotic value (K) in the curves with increasing light intensity was smaller for the CG genotype than for the CC genotype. The difference between genotypes is also observed in the distribution of the parameter a associated with the initial value of the curve and the intrinsic growth rate b, showing similar patterns. For qN, genetic control is predominantly observed in the asymptotic value K and the parameter b, which determines the position of the inflection point of the function. Notably, the genotype CG exhibits a larger K value, while the b value is smaller. The mean values of parameter a, which is associated with the relative growth rate of trees across different genotypes, exhibited minimal variation. However, the samples from genotype CG showed a highly concentrated distribution, whereas genotype CC samples were concentrated toward the higher end of the mean value.
The average values of the disturbance rate and disturbance scale coefficient in the environmental disturbance term, which are of particular interest, are comparable between the two genotypes. However, it is worth noting that the parameter distribution of the CC genotype is more dispersed for qP. This finding suggests that Populus simonii with the CC genotype at SNP 2309065 may exhibit a higher degree of responsiveness to environmental fluctuations. In terms of the parameters of phenotype qN, there is a more pronounced difference in the parameter distribution between the two genotypes. Specifically, the parameter c shows greater variability and tends to have larger values in the CC genotype samples. Conversely, the parameter d exhibits higher variability and tends to have smaller values in the CC genotype samples. The results of the parameter analysis indicated that Populus simonii with the CG genotype at SNP 2309065 exhibited greater stability in response to changes in the light environment.