A natural population consisting of 264 Chinese soybean accessions, including 212 improved varieties and 52 landraces, was utilized in this study. Genome-wide association study of the landrace panel and the cultivated panel was conducted with 2,597,425 SNPs. The particular information has been presented in our previous research (Zhang et al., 2021).

The shading stress was simulated using shade nets that reduced light by 30%, and the results were compared to normal conditions with natural light. The study took place in Nantong (32°1’N, 120°52’E), Jiangsu Province, China. Soybean germplasms were planted in June and harvested in October of both 2022 and 2023. Each soybean germplasm material was grown in 3 rows, with 10 holes per row, and each row is filled with 20-25 plants. After harvesting, six traits (first pod height, plant height, pod number per plant, grain weight per plant, branch number and main stem node number) were measured based on the Descriptors and Date Standard for Soybean (Glycine spp.) (Qiu and Chang, 2006). The shade tolerance coefficient (STC) for each trait was used as an evaluation indicator, calculated using the following formula:

In which, y ¯ i j ( T r e a t ) and y ¯ i j ( C K ) represent average observed value of genotype i (i=1, 2, 3…264) on the trait j (j=1, 2, 3) with or without shade treatment (Li et al., 2014).

Standardize the STC of each genotype for each trait using the subordinate function value (SFV) (scaled to the interval ([0,1]) by the following formulas:

In which, min (STCij) and max (STCij) represent the j (j=1, 2, 3) trait minimum and maximum of genotype i (i=1, 2, 3… 264), respectively; W_j represents the importance or weight of the j trait among all composite indicators, where pj is the contribution rate of the j trait for each soybean genotype. The D value is the comprehensive evaluation score of shade tolerance for each soybean genotype under shade stress conditions, obtained by assessing the comprehensive indicators. Here, Xj represents the j trait (Li et al., 2014; Du, 2023).

Based on the calculated Average Subordinate Function Value (ASFV), the data for all genotypes under the current trait are evenly divided into five categories. The grouping criteria for shade tolerance in these five categories are determined, with each genotype being classified according to its shade tolerance level. The higher the ASFV, the stronger the shade tolerance of the genotype.

The population resequencing data utilized in this study was previously reported in our earlier research. In brief, high-density map includes 2,597,425 single nucleotide polymorphisms (SNPs) from the landrace and cultivated accessions, with a linkage disequilibrium (LD) decay range of 120 kb (Zhang et al., 2021). For each year of the study, ten plants were selected for measurement. Genome-Wide Association Studies (GWAS) were conducted using the GAPIT package based on R software and a mix linear model (MLM) were employed.

Genotyping was performed using three sets of primers (F1, F2, and R) specifically designed for KASP markers, as detailed in Supplementary Table S3 . These primers were designed using the Primer-Blast tool available on the NCBI website (https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_ LOC=BlastHome). Genomic DNA was extracted using the 2×CTAB method (Jia et al., 2024). PCR amplification was carried out using the KASP V4.0 2×Mastermix (JasonGen, China), following the reagent’s instructions. The amplified DNA was then analyzed using a Quantitative Real-Time PCR System (ABI Quant Studio 5).

The software IBM SPSS 20 was used for descriptive statistics and analysis of variance (ANOVA) (IBM, Armonk, NY, USA). Correlation analyses were performed using Origin software (Origin Lab, USA). The frequency distributions of six traits for the soybean accessions in both years were calculated by Microsoft Excel 2016.

To evaluate shade tolerance, we cultivated a natural soybean population of 264 accessions under conditions with or without 30% shade treatment in Nantong during 2022 and 2023. Six agronomic traits (first pod height, plant height, main stem node number, pod number per plant, grain weight per plant, and branch number) were measured across the 264 accessions over two years. Overall, the average first pod height and plant height in 2022 (E1) and 2023 (E2) under shade treatment were higher than those under normal light conditions. In contrast, the average main stem node number, pod number per plant, grain weight per plant, and branch number showed varying degrees of decline ( Supplementary Table S1 ). To evaluate the shade tolerance of the soybean population, STC for six traits was calculated, and descriptive statistical analysis were performed to the 264 accessions from 2022 and 2023 ( Table 1 ). The STC values for six traits were defined as STC1-6, respectively. As shown, the mean STC values for these six traits in 2022 and 2023 did not exhibit significant differences ( Table 1 ), with heritability (h²) values of 43.64%, 40.15%, 44.57%, 36.52%, 40.61% and 38.99%, respectively ( Table 1 ).

An analysis of variance (ANOVA) was conducted on the six traits across the 264 accessions in 2022 and 2023, revealing significant differences among genotypes, stress treatments, and different environments ( Table 2 ). To explore the correlation among the six traits in 2022 and 2023 for the soybean population, a correlation analysis was conducted ( Figure 1 ). The results indicated that the STC1 and STC2, STC2 and STC3, STC3 and STC4, STC4 and STC6 showed significant positive correlation in 2022 and 2023 ( Figure 1 ). While, STC1 and STC4 in 2022, STC1 and STC4, STC5 in 2022, STC2 and STC5 in 2023 exhibited significant negative correlation ( Figure 1 ). Between two years, only STC4 in 2022 and STC2 in 2023, STC1 in 2022 and STC5 in 2023 shown negative correlation, STC5 in 2022 and STC4 in 2023, STC4 in 2022 and STC5 in 2023 shown significant positive correlation, respectively. But other traits between two years shown weak significant correlation ( Figure 1 ).

In 2022, the ASFV thresholds for different levels of shade tolerance were as follows: high shade tolerance was above 0.68, shade tolerance ranged from 0.47 to 0.68, moderate shade tolerance ranged from 0.38 to 0.47, shade sensitivity ranged from 0.29 to 0.38, and high shade sensitivity was below 0.29 ( Figure 2A ). In 2023, the thresholds were slightly adjusted: high shade tolerance was above 0.70, shade tolerance ranged from 0.56 to 0.70, moderate shade tolerance ranged from 0.49 to 0.56, shade sensitivity ranged from 0.41 to 0.49, and high shade sensitivity was below 0.41 ( Figure 2B ). Over the two years, moderate shade tolerant soybean germplasm was the most prevalent, comprising approximately 39% and 41% of the total population. Shade sensitive germplasm ranked second after moderate shade tolerant germplasm. High shade tolerant germplasm was relatively rare, accounting for 0.76% and 1.52% of the total soybean population in 2022 and 2023, respectively ( Figure 2A ).

In summary, this study identified a total of five high shade tolerant soybean germplasms over the two years, with NPS044 being selected in both years ( Table 3 ). These high shade tolerant materials offer a valuable foundation for further research into the genetic mechanisms underlying soybean shade tolerance and serve as important experimental materials for future breeding programs aimed at enhancing shade tolerance in soybeans.

To pinpoint key genomic loci responsible for shade tolerance in soybeans, we conducted GWAS on six traits across 264 soybean accessions under control and shade conditions, as well as STC of the six traits for the years 2022 and 2023 ( Figures 3 – 6 ; Supplementary Figures S3 , S4 ). The resulting frequency distribution maps and density curves indicated that the phenotypic data for the six traits followed a continuous distribution. This suggests that the natural soybean population in our study harbors rich genetic variation, making it well-suited for further GWAS analyses ( Supplementary Figures S1 , S2 ).

Over the course of two years, we identified a total of 733 significant SNPs associated with STCs of six traits ( Table 4 ). Specifically, for the STC of first pod height, 28 SNPs were detected. In 2022, 24 SNPs were associated with the STC of first pod height, while in 2023, 4 SNPs were identified. Notably, S11_19943066 was significant under control conditions, whereas S11_19738980 was significant under shade, with these two SNPs being approximately 204 kb apart ( Figures 3A, B ; Table 4 ). In 2023, S15_50935714 was significant under shade treatment, while S15_51517560 showed a significant correlation with the STC, with a distance of approximately 582 kb between these two SNPs ( Figures 3E, F ; Table 4 ).

In 2022, 29 SNPs were linked to the STC of plant height, while in 2023, 9 SNPs were identified. The SNP S15_15344441 showed significant association under control conditions, and S15_19750516 was significantly correlated with the STC ( Figures 4A, B ; Table 4 ). Moreover, S19_45102497 and S19_45149787 were significantly associated with plant height under control and shade conditions, respectively, with a distance of approximately 47 kb between them. The SNP S14_4288867 was significantly correlated with both plant height and STC. Additionally, S15_36102679 and S15_36134614 were significantly associated with soybean plant height under shading conditions and STC, with these two SNPs being approximately 32 kb apart ( Figures 4D–F ; Table 4 ).

For the main stem node number, two SNPs showed significant correlation with the STC. The SNPs S04_11807969 and S14_4288867 were notably correlated with both the main stem node number under shading and the STC. Additionally, 2 SNPs were significantly correlated with the STC in 2023 ( Supplementary Figure S3 ; Table 4 ).

For the trait of pod number per plant, 13 SNPs were identified as significantly correlated with the STC of pod number per plant in 2022. Additionally, 31 SNPs were significantly associated with the STC ( Supplementary Figure S4 ; Table 4 ).

Regarding grain weight per plant, a total of 105 SNPs in 2022 and 190 SNPs in 2023 were significantly correlated with the STC. Specifically, the SNP S15_1302099 showed a significant correlation with both single plant grain weight and the STC of soybean under shading. The SNP S10_43122500 was significantly associated with grain weight per plant under shading conditions, while S10_43121286 was significantly correlated with the STC, with these two SNPs being approximately 1 kb apart. Notably, the SNP associated with the STC for grain weight per plant had the highest explanatory power, with a -log₁₀(p) value of 9.28 and a phenotype explanatory rate of 18.02% ( Figure 5 ; Table 4 ).

A total of 324 SNPs were found to be significantly correlated with the STC of branch number. Specifically, 227 SNPs in 2022 and 97 SNPs in 2023 were significantly associated with the STC ( Table 4 ). In 2023, the SNP S18_55349193 showed significant correlation with branch number under control conditions, while S18_55354172 was significantly correlated with the STC, with these two SNPs being approximately 5 kb apart ( Figure 6 ).

To explore the phenotypic effects of allelic variations in significant SNPs, a haplotype analysis was conducted on the SNPs with the highest threshold detected for shade tolerance during the mature stages of 2022 and 2023. This analysis revealed a total of 4 SNPs showing significant differences between each genotype. For instance, the SNP S18_1766721 exhibited an allelic variation from A to G. The STC of first pod height was significantly higher in germplasm carrying the S18_1766721-G allele compared to those with the S18_1766721-A allele ( Figure 7A ). Another example includes the S09_48870909, which has a G/T allelic variation ( Figure 7B ). The nucleotide change at position S19_49517336 involves a substitution from G to A. Soybean with the S19_49517336-G allele exhibit a significantly higher average STC of pod number per plant compared to those with the S19_49517336-A allele ( Figure 7C ). For the S18_3429732, the allelic variation consists of A and G. Soybean germplasm with the S18-3429732-A allele has a significantly higher STC for average grain weight per plant compared to those with the S18-3429732-G allele ( Figure 7D ). The phenotypic variation explain rate of S18_1766721, S09_48870909, S19_49517336, and S18_3429732 are 8.32%, 9.01%, 8.93%, 8.73%, respectively. The favorable alleles ratio in the population of S18_1766721-G, S09_48870909-T, S19_49517336-G and S18_3429732-G were 37.9%, 34.5%, 87.5%, and 13.6%, respectively ( Table 5 ). And the HST germplasms NPS044, NPS060, NPS151, NPS187 and NPS254 each contained three, three, two, four and one favorable alleles ( Table 5 ). Also, we developed KASP markers for these four SNPs ( Table 6 ). The designed molecular markers effectively differentiate between these two genotypes ( Figures 7E–H ).

The LD of this population is 120 kb (Zhang et al., 2021). Therefore, we examined candidate genes within a 120 kb range upstream and downstream of SNPs significantly associated with soybean STC across six traits. Utilizing functional annotation information from the soybean genome, we identified four candidate genes significantly linked to soybean shade tolerance ( Table 7 ). Glyma.18G024000 associated with S18_1766721, encodes a trichome birefringence-like 33 protein. The gene Glyma.09G271100 linked to S09_48870909, encodes a protein from the auxin efflux carrier family. Glyma.19G248900, associated with S19_49517336, encodes an ethylene response factor 1. Lastly, the gene Glyma.18G040700, related to S18_3429732, encodes a MYB domain protein 43.

In natural environments, plants are often subjected to shade tolerance. Shade tolerance is essential for soybeans, especially in intercropping or relay cropping systems. When soybeans experience shade stress, their plant height and first pod height will be elongated and the branches, number, grain weight per plant, pod number and nodes number will be reduced, which posed a huge threat to soybean production (Yang et al., 2015; Raza et al., 2020; Table 1 ).

Previous studies have demonstrated that a 15% reduction in light is considered weak shading, whereas 60% shading often leads to lodging in most varieties, indicating excessive shading. However, at 30% shading, the proportion of lodging varieties and the coefficient of phenotypic variation are sufficient to meet the requirements for shade tolerance identification (Sun et al., 2017; Zhang, 2021). Therefore, this study employed a 30% light reduction to simulate shade treatment. As a result, under 30% shading, all the six traits exhibited more pronounced phenotypic changes, and the result all present normal distribution ( Supplementary Table S1 ; Supplementary Figures S1 , S2 ).

ASFV has been widely utilized for assessing crop resistance to various stressors, including salt, drought and shade tolerance (Zhao et al., 2023). In this study, ASFV was applied to evaluate the shade tolerance of soybean ( Figure 2 ). Five germplasms exhibiting high shade tolerance were identified: NPS044, NPS060, NPS151, NPS187 and NPS254. Notably, NPS044 showed consistent results across two years ( Table 3 ). Chen et al. (2003) measured various parameters such as STC of biological yield during pod setting, plant height, minimum pod height, pod number per plant, grain number per plant, grain weight per plant, and 100 grain weight, and calculated the ASFV of soybean varieties. Similarly, Huang et al. (2012) employed a comprehensive STC across nine indicators, including standard pod number, standard pod weight, 100 grain weight, plot yield, plant height, main stem node number, number of ingle grain pods number per plant, single plant pod weight per plant, and standard pod length, to determine soybean shade tolerance. Li et al. (2014) developed a mathematical model for evaluating soybean shade tolerance using stepwise regression and identified seven key indicators: main stem node number, branch number, internode length, lodging resistance, pod number per plant, 100-grain weight, and grain weight per plant. Wu et al. (2015) suggest that the rapid identification and prediction of shade tolerance in soybean seedlings can be achieved by measuring leaf dry weight, stomatal conductance, plant height, and maximum fluorescence yield under dark conditions. Tang et al. (2022) used traits such as first pod height, stem node number, pod number per plant, and grain number per pod to evaluate shade tolerance. In this study, six traits - first pod height, plant height, pod number per plant, grain weight per plant, branch number and main stem node number – were measured using STC as the indicator to evaluate shade tolerance in 264 soybean accessions. These traits are reliable for identifying key loci and genes associated with shade tolerance in soybeans.

A total of 733 SNPs were identified as being associated with the STC of six traits over two years ( Figures 3 – 6 ; Supplementary Figures S3 , S4 ; Table 4 ). Due to the significant influence of environmental factors on these traits ( Supplementary Figure S1 ), we didn’t co-locate any significant loci between two years. More environments may need to be added.

Based on GWAS, four SNPs S18_1766721, S09_48870909, S19_49517336 and S18_3429732 were selected for further study. Their phenotypic explanation rate ranged from8.32% to 9.01%, which can be used for soybean genome selection breeding. Four candidate genes associate with the four SNPs were identified. Glyma.18G024000, associated with S18_1766721, encodes the protein Trichome birefringence-like 33 (TBL33). Members of the TBL family, previously characterized, are localized in the Golgi apparatus and function as polysaccharide O-acetyltransferases catalyzing the O-acetylation of specific cell wall polymers (Stranne et al., 2018; Sinclair et al., 2020; Lunin et al., 2020). TBL proteins have been reported to play roles in biotic (disease, herbivore) and abiotic resistance (salt, drought and freezing) (Xiong et al., 2013; Gao et al., 2017; Sun et al., 2020). Glyma.09G271100 is an auxin efflux carrier family protein, known as PIN -like (PILS) which plays a crucial role in auxin signaling (Bogaert et al., 2022; Feraru et al., 2022; Waidmann et al., 2023). Numerous studies have demonstrated that auxin plays pivotal roles in integrating responses to abiotic stresses such as temperature, water, light and salt and in controlling downstream stress responses (Iglesias et al., 2018; Waadt et al., 2022; Xie et al., 2022; Jing et al., 2023). Organ-specific transcriptome analysis has revealed that shade induces a set of auxin-responsive genes, such as SMALL AUXIN UPREGULATED RNAs (SAURs) and AUXIN/INDOLE-3-ACETIC ACIDs (AUX/IAAs) (Nguyen et al., 2023). In the initial response to shade signals, auxin biosynthesis, transport, and sensitivity are rapidly activated, promoting cell elongation in hypocotyls and other organs (Ma and Li, 2019). Glyma.19G248900 associated with S19_49517336, encodes an ethylene response factor 1 (ERF1). GmERF3 has been reported to positively regulates resistance to virus, high salinity and dehydration stresses (Zhang et al., 2009; Liu et al., 2024). Ethylene is known to play a crucial role in mediating plant adaptations to environmental conditions (Waadt et al., 2022). Recent studies have shown that shade stress can induce ethylene biosynthesis, accelerating soybean senescence and hindering nitrogen remobilization (Deng et al., 2024). ERFs are significant in enhancing flood tolerance in rice (Xu et al., 2006), where ethylene accumulation in submerged tissues induces the expression of ERFs such as SNORKEL1 and SNORKEL2, which are major QTLs associated with deepwater internode elongation (Hattori et al., 2009). Glyma.18G040700 related to S18_3429732, encodes MYB domain protein 43. This protein plays a critical role in various aspects of plant growth and development, including secondary metabolic regulation, responses to hormones and environmental factors, cell differentiation, organ morphogenesis, and cell cycle regulation (Li et al., 2019a). The homolog AtMYB43 has been reported to be involved in regulating tolerance to cadmium and freezing (Zheng et al., 2022, 2023). In summary, these four genes Glyma.18G024000, Glyma.09G271100, Glyma.19G248900, and Glyma.18G040700 may be involved in soybean responses to shade tolerance.

This study developed four KASP markers based on SNPs associated with soybean STC obtained from GWAS. These markers have been successfully used for genotyping ( Table 7 ; Figure 7 ). Specifically, S18_1766721 is associated with the STC of first pod height, S09_48870909 with the STC of plant height, S19_49517336 with the STC of pod number per plant, and S18_3429732 with the STC of branch number. These markers are valuable tools for identifying shade-tolerant soybean germplasms and can enhance the efficiency and accuracy of selection in molecular marker-assisted breeding. However, KASP markers for the STC of node number and grain weight were not developed. This gap may be due to the influence of multiple factors, suggesting that further efforts and research are needed to identify effective markers for these traits.