The NJRINP contains 284 lines derived via single seed descent. All the 284 RILs along with their parents were planted in two environments viz. Jiangpu Experimental Station, Nanjing, Jiangsu Province (Latitude 33°03′ N; Longitude 118°63′ E) in 2011 (JP2011) and Baima Experimental Station, Nanjing, Jiangsu Province (Latitude 31°62′ N; Longitude 119°18′ E) in 2020 (BM2020). Each line was planted in one-row plot (length × width, 1.5 × 1 m). Field management followed standard conditions in each location.

The third leaf from the top of each stem was taken from three plants in the field at V6 stage, then put in icebox and transported to the laboratory. Samples were dissected between the main vein and lateral vein near the base of the middle-leaflet of trifoliolate with 8 mm diameter puncher (avoiding primary veins) (Xing et al., 2013). Then the leaf discs were used to take photographs with an area of 12 mm⁻² under a Leica stereo microscope. The software ImageJ¹ was used to generate PL and PD. PD was converted from 12 mm⁻² to an area of 10 mm⁻² as the final density. As for PL, the pubescence on leaf surface was divided into two types: long and short, the length of three representative hairs of each type was measured, and average length was calculated by the weighted average method as

Where x¯is a weighted average of PL, x¯l is the average length of three representative long pubescence, f_l is the number of long pubescence, x¯sis the average length of three representative short pubescence, and f_s is the number of short pubescence.

R software was used to draw the frequency distribution of phenotypic data. The descriptive statistics, such as mean, maximum and minimum, coefficient of variation (CV), correlation analysis, and ANOVA of traits were calculated using SAS software (SAS Institute, 2010. SAS/STAT software version 9.2. SAS Institute Inc., Cary, NC, United States). The broad-sense heritability (h²) for individual environments (Eq. 1) and combined environments (CE; Eq. 2) were computed following the formula proposed by Nyquist and Baker (1991).

Where σg2 is the genotypic variance,σp2 is the phenotypic variance (PV), σg⁢e2is the genotype by environment interaction variance, σe2is the error variance, n is the number of environments, and r is the number of replications.

In the present study, a high-density genetic linkage map was constructed using restriction site-associated DNA sequencing (RAD-seq) (Wang et al., 2016). Briefly, restriction enzymes were used to digest the purified genomic DNA firstly, then ligated digested products with P1 adapter by T4 DNA ligase. Every 24 RILs were collected together and randomly sheared ultrasonically and used a purification kit to purify DNA fragments. Next, the fragment end was repaired with a Quick Blunting kit (NEB). Finally, the collected fragments were enriched by PCR amplification and purified by a QIAquick PCR purification kit. In addition, standardized samples were sequenced on HiSeq 2000 instruments. The soybean genome sequence (G. max, Wm82.a1. v1) was used as a reference to predict digestion sites. A total of 5,728 bin markers were obtained from 89,680 single nucleotide polymorphisms, spanning a total genetic distance of 2,204.6 cM with an average distance of 0.4 cM between neighboring bins. The linkage map of bin markers was constructed for the RIL population using R with the package LinkageMapView.

Two QTL mapping models were adopted to map additive effect QTLs in the present study to discover the genetic basis of pubescence development. Firstly, CIM was implemented in WinQTLCart 2.5 software with a 10 cM window at a walking speed of 1 cM to map additive effect QTLs. The log of odd (LOD) threshold was determined by 1,000 permutation tests for each trait with an experimental-wise error rate of P = 0.05 to determine whether the QTL was significant. The QTLs detected with overlapping or closely linked confidence intervals (CIs) in different environments were recognized as the same QTL.

Secondly, QTL Network v2.0 software with MCIM model was used to map additive effect QTLs with the critical F-value calculated with 1,000 permutation tests. In addition, the QTL effects were estimated using the Markov Chain Monte Carlo (MCMC) method with 20,000 Gibbs sampler iterations. The significance level configuration of candidate interval selection, putative QTL detection, and QTL effects were calculated with an experiment-wise type I error under α = 0.05. The above analyses were done for individual environments (JP2011 and BM2020), averages from JP2011 and BM2020 were designated as the CE.

