Detection of Candidate Genes and Development of KASP Markers for Pod Length and Pod Width by Combining Genome-Wide Association and Transcriptome Sequencing in Vegetable Soybean

Abstract

Vegetable soybeans are one of the most important vegetable types in East Asia. The yield of vegetable soybeans is considerably influenced by the size of their pods. To facilitate the understanding of the genetic basis of the pod length and width in vegetable soybeans, we conducted a genome-wide association study (GWAS) and transcriptome sequencing. Four quantitative trait loci, namely, qGPoL1, qGPoL2, qGPoW1, and qGPoW2, were mapped via GWAS analysis. Through the integration of gene function annotation, transcriptome sequencing, and expression pattern analysis, we identified Glyma.06G255000 and Glyma.13G007000 as the key determinants of the pod length and width in vegetable soybeans, respectively. Furthermore, two kompetitive allele-specific polymerase chain reaction (KASP) markers, namely, S06-42138365 (A/T) and S13_628331 (A/T), were developed and effectively validated in 27 vegetable soybean accessions. Overall, our research identified genes that regulate the pod length and width and determined KASP markers for molecular marker-assisted selection breeding. These findings have crucial implications for the improvement of soybean crops and can contribute to the development of efficient breeding strategies.

1. Introduction

Vegetable soybean (Glycine max (L.) Merr.) is a type of soybean cultivated for its suitability for vegetable consumption [1]. This crop is typically harvested during the R6 to R7 growth stages, when pods are plump and have reached 80–90% of their full size, and the pods and seeds display a vibrant emerald green color [2,3]. Vegetable soybeans are a delicate and nutritious food and abundant in proteins, unsaturated fatty acids, vitamins, plant fiber, plant hormones, and various minerals, such as calcium, iron, and phosphorus. These nutrients are easily absorbed and utilized by the human body and play a crucial role in the regulation of one’s dietary structure and enhancement of nutritional status [4]. Vegetable soybean varieties have high requirements for the appearance and quality of pods and seeds and typically require a larger size. The size of the pods, which is measured by their length and width, is a crucial factor in determining their storage capacity and serves as an important indicator for the evaluation of the overall appearance and yield quality of fresh soybeans [2].

Thus far, numerous quantitative trait loci (QTLs) associated with pod size have been identified [5]. Xie et al. utilized soybean chromosome segment substitution lines to map the QTLs responsible for soybean two-seed pod length (TSPL) and two-seed pod width (TSPW) [6]. They identified three QTLs for TSPL and five QTLs for TSPW. Subsequent investigation indicated that Glyma.11G051600 may regulate TSPL and TSPW [6]. Kulkarni et al. utilized the recombinant inbred line (RIL) population generated from “PI 483463” and “Hutcheson” to conduct a QTL analysis of the soybean 100-grain weight. Fourteen QTLs were identified, and they explained 3.83–12.23% of the total phenotypic variance [7]. The RIL population (“TU” × “Musica”) was utilized to identify 35 QTLs associated with pod morphological traits.

Genome-wide association analysis (GWAS) is a strategy based on linkage disequilibrium used for locating QTLs [8]. Given the advantages of a low false positive rate, high resolution, and high throughput [9], GWAS, which has been widely used in soybeans with different quantitative traits, can simultaneously locate different traits in various groups. Fang et al. identified 245 significant genetic loci and 84 agronomic traits with a panel consisting of 809 soybean accessions by GWAS [10]. Liang et al. performed a GWAS study and identified a major QTL named GmSW17 (Seed Width 17), which determines soybean seed width/weight [11]. Three QTLs associated with low-temperature tolerance were identified via GWAS through an association panel comprising 260 soybean accessions, and fifteen candidate genes were detected through haplotype and gene expression analyses [12]. However, thus far, studies reporting GWAS research related to vegetable soybean pods are limited. Li et al. identified 4 SNPs associated with the 100-pod fresh weight and 15 SNPs associated with the 100-seed fresh weight through GWAS analysis [13].

Gene Ontology (GO) enrichment analysis has revealed the considerable enrichment of phenylpropanoid metabolic process and auxin response genes, which indicates that these pathways contribute to the regulation of pod morphology [14]. Nevertheless, the molecular mechanisms governing the size of vegetable soybean pods remain to be fully elucidated.

