K099 (shorter primary root) and Fendou 16 (longer primary root) are parental genotypes of the mapping population used to identify the qRL16.1 a primary root length locus (Chen et al., 2021). Residual heterozygous line-derived near-isogenic lines (NILs): Near-isogenic lines differing at qRL16.1, which were selected from the same recombinant inbred line (RIL) but heterozygous at qRL16.1 (Chen et al., 2021). NIL-F had the Fendou 16 homozygous genotype, and NIL-K had the K099 homozygous genotype at qRL16.1. Advanced backcross lines BC4-F and BC4-K: BC4-F had the Fendou 16 homozygous genotype and BC4-K had the K099 homozygous genotype at the mapped qRL16.1 region in the background of K099 (Chen et al., 2021). A total of 76 germplasm accessions from a mini-core collection of the world soybean germplasm (Kaga et al., 2012), together with K099 and Fendou16, were used for phenotyping and analysis of the association of root traits with identified polymorphisms.

Root samples of NIL-K and NIL-F were collected from 7-day-old seedlings. Each sample comprised roots collected from five individual plants, and three biological replicates were analyzed. Total RNA was extracted from each sample and cDNA was synthesized with adapters and sequenced using the Illumina HiSeq 2500 analyzer at BGI Technologies (Shenzhen, China). The adapter reads and low-quality reads (bases with quality values ≤ 5) were filtered out of the raw data to obtain high-quality (clean) reads, which were mapped to the soybean reference genome (G. max Wm82.a2.v2) downloaded from Phytozome (Goodstein et al., 2012) using HISAT v2.0.4 (http://www.ccb.jhu.edu/software/hisat). Significant differentially expressed genes (DEGs) were obtained with fold changes ≥ 2.0 or ≤ −2.0 for up- or down-regulated genes and adjusted P-value ≤ 0.05.

Genomic sequences containing coding sequences and upstream promoter regions of two DEGs, Glyma.16g108500 and Glyma.16g112400, were sequenced from both parents, i.e., K099 and Fendou16, using Sanger sequencing. For real-time PCR based expression analysis of DEGs, primary root tip and lateral root tissues from K099, Fendou16, NIL-K, NIL-F, BC4-F, and BC4-K were collected from 3-week-old plants and frozen in liquid nitrogen. RNA isolation from three biological replicates of each sample was performed using the RNeasy plant mini kit (Qiagen, Germany). The isolated RNA was treated with DNAse I using the RapidOut DNA Removal Kit (Thermo Scientific, USA). Purified RNA was used for cDNA synthesis using the PrimeScript RT master mix (TAKARA BIO, Japan). Comparative quantitative real-time PCR (RT-PCR) was performed using Powerup SYBR Green Master mix (Applied Biosystems, USA) for selected candidate genes. GmActin (Glyma.18g290800) was used as an internal reference gene. The sequences of primers used for RT-PCR are given in Table 1 . The 2^-ΔΔCT method was used for analyses of relative gene expression levels (Livak & Schmittgen (2001).

A total of 76 accessions from the world soybean mini-core collection of Japan, along with K099 and Fendou16, were used for phenotyping of primary root length, total root length, root volume, surface area, and root tips using hydroponic culture (Chen et al., 2021). The roots of 15-day-old plants (n = 4) were scanned and analyzed using WinRHIZO software (Regent Instruments Inc., Canada). Phenotyping data from these accessions were used for the marker-trait association analysis. The top four accessions for longer primary roots and the top four accessions for shortest primary roots were used for the validation of four differentially expressed genes.

Accessions of the world soybean mini-core collection of Japan, together with K099 and Fendou16, were genotyped with one SNP (RL1SNP968_T/C) and one Indel (RL1Indel363_–/AA) marker selected from the promoter region of Glyma.16g108500. Allele-specific PCR primers were used for the amplification of specific SNP or Indel alleles and were separated on 8.0% polyacrylamide gel. The DNA band pattern was visualized using a Pharos FX™ Molecular Imager (Bio-Rad, Tokyo, Japan) instrument for scoring alleles. The Welch’s t-test was used to test the (null) hypothesis that two populations formed by alternate alleles have equal means.

In order to identify CEP-like genes in soybean, protein sequences from 16 known CEP genes of Arabidopsis thaliana were used as queries in a BLASTP search against the Williams 82 genome Wm82.v2 in SoyBase (www.soybase.org). The 14 homologs identified in the soybean genome were used again for iterative BLASTP to identify more CEP-like genes in the soybean genome. Thus, 23 more matching gene sequences were identified. The 37 matching sequences were manually filtered for proteins with an N-terminal signal peptide sequence for the secretory pathway and a conserved C-terminal ‘CEP-like’ sequence motif ‘FRPToPGHSPGVGH’ (at 70% consensus cutoff), using multiple sequence alignment by CLUSTALW and motif prediction by MEME Suite tool (Bailey et al., 2015). Finally, 25 CEP-like genes were selected and used for phylogeny and motif analysis. MEGA11 software was used to construct the phylogeny of 25 CEP-like gene protein sequences (Tamura et al., 2021). The neighbor-joining clustering method with 10,000 bootstrap iterations was used to build a consensus phylogenetic tree. The MEME Suite tool was used for consensus motif detection with the following settings: motif site distribution as any number of sites per sequence, motif width 6 to 50, minimum site for motif 2, and maximum number of motifs of 10 (Bailey et al., 2015).

The conserved regions of CEP for Glyma.16g108500 and its homolog Glyma.16g108700 were predicted, and putative post-translational modifications of proline amino acid residues, that is, hydroxyproline, were included for synthesis ( Table 2 ). Synthetic peptides Glyma.16g108500 (CEP-RL1) and Glyma.16g108700 (CEP-RL2) were synthesized artificially at Eurofins, Japan. The seeds of the short-rooting genotype K099 were sown in vermiculite and young seedlings 4 days after sowing were used for the synthetic peptide assay. The 4-day seedlings were transferred to 0.5× Hoagland solution (pH 6.5). For each treatment, eight seedlings were placed in a plastic box (13 cm × 8.5 cm × 6 cm) containing 500 ml Hoagland solution. The outer surface of the box was covered with aluminum foil to prevent root exposure to the light. Following concentrations were used for the treatment in Hoagland solution: i) Control: No CEP; ii) CEP-RL1: 1 µM; iii) CEP-RL1: 2 µM; iv) CEP-RL2: 1 µM; v) CEP-RL2: 2 µM. The seedlings were incubated in an incubation chamber kept at 22° C night and 24°C day temperatures and supplied with 14 h of fluorescent light. Roots were harvested after 10 days of peptide treatment and scanned for recording of root traits using WinRHIZO (Regent Instruments Inc., Canada). The recorded root traits were primary root length, total root length, and number of root tips.