The physical position of two flanking markers of major QTLs can be obtained by mapping sequencing data to Wm82.a1. v1. Both the model genes and annotation information within the physical genomic interval of major QTLs were obtained from SoyBase.² The expression data of model genes were downloaded from SoyBase³ and Phytozome.⁴ The genes that expressed in young leaf were further classified according to Protein Analysis THrough Evolutionary Relationships (PANTHER).⁵ Based on functional annotations, PANTHER analysis, and available literatures, some genes were selected as candidate genes and their relative expression levels available on SoyBase and Phytozome were heatmapped using TBtools (Chen et al., 2020).

To perform qRT-PCR, total RNA of leaf samples of two parents (NN 86-4 and PI 342618B) were isolated using RNA-prep Pure Plant Kit (TIANGEN DP-432, China) and full-length cDNA was reverse transcribed using a cDNA synthesis kit (Vazyme, R223) according to the protocol of the manufacturer. qRT-PCR was performed using ChamQ SYBR qPCR Master mix (Vazyme Q311) on Roche LightCycler 480 II. The housekeeping gene GmActin11 was used as the reference. Three biological replicates were conducted for each analysis. The relative expressions of selected genes were computed using a 2^–△CT method (Livak and Schmittgen, 2001). The primer sequences for qRT-PCR are listed in Supplementary Table 1. The cluster analysis was performed using the Neighbor-Joining method in MEGA6 (Tamura et al., 2013). The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1,000 replicates) are shown next to the branches.

Common cutworm pupa stock was obtained from Soybean Research Institute of Nanjing Agriculture University. Third-instar larvae with uniform size were used for the experiment as described by Xing et al. (2017). The third leaf from the top of the stem of four kinds of pubescence morphology [long and sparse (LS); short and sparse (SS); long and dense (LD); and short and dense (SD)] was used to feed CCW larvae (five larvae per replication), and 10 different representative lines were selected for each pubescence morphology (PL/PD class). This experiment was conducted with three biological replicates per line. The initial larval weight was recorded and measured 3 days after forcibly feeding. The increased larval weight was subjected to ANOVA as a resistance indicator and multiple comparisons of different pubescence morphology (LS, SS, LD, and SD) were conducted via the least significant difference at 5% probability in SAS software.

PL and PD of 284 RILs and their parents in JP2011, BM2020, and CE are presented in Figure 1. PL and PD of the male parent PI 342618B were ranged 0.18–0.33 mm and 54.53–90.83 hairs 10 mm⁻², respectively, which were shorter and denser than that of the female parent NN 86-4 (0.45–0.58 mm and 17.22–21.67 hairs 10 mm⁻²) (Table 1 and Figure 1). The variation between two parents offered a broader genetic resource among the RILs for quantitative trait analysis. The mean value of some RILs exceeded two parents in both directions, indicating that RILs showed transgressive segregation in PL and PD (Figure 1). The phenotypic variation of PL and PD among RILs showed continuous distribution, suggesting both two traits are controlled by multiple genes, and thus suitable for QTL mapping.

In addition to the above, ANOVA and correlation analysis suggested that both traits (PL and PD) were influenced by environment and genotype by environment interaction (G × E) (Table 1, Supplementary Table 2, and Figure 1). The high h² of the two traits (PL and PD) indicated that either of the traits is largely regulated by genetic factors. The correlation coefficients (r) of the same trait in different environments were ranged from 0.40 to 0.49, but a negative correlation was observed between PL and PD (r = −0.17 to −0.36; Supplementary Figure 1).

The bin-marker distribution on each chromosome is shown in Supplementary Table 3 and Supplementary Figure 2. A total of 16 QTLs comprising nine and seven for PL and PD, respectively, were identified by the CIM model with LOD (3.4–14.0) and phenotypic variation explained (R²) (4.2%–15.8%) (Table 2). The highest number of seven QTLs (qPL-12-1, qPL-12-2, qPL-12-3, qPL-12-4, qPD-12-1, qPD-12-2, and qPD-12-3) was mapped on Chr12 followed by three QTLs (qPL-1-1, qPL-1-2, and qPD-1-1) on Chr01, two QTLs (qPD-11-1, qPD-11-2) on Chr11 (Table 2 and Figure 2). These results suggested that PL and PD are largely controlled by Chr12, Chr01, and Chr11.