Marker-assisted selection (MAS) uses molecular markers to guide breeding. This method involves the use of molecular markers to track genes or gene segments associated with a trait of interest to aid in the selection of individuals with superior traits [15]. The advancement of biotechnology has led to an increase in the research and development of various types of molecular markers. Kompetitive allele-specific polymerase chain reaction (KASP) marker is a novel approach to single-nucleotide polymorphism (SNP) typing. This method relies on allele-specific amplification and high-sensitivity fluorescence detection. The genotype of the amplified product can be verified through fluorescence detection. The key advantage of KASP is its low cost and high throughput, which allow for precise allelic genotyping at SNPs and insertion/deletion (InDel) loci through specific base pair matching at the ends of primers [16]. This method has been extensively applied in the fields of crop molecular marker-assisted breeding, germplasm resource identification, genetic map construction, and seed purity identification [17].

To elucidate the molecular mechanism of pod development in vegetable soybean, we conducted a GWAS and transcriptome analysis and identified candidate genes related to the soybean pod length and width. Subsequently, KASP molecular markers were developed. The results of this study will lay a good foundation for the molecular design breeding of soybeans.

2. Materials and Methods

2.1. Plant Materials and Phenotyping

A panel of 301 soybean accessions was used for genome-wide association analysis (Table S1). These accessions were collected from various provinces, including Jiangsu, Anhui, and Heilongjiang, etc. All accessions were cultivated at three locations (Xuanwu, Liuhe, and Lishui) in Jiangsu Academy of Agricultural Sciences Nanjing, China, during 2022 and 2023. The latitude of the three experimental locations are as follows: Xuanwu (118°50′ E, 32°03′ N, 26 m), Liuhe (118°34′ E, 32°37′ N), and Lishui (119°02′ E, 31°35′ N). The average altitudes of Xuanwu, Liuhe, and Lishui were 19, 20, and 39 m, respectively. All experiments followed a randomized block design with three replicates, and twelve plants were planted per replicate for each accession.

For collection and recording of the data on soybean pods, we investigated the length and width of the pods uniformly at the R6 stage. Then, 20 two-seed pods were randomly selected from each line and measured for their length and width. Statistical parameters, including the mean and standard deviation (SD), of soybean accessions were calculated using Microsoft Excel. In this study, 6 phenotype data are available for one soybean accession, and another accession corresponding to one phenotype data point was required for GWAS analysis. Therefore, we calculated the best linear unbiased estimate (BLUE) value as a phenotype value. To calculate the BLUE value and h² using the lme4 package (version 1.1-35.5), we calculated the Spearman coefficient for pod length and width traits and visualized them using the corrplot package (version 0.92) [18].

2.2. Genotypic Analysis and GWAS