Two time stages, 4-day, and 7-day, were selected for gene expression analysis of Glyma.16g108500 and Glyma.16g108700 in parental genotypes K099 and Fendou16. Tissues were collected at 4 days after sowing, from the shoot apex (2 mm) and primary root tip (2 cm), and at 7 days after sowing, from shoot apex (2 mm), the primary root tip (2 cm), and lateral root tips (2 cm), for tissue-specific gene expression analysis.

RNA-Seq analysis of the root transcriptome of two NILs identified total 317 genes upregulated and 192 downregulated in NIL-F compared to NIL-K ( Supplementary Figure 1A ). Differential expression analysis for genes from the genomic region of primary root length QTL qRL16.1 (Chr.16:27158682.28953489), among two NILs, identified only two genes, Glyma.16g108500 and Glyma.16g112400, showing a significant fold difference in their expression levels ( Supplementary Table 1 ). The gene expression analysis by RT-PCR identified that the expression of Glyma.16g108500 was 3.0 fold higher in Fendou16 tissues compared to K099 tissues, and it was 4.47 fold higher in NIL-F tissues than in NIL-K tissues. Conversely, the expression level of Glyma.16g112400 was 18.11 fold lower in Fendou16 tissues compared to K099 tissues, and it was 7.21 fold lower in NIL-F tissues than in NIL-K tissues ( Supplementary Figures 1B, C ). The genomic sequence of two DEGs, Glyma.16g108500 and Glyma.16g112400, including a 2-kb promoter region, was sequenced from both parents, Fendou16 and K099. No sequence differences were identified in the coding DNA sequence of both genes between the two parents. Upstream 2-kb sequencing of Glyma.16g108500 identified two base insertions at the 363 position (–/TA), a single nucleotide polymorphism (SNP) at -968 (T/C), and a single base insertion at the 1038 position (A) in Fendou 16. The 2 kb upstream sequencing of Glyma.16g112400 did not detect any sequence difference between the two parental genotypes.