For PL, nine QTLs were identified on five chromosomes (Chr01, Chr03, Chr04, Chr12, and Chr14) with phenotypic variation ranging from 4.2 to 9.9% and LOD ≥ 3.4 (Table 2 and Figure 2). Among these, QTLs qPL-1-2, qPL-12-1, and qPL-12-3 were detected in JP2011 and CE with LOD values of 5.3–6.4, 6.8–7.9, and 4.6–6.8, R² values of 6.8–7.7%, 8.8–9.8%, and 5.8–8.2%, respectively (Table 2). qPL-1-1 was detected in JP2011 and accounted for an R² value of 9.9%. The remaining five QTLs have relatively smaller R² (4.2–8.8%). Aside from the additive effect of qPL-12-3 and qPL-12-4 that are negative, alleles for increased PL emanated from wild soybean PI 342618B, while others obtained increased PL additive effect from NN 86-4 (Table 2).

The QTLs related to PD were identified on four chromosomes (Chr01, Chr08, Chr11, and Chr12) with LOD scores ranging from 4.1 to 14.0, which could explain phenotypic variation from 4.5 to 15.8% (Table 2 and Figure 2). Among them, qPD-1-1 could be detected in both environments and CE with LOD_vs_R² values of 7.4–11.5_vs_9.1%–12.1%, qPD-11-1, qPD-11-2, and qPD-12-2 could be detected in one environment and CE with LOD_vs_R² values of 4.2–6.7_vs_5.0–6.7%, 4.1–5.6_vs_4.8–5.7%, and 14.0_vs_15.4–15.8%, respectively. Besides, qPD-12-3 was detected in only one environment (JP2011) but caused the highest R² of 13.2%. The remaining two QTLs (qPD-8-1 and qPD-12-1) were detected in single environment and could explain relatively lower phenotypic variation (Table 2). Interestingly, except qPD-8-1 with a positive additive allele that increased PD allele from NN 86-4, all additive effects of other QTLs were negative with alleles for increasing PD from PI 342618B, suggesting PD is a domestication related trait which may have been lost or reduced during domestication. Hence the wild soybean PI 342618B contains beneficial alleles which could be exploited to increase PD in its domesticated progenies (Table 2).

By comparing QTLs for PL and PD, qPL-12-1 and qPD-12-2, qPL-12-2 and qPD-12-3 overlapped, respectively. These four QTLs were flanked by markers bin3269-bin3276, bin3272-bin3277, bin3282-bin3285, and bin3284-bin3288, respectively (Table 2 and Figure 2). The increased pubescence traits alleles of PL and PD were from different parents, hence they could be the same QTL with pleiotropic effect, and explain the negative correlation between PL and PD.

In all, one additive QTL (qPL-16-1) of PL and four additive QTLs (qPD-1-2, qPD-2-1, qPD-11-1, and qPD-12-2) of PD were detected by MCIM implemented in QTL Network v2.0 software.

For PL, qPL-16-1 could cause phenotypic variation of 3.0%. It is a novel locus detected for the first time with an increased PL additive effect from NN 86-4 (Table 3). A total of four QTLs were detected related to PD accounted for 0.8%–12.8% phenotypic variation, among which qPD-12-2 could explain the phenotypic variation of 12.8%. All the QTLs inherited their increased PD alleles from PI 342618B (Table 3), supporting our earlier assertion that pubescence may be one of the domestication syndrome traits in soybean. qPL-12-1 and qPD-12-2 by CIM were overlapped with qPD-12-2 by MCIM (Tables 2, 3), hence this locus was considered as the major QTL for pubescence development in the present panel. In addition, qPD-11-1 was detected by CIM and MCIM, respectively (Tables 2, 3). Therefore this region was considered as major QTL for regulating PD in this population.