The genotypic data on 301 soybean accessions were described in our previous study [19]. Briefly, DNA libraries were sequenced using an Illumina sequencing platform by Genedenovo Biotechnology Co., Ltd. (Guangzhou, China). BWA 0.7.1 was used to contrast high-quality clean reads with the reference genome [20]. The Glycine_max_v4.0 version of Williams 82 (https://www.ncbi.nlm.nih.gov/assembly/GCF_000004515.6/, accessed on 10 February 2023), was employed as a reference genome. SNPs and InDel variants were detected by utilizing GATK software v4.1 [21]. GWAS analysis involved 2777022 filtered high-quality SNPs spanning the entire genome with a minimum allele frequency (MAF) of ≥0.05. A mixed linear model (MLM) in GAPIT 3.0 was utilized to identify the association signals in relation to pod length and width [22]. The GWAS was performed with a modified Bonferroni correction (p < 1/2777,022), that is, the threshold of −log10 (P) was approximately 6.4 [23]. Pairwise linkage disequilibrium (LD) decay was determined using PopLDdecay software (Version 3.42), and the mean LD across all chromosomes decayed to r² = 0.10 within approximately 100 kb [24]. Manhattan and quantile–quantile (QQ) plots were generated using the “CMplot” package (version 4.4.1) in R environment [25]. If the physical distance between neighboring significant SNPs was less than 100 kb, the SNPs were grouped into a single QTL. The leading SNPs comprised those with the minimum p-value within the QTL.

2.3. RNA-Seq and Data Analysis

Four vegetable soybean varieties with extremely different pod lengths and widths were selected for subsequent RNA-seq analysis. Pods without seeds at the R5 stage of M3A00102412 (60.60, 12.02), GNY001 (19.69, 3.99), M3A00100063 (56.33, 13.84), and ZDD4867 (39.69, 9.82) were collected, with three biological replications, respectively. Total RNA was extracted using a TRIzol Kit (Invitrogen, Waltham, MA, USA). Sequencing libraries were constructed in accordance with the manufacturer’s instructions. Then, libraries were sequenced on the Illumina HiSeq2500 platform. GO enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed using OmicSmart platform (https://www.omicsmart.com, accessed on 19 December 2023).

2.4. RNA Extraction and Reverse Transcription Quantitative Real-Time PCR (RT-qPCR)

Around 0.1 g soybean tissue was ground into a powder in liquid nitrogen. A Total RNA Miniprep kit (Axygene, Shanghai, China) was used to extract total RNA following the manufacturer’s instruction. cDNA was synthesized from total RNA using a HiScript IV RT SuperMix for qPCR (Vazyme, Nanjing, China). RT-qPCR was performed using SsoAdvanced Universal SYBR Green Supermix (Bio-rad, Hercules, CA, USA) and an Applied Biosystems CFX Duet (Bio-rad, Hercules, CA, USA). The primers used for RT-qPCR are listed in Table S4. GmActin was used as an internal control. Gene expression level was calculated as described by Wang et al. [26].

2.5. Development of KASP Marker

KASP-polymerase chain reaction (PCR) amplification primers were designed using primer BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/ accessed on 11 October 2023) based on the SNP loci that showed a significant correlation with soybean pod length and width. Each pair of primers consisted of two specific forward primers, F1 and F2, and a generic reverse primer, R (Table S4). PCR amplification was performed as described by Jia et al. [27]. The genotype of the amplified product can be verified through fluorescence detection.

3. Results

3.1. Phenotype Description of Pod Size in the Association Panel

In this study, the phenotypic observations of the pod length and width exhibited continuity and a normal distribution. The mean pod length and width were calculated using the data collected across six different growth environments (Xuanwu in 2022, Lishui in 2022, Liuhe in 2022, Xuanwu in 2023, Lishui in 2023, and Liuhe in 2023). The results show that the mean pod length was 41.09 mm, and the mean pod width was 10.09 mm (Table S1). The coefficient of variation (CV) ranged from 13.75 to 17.21 (

Table 1. Statistical analysis of pod length and width in soybean accessions.

Traits	Site	Year	Mean ± SD	Kurtosis	Skewness	Min	Max	CV (%)
Pod length	Xuanwu	2022	41.93 ± 6.25	0.74	0.80	28.13	66.12	14.91
	Lishui		40.69 ± 6.75	0.41	0.71	26.12	62.05	16.59
	Liuhe		40.75 ± 6.29	1.58	0.62	23.82	67.84	15.44
	Xuanwu	2023	42.87 ± 6.11	1.38	0.30	18.50	62.97	14.25
	Lishui		42.70 ± 6.46	1.04	0.39	19.96	63.65	15.13
	Liuhe		42.74 ± 7.00	0.53	0.46	25.61	66.84	16.38
Pod width	Xuanwu	2022	10.14 ± 1.49	0.34	0.53	7.16	15.31	14.69
	Lishui		9.90 ± 1.56	0.03	0.05	5.20	14.29	15.76
	Liuhe		10.04 ± 1.38	0.33	0.38	6.78	15.29	13.75
	Xuanwu	2023	10.17 ± 1.62	−0.13	0.47	6.69	14.99	15.93
	Lishui		10.23 ± 1.65	−0.53	0.34	6.53	14.08	16.13
	Liuhe		10.11 ± 1.74	0.02	0.40	5.61	14.90	17.21

). The results indicate that the phenotype was relatively stable, and the traits of the pod length and width were regulated under genetic factors. These phenotypic data are suitable for GWAS analysis. The broad-sense heritability (h²) values of the pod length and width were 0.82 and 0.86, respectively, which indicate that genetic factors play a predominant role in determining the phenotypic expression of these traits (Table 1). Pearson’s correlation analysis revealed significant and strong correlations between the pod length measurements across various locations and years and the pod width measurements. No correlations were observed between the pod length and width (Figure 1).

3.2. Genome Resequencing and Population Structural Analysis

All 301 soybean accessions were resequenced at a sequencing depth of 12.09× [19]. After filtering the SNPs with an MAF greater than 0.05 and a missing rate less than 0.20, we obtained 2777022 SNPs for subsequent GWAS analysis (Figure 2). The average LD decay distance of the panel was approximately 100 kb, at which point r² dropped to 0.1. The phylogenetic analysis showed that the soybean population used in this study was divided into four subgroups (Figure S1), and combined with the results of the ADMIXTURE analysis, we selected k = 4 for the subsequent GWAS analysis.

3.3. Genome-Wide Association Studies

BLUE values from a two-year period were utilized for the GWAS analysis in this study. A GWAS was performed using an MLM method considering the genetic relatedness (K matrix). A total of 85 SNPs (−log₁₀P > 5.4) significantly correlated with the pod length and width were identified. Of these SNPs, 49 were significantly associated with the pod length, and 36 were significantly associated with the pod width (Figure 3 and Tables S2 and S3). These SNPs were located on chromosomes 6, 9, 13, and 17 (Figure 3), which can be combined into four QTLs. Two QTLs, namely, qGPoL1 and qGPoL2, are associated with the pod length. These loci are located on chromosomes 6 and 17, with leading SNP p-values of 1.65 × 10⁻⁶ and 4.4 × 10⁻⁸, respectively. The genomic regions of qGPoL1 and qGPoL2 span 846.08 and 257.24 kb, respectively (Figure 3A,C). Similarly, two QTLs (qGPoW1 and qGPoW2) associated with the pod width were identified on chromosomes 9 and 13, with leading SNP p-values of 1.01 × 10⁻⁶ and 1 × 10⁻⁸, respectively. The genomic region of qGPoW1 and qGPoW2 cover 215.02 and 200.04 kb, respectively (Figure 3B,D).

3.4. Analysis of Haplotype and Development of KASP Markers

To validate the efficacy of SNPs, we performed a haplotype analysis on the SNP sites strongly associated with the pod length and width. After observation, an allele variation of A/T was observed at the SNP site S06_42138365. The average pod length for S06_42138365-A was 46.99 mm, and that for S06_7373336-T was 37.50 mm (Figure 4A).

Another locus significantly associated with the pod length was S17_20494877, where the allelic variation was G/A. The average pod length for G was 39.67 mm, and it was 40.85 mm for A (Figure 4B). For the pod width, an allele variation of G/A was detected at SNP site S09_63148. The average pod width for S09_63148-G was 8.74 mm, and that for S09_63148-A was 9.85 mm (Figure 4C). Meanwhile, in another SNP site (S13_628331) with an A/T allele variation, the average pod width was measured as 9.29 mm for S13_628331-T and 12.59 mm for S13_628331-A (Figure 4D).

These results indicate that the S06_42138365 and S13_628331 loci are significantly associated with the pod length and width phenotype in vegetable soybeans, respectively. To enhance the utilization of these loci in production, we developed KASP markers for SNP sites S06_42138365 (A/T) and S13_628331 (A/T), which showed significant correlations with the pod length and width in vegetable soybean (Table S4). A total of 27 vegetable soybean germplasms were genotyped using the KASP-labeled primers designed for S06-42138365 (A/T) and S13_628331 (A/T). The two KASP markers showed significant differentiation between the various vegetable soybean lines (Figure 5).

3.5. Transcriptome Analysis

We selected pods from soybean accessions at the R5 stage that exhibited significant differences in pod length and width (Figure 6A,B, respectively). Subsequently, we performed a transcriptome analysis using three biological replicates for each variety (Table S5). In the short- and long-pod groups, 3018 differentially expressed genes (DEGs) were identified. In the long-pod group, 1602 genes were up-regulated, and 1416 genes were down-regulated compared with the short-pod group (Figure 6C). A total of 1815 DEGs were identified in the wide- and narrow-pod groups. In the wide-pod group, 991 genes were up-regulated, and 824 were down-regulated compared with the narrow-pod group (Figure 6D). These up- or down-regulated genes, including kinases and transcription factors, such as ubiquitin carboxy terminal hydrolase, BHLH domain containing protein, mitogen activated protein kinase kinase kinase 17, ubiquitin conjugating enzyme E2 and E3 ubiquitin protein ligase, and so on, are involved in multiple functions.

The results of the GO classification reveal that the DEGs were primarily enriched in molecular functions, followed by biological processes and cellular component terms (Figure 7A,B and

Table 2. Top 20 GO terms in Figure 7A,B.

Group	ID	Descrption
S vs. L	GO:0043531	ADP binding
	GO:0016762	xyloglucan:xyloglucosyl transferase activity
	GO:0048046	apoplast
	GO:0006073	cellular glucan metabolic process
	GO:0044042	glucan metabolic process
	GO:0015197	peptide transporter activity
	GO:0042887	amide transmembrane transporter activity
	GO:0016758	transferase activity, transferring hexosyl groups
	GO:0008194	UDP−glycosyltransferase activity
	GO:0003824	catalytic activity
	GO:1904680	peptide transmembrane transporter activity
	GO:0005618	cell wall
	GO:0046527	glucosyltransferase activity
	GO:0035251	UDP−glucosyltransferase activity
	GO:0016757	transferase activity, transferring glycosyl groups
	GO:0030246	carbohydrate binding
	GO:0015238	drug transmembrane transporter activity
	GO:0090484	drug transporter activity
	GO:0005260	channel-conductance-controlling ATPase activity
	GO:0001871	pattern binding
N vs. W	GO:0043531	ADP binding
	GO:0055114	oxidation−reduction process
	GO:0003824	catalytic activity
	GO:0016491	oxidoreductase activity
	GO:0004800	thyroxine 5′−deiodinase activity
	GO:0005506	iron ion binding
	GO:0019203	carbohydrate phosphatase activity
	GO:0016773	phosphotransferase activity, alcohol group as acceptor
	GO:0016301	kinase activity
	GO:0004672	protein kinase activity
	GO:0016772	transferase activity, transferring phosphorus−containing groups
	GO:0004805	trehalose−phosphatase activity
	GO:0044699	single-organism process
	GO:0005976	polysaccharide metabolic process
	GO:0046351	disaccharide biosynthetic process
	GO:0006468	protein phosphorylation
	GO:0044706	multicellular organism process
	GO:0036094	small molecule binding
	GO:0008509	anion transmembrane transporter activity
	GO:0015197	peptide transporter activity

). The KEGG classification indicated the significant enrichment of DEGs in pathways such as metabolism and environmental information processing (Figure 7C,D and

Table 3. Top 20 pathways in Figure 7C,D.

Group	ID	Descrption
S vs. L	ko00410	beta−Alanine metabolism
	ko04075	Plant hormone signal transduction
	ko00250	Alanine, aspartate and glutamate metabolism
	ko00909	Sesquiterpenoid and triterpenoid biosynthesis
	ko00903	Limonene and pinene degradation
	ko00620	Pyruvate metabolism
	ko00280	Valine, leucine and isoleucine degradation
	ko01110	Biosynthesis of secondary metabolites
	ko00650	Butanoate metabolism
	ko00340	Histidine metabolism
	ko00071	Fatty acid degradation
	ko00053	Ascorbate and aldarate metabolism
	ko00380	Tryptophan metabolism
	ko00592	alpha−Linolenic acid metabolism
	ko00430	Taurine and hypotaurine metabolism
	ko00010	Glycolysis/Gluconeogenesis
	ko00910	Nitrogen metabolism
	ko00561	Glycerolipid metabolism
	ko00310	Lysine degradation
	ko00591	Linoleic acid metabolism
N vs. W	ko00591	Linoleic acid metabolism
	ko01110	Biosynthesis of secondary metabolites
	ko00196	Photosynthesis—antenna proteins
	ko01100	Metabolic pathways
	ko04075	Plant hormone signal transduction
	ko00053	Ascorbate and aldarate metabolism
	ko04712	Circadian rhythm—plant
	ko00909	Sesquiterpenoid and triterpenoid biosynthesis
	ko02010	ABC transporters
	ko00250	Alanine, aspartate and glutamate metabolism
	ko00592	alpha−Linolenic acid metabolism
	ko00906	Carotenoid biosynthesis
	ko00910	Nitrogen metabolism
	ko00650	Butanoate metabolism
	ko00908	Zeatin biosynthesis
	ko00564	Glycerophospholipid metabolism
	ko00430	Taurine and hypotaurine metabolism
	ko00565	Ether lipid metabolism
	ko00500	Starch and sucrose metabolism
	ko00350	Tyrosine metabolism

).

3.6. Expression Pattern Analysis of Candidate Genes

Candidate gene screening and function prediction were performed within a range of 100 kb upstream and downstream of SNP sites significantly associated with the pod length and width of vegetable soybeans. Concerning the gene function annotation information on the soybean genome or the functional annotation of homologous genes in Arabidopsis thaliana, we identified 17 candidate genes significantly associated with the pod length and width (

Table 4. Candidate genes associated with pod length and width identified by GWAS.

Trait	Gene Id	Chr	Start (bp)	End (bp)	Description	Homologues in Arabidopsis	Symbols	Associated QTL	Range (bp)
Pod length	Glyma.06g254000	6	42,631,245	42,637,778	Arf-GAP domain-containing protein	AT5G54310.1	NEV, AGD5	qGPoL1	227,423
	Glyma.06g254200	6	42,644,161	42,649,221	Proteasome subunit alpha type	AT5G35590.1	PAA1
	Glyma.06g254400	6	42,705,376	42,706,251	Myb_DNA-bind_3 domain-containing protein	AT4G02210.1	MYB3
	Glyma.06g255000	6	42,821,205	42,830,287	Ubiquitin carboxyl-terminal hydrolase	AT3G20630.1	UBP14, TTN6, PER1
	Glyma.06g256800	6	43,298,448	43,299,529	heparanase-like protein 1	AT5G07830.1	GUS2
	Glyma.17g172400	17	17,415,802	17,419,290	BHLH domain-containing protein	AT4G36930.1	SPT	qGPoL2	4,080,227
	Glyma.17g173000	17	17,871,127	17,872,552	Mitogen-activated protein kinase kinase kinase 17	AT4G36950.1	MAPKKK21
	Glyma.17g173100	17	17,944,623	17,946,054	Transcription factor HEC2	AT3G50330.1	HEC2
	Glyma.17g175400	17	18,426,656	18,427,743	Uncharacterized protein	AT2G15680.1	/
	Glyma.17g175700	17	18,450,866	18,452,903	Transcription factor HEC2	AT1G49620.1	KRP7, ICN6, ICK5
Pod width	Glyma.09g001100	9	87,106	89,636	RING finger and U-box domain-containing protein isoform 1	AT2G44410.1	C3HC4_3	qGPoW1	215,024
	Glyma.09g001200	9	91,019	91,684	uncharacterized protein	AT5G61510.1	ADH
	Glyma.09g001500	9	117,912	120,188	Uncharacterized protein	AT3G60460.1	DUO1
	Glyma.09g002600	9	217,790	224,912	Ethylene receptor	AT1G66340.1	ETR1, EIN1, ETR
	Glyma.13g007000	13	2,061,623	2,062,337	Ubiquitin-conjugating enzyme E2	AT3G08690.1	UBC11	qGPoW2	200,041
	Glyma.13g008300	13	2,495,461	2,501,101	Uncharacterized protein	AT3G11220.2	PAXNEB
	Glyma.13g009100	13	2,637,797	2,639,471	E3 ubiquitin-protein ligase	AT2G04240.1	XERICO

).

We analyzed the expression levels of the 17 candidate genes, which were identified through GWAS analysis, at the R5 stage. We observed that the expression levels of Glyma.06G255000 and Glyma.17G173000 exhibited significant differences between the short- and long-pod groups (Figure 8A). Glyma.06G255000 encodes a ubiquitin carboxyl terminal hydrolase family protein, and Glyma.17G173000 encodes mitogen-activated protein kinase kinase kinase 17. Glyma.13G007000 encodes the ubiquitin conjugating enzyme E2. The expression level of Glyma.13G007000 in the wide-pod accessions was significantly higher than that in the narrow-pod accessions (Figure 8B).

We hypothesized that these three genes may be involved in the regulation of the pod length and width in vegetable soybeans. Consequently, we conducted further analysis on the expression levels of these genes in various soybean tissues [28]. The results show that Glyma.06G255000 was preferentially expressed in developing pods and seeds, and Glyma.13G007000 was mainly expressed in developing pods only (Figure 9A,B). However, the expression levels of Glyma.17G173000 were relatively higher in the roots than in other tissues (Figure 9C) but were lower in the pods and seeds, which suggests that they may play a relatively minor role in pod development.

4. Discussion

In soybean functional genomics research, a GWAS is often used for gene mapping. Through a GWAS analysis combined with multi-omics data analysis, a gene controlling soybean grain thickness and size, that is, GmST05 (Seed Thickness 05), was identified [29]. Dt2, the main control gene for soybean branch number, was identified via a GWAS analysis [30]. A major QTL for soybean plant height regulation, that is, PH13 on chromosome 13, was colocated through a GWAS analysis and transcriptome association analysis [31]. Liang et al. identified GmSW17 (Seed Width 17), which controls soybean grain width, through GWAS analysis [11]. These studies mainly focused on yield-related traits of soybean for grain use. However, attempts with vegetable soybean, especially its pod traits, have been limited. García-Fernández et al. (2021, 2023) identified more than 60 QTLs associated with pod morphology in common bean through GWAS and QTL mapping. They also revealed two major loci with epistatic effects located on chromosomes Pv01 and Pv06, which control pod morphology [14,32]. These findings indicate the complex genetic regulation of pod morphology.

The pod size is closely related to the yield of vegetable soybeans. Fresh soybean varieties with large pods must be screened to increase the vegetable soybean yield and improve the economic benefits. In this study, we performed a GWAS analysis on the soybean pod length and width and successfully identified four loci significantly associated with pod length and width traits in vegetable soybeans (Figure 1, Figure 2 and Figure 3). Through further transcriptome and RT-qPCR analyses, our findings suggest that Glyma.06g255000 is involved in the regulation of pod length and Glyma.13g007000 in the regulation of pod width (Figure 6, Figure 7, Figure 8 and Figure 9).

According to the annotation on the Soybase website (https://www.soybase.org/, accessed on 19 September 2024), Glyma.06g255000 encodes ubiquitin carboxyterminal hydrolase (Table 4). Our findings indicate that in groups of long-pod vegetable soybeans, the expression level of Glyma.06g255000 is lower than that in short-pod vegetable soybeans (Figure 8A). This finding suggests that Glyma.06g255000 plays a negative regulatory role in pod length. These findings are consistent with the results obtained in Arabidopsis, where the mutation of the homologous gene AtUBP14 resulted in organ enlargement [33], which further validated the accuracy of our results. Glyma.13g007000 is annotated as encoding ubiquitin conjugating enzyme (Table 4). The expression pattern analysis indicated that Glyma.13g007000 is highly expressed in wide-pod resources (Figure 8B), which suggests a positive regulatory role of the gene in pod width.

Glyma.06g255000 and Glyma.13g007000 participate in protein ubiquitination, which suggests that the ubiquitination pathway contributes to the regulation of the pod length and width of vegetable soybeans. Ubiquitination is a common post-translational modification. Protein ubiquitination plays a crucial role in the regulation of organ size [34,35]. In Arabidopsis, DA1 and DA2 regulate organ size through the ubiquitin pathway [35,36,37]. In rice, GW2 encodes an E3 ubiquitin ligase that negatively regulates grain width by promoting the ubiquitination of WG1 and GW9 [38,39,40]. Our findings provide a novel perspective for the comprehension of the regulatory network underlying soybean pod development. However, the molecular mechanism of pod length and width for vegetable soybeans remains unclear. Hence, more functional genes responsible for pod length must be cloned, and their regulatory mechanisms must be investigated.

Traditional crop breeding practices typically involve the observation of the traits of a large number of individuals within a population to select the superior offspring. This process is not only arduous and time consuming but also ineffective. The rapid advancement of molecular biotechnology has led to the emergence of molecular MAS as an effective tool for modern molecular breeding. KASP is a recently developed high-throughput genotyping technique. This process is mainly based on SNPs and has great potential for application in fields, such as crop trait improvement, due to its advantages of high throughput, low cost, and strong operability [16,17]. Currently, KASP markers are extensively utilized in molecular MAS breeding, genetic map construction, and seed purity identification in crops, such as rice, wheat, maize, and Chinese cabbage [41,42,43,44]. In this work, we developed KASP molecular markers based on significant SNPs. These markers can effectively distinguish different genotypes (Figure 5) and provide new reliable markers for molecular MAS breeding of soybeans. This development holds potential for increasing the efficiency and precision of soybean breeding programs.

However, further research is necessary to validate the functions of these genes and evaluate their efficacy in practical applications. Furthermore, we propose the completion of gene editing and other experiments to generate novel varieties of ideal vegetable soybeans. By undertaking these comprehensive studies, we anticipate making significant contributions to the genetic improvement and quality enhancement of fresh soybeans.

5. Conclusions

Our findings provide a robust toolkit for the molecular breeding of vegetable soybean. Through a GWAS analysis, transcriptome analysis, and other experimental methods, we identified Glyma.06G255000 and Glyma.13G007000 as possible regulators of pod length and width, respectively. Moreover, we have developed and validated KASP molecular markers for the screening of soybeans with different pod lengths and widths. These results deepen our understanding of the regulatory mechanisms underlying soybean pod development.