According to gene annotations assigned in SoyBase (www.soybase.org), Glyma.16g108500 was annotated as an uncharacterized protein, root development, and DNA-directed RNA polymerase II subunit RPB1-like protein, while Glyma.16g112400 was annotated as an uncharacterized protein. Furthermore, the analysis of protein homology in the InterPro protein database identified that Glyma.16g108500 encodes a C-terminally encoded peptide (CEP), a C-terminally conserved secreted protein in Arabidopsis (https://www.ebi.ac.uk/interpro/entry/InterPro/IPR033250/). Based on protein homology, two more paralogs of Glyma.16g108500 were also identified in close proximity to Glyma.16g108500 in the genomic region of qRL16.1. These two paralogs, Glyma.16g108700 and Glyma.16g108400, showed 72% and 80% homology of the protein sequence with Glyma.16g108500, respectively.

The two genes Glyma.16g108500 and Glyma.16g112400 were validated for differential expression by real-time PCR in primary root tip of two parents, a pair of contrasting NILs, and a pair of contrasting BC₄ backcross lines. Expression analysis identified that Glyma.16g108500 showed higher expression in longer rooting genotypes, whereas Glyma.16g112400 showed higher expression in shorter rooting genotypes ( Figure 1 ). In the case of Glyma.16g108500, although the expression fold difference in the contrasting alleles NILs and BC₄s was similar to that of the parents, the level of expression in both pairs of contrasting lines was slightly reduced compared with that in the two parents. The analysis of gene expression for two Glyma.16g108500 paralogs, Glyma.16g108700 and Glyma.16g108400, among three pairs of contrasting genotypes identified that Glyma.16g108700 and Glyma.16g108400 also showed higher expression in longer rooting genotypes; however, the level of expression fold difference was relatively higher in the case of Glyma.16g108700 ( Figure 1 ).

The descriptive statistics of different root traits in 76 soybean accessions, Fendou 16 and K099, at 15 days of growth are given in Table 3 . The primary root length ranged from 16.3 cm in CHUUHOKU 2 to 51.45 cm in BARITOU 3A, with a mean root length of 34.4 cm. The total root length varied from 142.8 cm to 606.4 cm with a mean value of 402.7 cm.

The five high and five low contrasting accessions for primary root length and total root length traits along with Fendou 16 and K099 were again evaluated. In this confirmation experiment, similar results were obtained for short- and long-rooting genotypes. Among the top five accessions for primary root length, BARITOU 3A showed the highest primary root length of 53.1 cm, followed by RINGGIT (51 cm) and Fendou 16 (49.2 cm). From this experiment, the top four accessions for longer primary root and the four best accessions for short primary root were selected for validation of four differentially expressed genes ( Table 4 , Figure 2 ).

The four differentially expressed genes identified among three contrasting pairs of high and low root length lines from the genomic region of qRL16.1 were further analyzed for their expression in the top four long rooting and top-four shortest rooting germplasm selected from the world core collection of soybean. Gene expression analysis using real-time PCR showed that Glyma.16g108500 and Glyma.16g108700 had relatively higher expression in longer rooting genotypes than in shorter rooting accessions, except for accession Heukdaelip ( Figure 3 ). However, Glyma.16g108400 and Glyma.16g112400 did not show differential expression between the two contrasting groups in accordance with root length and expression pattern observed in the NILs, thus demeaning their role in the regulation of the development of primary root in soybean ( Figure 3 ).

Coding region of both DEGs, Glyma.16g108500 and Glyma.16g108700, was sequenced but no polymorphism was detected between K099 and Fendou 16. The two polymorphisms, SNP at -968 (RL1SNP968_T/C) and Indel at -363 (RL1INDEL363_–/TA), identified in the upstream sequence of Glyma.16g108500, were analyzed in the germplasm to test the association of allelic variation present in the proximity of two DEGs identified in contrasting germplasm. Genotyping of the RL1SNP968_T/C and RL1INDEL363_–/TA markers showed identical marker genotypes among 78 genotypes with 36% and 64% allele frequencies for the K099 and Fendou16 alleles, respectively. The Welch t-test showed significant differences in primary root length and total root length between two populations formed by alternate alleles of RL1INDEL363_–/TA ( Table 5 , Figure 4 ).

A total of 25 CEP-like genes were identified that contain an N-terminal signal peptide and the CEP-like conserved motif ‘FRPToPGHSPGVGH’ (at consensus cutoff of 70%). These 25 CEP-like genes were distributed on seven chromosomes. Except for Chr.06 and 09, other CEP-like genes were located in tandem, suggesting that their copy number increased by tandem duplication. The highest tandem copies (eight) were present on Chr. 17. Phylogeny analysis using the neighbor-joining clustering method showed that most CEP-like proteins were clustered together in group I, while two CEP-like proteins from Chr.01 and Chr.09 were grouped in cluster II. The two CEP-like proteins of Glyma.16g108400 and Glyma.16g108500 from Chr.16 were clustered in group III, whereas the CEP protein from Glyma.16g108700 was most divergent ( Figure 5 ). Consensus motif analysis using MEME Suite identified consensus motif distribution patterns among 25 CEP proteins with an N-terminal signal peptide motif having only one site for each sequence, while a C-terminally encoded peptide conserved motif (Motif 1) was present as a single copy in 19 members and as two copies in six members ( Figure 6 ). Two conserved domains (VLWYTPND and MINIHPPPAIPRSPQ) are specifically found in three CEP-like genes of Chr.16 and two members of group II of the phylogenetic tree, that is, Glyma.01G152800 and Glyma.09G218000 ( Figure 6 ). Three proline residues were present in the 15 amino acid CEP domain of Glyma.16g108500, while only two proline residues were present in the CEP domain of Glyma.16g108700.

Because the analysis of gene expression in contrasting germplasm excluded Glyma.16g112400 and Glyma.16g108400 from the possibility of their positive role in regulating the development of primary root length, we focused on Glyma.16g108500 and Glyma.16g108700 for further studies using 15 amino acid synthetic peptide of the predicted CEP domain for these genes. Two different concentrations of two CEPs for these genes (CEP-RL1 and CEP-RL2) were applied in Hoagland solution-based hydroponic culture (pH 6.5) for the short-rooting genotype K099. The results of this experiment showed that the length of the primary root was significantly higher for CEP-RL1 at concentrations of 1 µM and 2 µM s, whereas for CEP-RL2, the length of the primary root was significantly higher at a concentration of 2 µM concentration, compared to the control ( Figure 7 ). Analysis of the effect of CEP treatment on total root length and root tips showed that at higher concentrations of each CEP (2 µM), total root length and root tips were significantly lower in CEP-RL1 treatment, but not in the case of CEP-RL2 treatment ( Supplementary Figure 2 ).

In order to observe the expression pattern of two CEP genes in meristematic tissues of parental genotypes, the shoot apex and primary root tip were analyzed 4-day after sowing. Lateral roots appeared properly after 6–7 days of sowing, and so the shoot tip, the primary root tip, and the lateral root tips were analyzed at the 7-day stage. The analysis of gene expression for Glyma.16g108500 in the shoot and the primary root tip of 4-day seedlings showed expression only in the primary root tip, where the expression was 2.7-fold higher in Fendou 16 than in K099 ( Figure 8A ). For Glyma.16g108700, expression was observed in the primary root, but the expression level was similar in both genotypes ( Figure 8A ). No expression was observed in shoot apex for either gene. These results indicated that the genes Glyma.16g108500 and Glyma.16g108700 are root-specific in the early stage of root development. In the primary root and lateral root tissues of 7-day seedlings, both genes showed expression and the expression level was higher in the primary root than in the lateral roots ( Figure 8B ). Fendou 16 showed a 3.2-fold higher expression of Glyma.16g108500 and 1.8-fold higher expression of Glyma.16g108700 in primary root than K099, but no differential expression was observed in lateral roots. In the shoot apex of 7-day seedlings, no expression was observed for Glyma.16g108500, but low expression was observed for Glyma.16g108700. These results also showed the involvement of these two genes in the early stage of primary root development.