For PL: qPL-1-2, qPL-12-1, and qPL-12-3 were detected in JP2011 and CE, which were considered as major QTLs. These three QTLs contain 131, 60, and 80 genes (total 271), respectively (Wm82.a1. v1). For PD, qPD-1-1 was mapped in both JP2011 and BM2020 and CE, while qPD-11-1 and qPD-12-2 were detected by both CIM and MCIM. Besides, qPD-12-3 with the LOD value of greater than 10 caused 13.2% phenotypic variation in PD and overlapped with Pubescence density 2-7 (Du et al., 2009), hence this could be considered as a major QTL. qPD-11-2 was not considered as a major QTL due to lower phenotypic variation. Within the genomic regions of qPD-1-1, qPD-11-1, qPD-12-2, and qPD-12-3, 70, 75, 41, and 39 model genes (total 225) were downloaded, respectively.

A total 247 of 271 genes for PL and 200 of 225 genes for PD were informatively annotated, respectively (Supplementary Tables 4, 5). Expression data of distinct tissues in soybean have been completed in previous research (Severin et al., 2010). A total 203 and 151 genes of PL and PD with expression in young leaf were selected for PANTHER analysis (Wm82.a2. v1) (Supplementary Tables 6, 7). For PL, 93 out of the 203 genes were included in PANTHER protein classes and involved in 14 pathways. Then 17 out of the 93 genes were considered as the candidate genes according to literatures and annotation information (Table 4). Furthermore, there were other four candidate genes that did not include in protein classes (Table 4). For PD, 70 out of the 151 genes were included in PANTHER protein classes and involved in nine pathways. A total of nine genes were considered as the candidate genes according to literatures and annotation information. Except for the above genes, three candidate genes were not included in protein classes (Table 4). Among above 21 and 12 candidate genes of PL and PD, respectively, Glyma.12g185500, Glyma.12g188600, and Glyma.12g188800 were responsible for both PL and PD, suggesting a possible pleiotropic effect of some candidate genes.

Most of these 30 (three genes are the same for both PL and PD) candidate genes with relative higher expression in young leaf or shoot apical meristem (SAM) (Supplementary Table 8 and Figure 3). Their expression was measured subsequently by qRT-PCR in the leaves of two parents: PI 342618B and NN 86-4 (Figures 4A,B). Eight and five genes of PL and PD, respectively, were expressed differentially by more than 2-fold between two parents. For PL, Glyma.01g198100 expressed more than 50 folds higher in NN 86-4 than PI 342618B (Figure 4A). For PD, Glyma.12g195900 has an expression level of more than 9-fold in PI 342618B compared with NN 86-4 (Figure 4B). Also, Glyma.12g195900 is the homologous gene of CYCU1 which could promote meristem cell division in Arabidopsis (Peng et al., 2014). Therefore, Glyma.01g198100 and Glyma.12g195900 were considered as the important candidate genes for PL and PD, respectively.

Glyma.01g240100, one of the 80 model genes of qPD-1-1, has been identified as Pd1 due to a T to C single nucleotide polymorphism (SNP) in the last exon (Liu et al., 2020). Only one G to A synonymous mutation was detected in the seventh exon between PI 342618B and NN 86-4 in the present study (Figure 4C and Supplementary Figure 3). Glyma.12G187200, a candidate gene of major QTL qPL-12-1 and qPD-12-2 with incompletely annotated information, was identified as Ps due to different copy numbers and expression level (Liu et al., 2020). The expression level of Glyma.12G187200 in the leaves of NN 86-4 was fivefold greater than in PI 342618B (Figure 4D). These results suggested that Glyma.01g240100 may not be the candidate gene of qPD-1-1, and Glyma.12G187200 was the candidate gene that not only controls PD of soybean leaves, but it also contributes to PL in the present study. Then a cluster analysis of the above candidate genes was conducted with their homologs, which have been known to be associated with pubescence development in several species (Supplementary Figure 4). Most candidate genes clustered together with their homologs, with seven were reliable (bootstrap values > 65%), such as Ps and Pd1 (Supplementary Figure 4). Among them, Glyma.01g195900 and Glyma.12G187200 (Ps) of PL and Glyma.12g194400, Glyma.12G187200 (Ps), Glyma.12g195200, and Glyma.12g195900 of PD were suggested as important candidate genes due to significant differences in the expression between parents. Based on qRT-PCR and cluster analysis, three (Glyma.01g195900, Glyma.01g198100, and Glyma.12G187200) and four (Glyma.12G187200, Glyma.12g194400, Glyma.12g195200, and Glyma.12g195900) important candidate genes were identified for PL and PD, respectively.

To determine the relationship between pubescence morphology and resistance to CCW, an antibiotic test was carried out using lines with four types of pubescence morphology and the typical photographs of pubescence are shown in Figures 5A–D. There was a significant difference between increased larval weight of feeding with leaves with SD pubescence (similar to the pubescence morphology of wild soybean PI 342618B) and leaves with the other three pubescence morphology types. The former (SD, viz. wild soybean pubescence morphology) had the strongest resistance to CCW, followed by LD, SS, and LS (Figures 5E–I). The increased larval weight had positive and negative correlations with PL and PD, respectively (Table 5).

Pubescence can form a physical barrier that protects plants against various biotic and abiotic stresses (Kennedy, 2003; Kang et al., 2016). The effectiveness of this preventive mechanism depends on the length and density of pubescence (Hua et al., 2021). Therefore, it is necessary to understand the genetic mechanisms associated with polygenic quantitative characters PL and PD. Although several QTLs related to PL and PD have been identified and reported over the past decades, most of them have large interval regions due to the small mapping population (<200 lines) and low-density genetic map based on simple sequence repeat (SSR) markers (Du et al., 2009; Oki et al., 2012). These may be difficult to be used in practical plant breeding and predicting probable candidate genes.

In the present study, two QTL mapping methods were used to complement and validate the results of each other and improve the accuracy of QTL mapping. A total of nine and seven QTLs were detected of PL and PD by CIM, respectively, while one and four QTLs were mapped by MCIM. By comparing the results of two methods, two QTLs (qPD-11-1 and qPD-12-2) related to PD were detected in both methods, and qPD-12-2 (Chr12) was the same as previously reported QTLs viz. Pubescence density 2-8, Pubescence density 3-2, and PD12-1 (Du et al., 2009; Oki et al., 2012; Xing et al., 2013), suggesting these loci play a role in pubescence development. Besides, the CIM method detected more QTLs with higher additive effect and R² value while the MCIM method mapped QTLs with narrow CIs. It was possible to miss some important loci if only the MCIM method was used, thus, it is better to use two methods.

By comparing QTLs detected in two methods with the reported ones in previous studies, qPL-12-2 (Chr12) with an interval of 35,364,334–35,611,487 bp overlapped with Pubescence length 1-2, which had a larger CI (35,108,089–36,780,375 bp) (Oki et al., 2012). qPL-12-3 (Chr12) was overlapped with PL12-1 (Xing et al., 2013). The remaining eight QTLs related to PL were detected for the first time in the present study (Tables 2, 3 and Figure 2). For PD, a total of nine QTLs were detected by CIM and MCIM. In addition to qPD-12-2 mentioned in the previous paragraph, qPD-12-3 (Chr12) was overlapped with Pubescence density 2-7 (Du et al., 2009); qPD-1-1 (Chr01) and qPD-2-1 (Chr02) were same to Pubescence density 3-1, PD1-1 (Oki et al., 2012; Xing et al., 2013) and Pubescence density 2-5 (Du et al., 2009), respectively, but in a narrower genomic region in the present study. The remaining five QTLs were detected for the first time. Most QTLs were novel indicating the distinct and abundant genetic architecture of pubescence in wild soybean. It also suggests the need to utilize more diverse parents to develop a mapping population to reveal the complex genetic basis of pubescence development in soybean and provide more valuable information for the gene identification related to pubescence development. In addition, the majority of QTLs identified in the present study were in small physical genomic regions, suggesting the importance of using a high-resolution genetic map for QTL detection and candidate gene exploration.

Although many studies had demonstrated that dense and long pubescence have higher resistance to abiotic stress (Turnipseed, 1977; Gunasinghe et al., 1988), the purpose of soybean breeding is not always to increase density and length of pubescence. It is important to keep PL and PD within a suitable range for the better growth and development of plants. In addition, it was found that there was a negative correlation between PL and PD in the present study, hence, materials with dense and long pubescence may not be easy to obtain. Our results demonstrated that soybean leaves with SD viz. wild soybean pubescence morphology instead of LS viz. cultivar soybean pubescence morphology had the stronger resistance to CCW thus the former can be used in soybean breeding.

It is of great significance for both theoretical research and breeding practice to identify the candidate genes of major QTL regions of pubescence traits in soybean. Many factors were identified to be related to trichome development in other species, providing the useful information to explore candidate genes of soybean pubescence development. The mechanism of Arabidopsis trichome development has been comprehensively explained (Shang et al., 2020). The core regulatory components are the R2R3-MYB/basic helix-loop-helix (bHLH)/WD complex (Ishida et al., 2008). R3-MYB negatively regulates trichome formation by competing with R2R3-MYB for binding to bHLH (Hülskamp et al., 1994; Schellmann et al., 2002; Esch et al., 2003). In tomato, bHLH TF (SlMYC1) and R2R3-MYB TFs (SlTHM1 and SlMYB52) play an important role in the formation of trichomes (Xu et al., 2018; Yuan et al., 2021). The fiber initiation and elongation are somehow similar to trichome development. In cotton, R2R3-MYB TF GhMYB109 is required for cotton fiber development (Pu et al., 2008) and the homologues of GhMYB109 in tomatoes might also participate in the regulation of trichome elongation (Hua et al., 2021). Therefore, MYB TF, bHLH TF, and WD play an extremely important role in both unicellular and multicellular trichome development.

Recently, actin filaments and microtubules were reported to play coordinated but distinct roles in the formation of tomato trichome (Chang et al., 2019). Both in Arabidopsis and tomatoes, mutations in genes of SCAR/WAVE complex could lead to distorted trichomes (Kang et al., 2016; Chang et al., 2019; Li et al., 2019). In soybean, GmNAP1 was involved in actin filament assembling during trichome and pavement cell development (Campbell et al., 2016; Tang et al., 2020). Thus, actin and microtubules were identified as having an undeniable role in trichome development in recent years.

Additionally, a set of C2H2 zinc finger TFs, such as GIS, GIS2, GIS3, ZFP5, and Hair (H) gene, were detected to be involved in trichome development in Arabidopsis and tomatoes, respectively (Gan et al., 2007; Zhou et al., 2011; Sun et al., 2015; Chang et al., 2018). AP2 and AP2/ERF TFs: TAR1 and Hairy Leaf 6 (HL6), play an important role in trichome development in A. annua (Tan et al., 2015) and rice (Sun et al., 2017; Shang et al., 2020). Besides, the WUS-type homeobox gene OsWOX3B was found to be required for macro-hair initiation and trichome development in rice (Zhang et al., 2012). Cyclins were involved in the transition of the cell cycle and function as positive regulators of cell proliferation in eukaryotes (Meijer and Murray, 2001) and a B-type cyclin gene, SlCycB2, plays key roles in trichome initiation in tomatoes (Gao et al., 2017). These results suggested that the above factors may be functional during trichome development.

A total of 22 and 13 candidate genes (together with Ps) were identified for PL and PD, respectively, based on PANTHER analysis, expression data, and literatures in the present study (Table 4). It should be noted that genes within the physical genomic interval that are not annotated and/or have no expression in young leaves may be ignored. A cluster analysis was conducted of these candidate genes and the homologs mentioned above (Supplementary Figure 4). The genes Ps and Pd1, which have known to be related to soybean PD, were clustered with homologous genes of other species, indicating the reliability of this analysis method. However, more study is needed for their functional validation.