Combining Fine Mapping, Whole-Genome Re-Sequencing, and RNA-Seq Unravels Candidate Genes for a Soybean Mutant with Short Petioles and Weakened Pulvini

A short petiole is an important agronomic trait for the development of plant ideotypes with high yields. However, the genetic basis underlying this trait remains unclear. Here, we identified and characterized a novel soybean mutant with short petioles and weakened pulvini, designated as short petioles and weakened pulvini (spwp). Compared with the wild type (WT), the spwp mutant displayed shortened petioles, owing to the longitudinally decreased cell length, and exhibited a smaller pulvinus structure due to a reduction in motor cell proliferation and expansion. Genetic analysis showed that the phenotype of the spwp mutant was controlled by two recessive nuclear genes, named as spwp1 and spwp2. Using a map-based cloning strategy, the spwp1 locus was mapped in a 183 kb genomic region on chromosome 14 between markers S1413 and S1418, containing 15 annotated genes, whereas the spwp2 locus was mapped in a 195 kb genomic region on chromosome 11 between markers S1373 and S1385, containing 18 annotated genes. Based on the whole-genome re-sequencing and RNA-seq data, we identified two homologous genes, Glyma.11g230300 and Glyma.11g230600, as the most promising candidate genes for the spwp2 locus. In addition, the RNA-seq analysis revealed that the expression levels of genes involved in the cytokinin and auxin signaling transduction networks were altered in the spwp mutant compared with the WT. Our findings provide new gene resources for insights into the genetic mechanisms of petiole development and pulvinus establishment, as well as soybean ideotype breeding.


Introduction
The geometrical and topological organization of components of various plant types and shapes defines the architecture of plants [1]. For aboveground plant parts, this mainly includes stem growth habit, branching pattern, plant height, internode length, petiole length, and leaf size as well as shape in soybean. Among these components, petiole length is an important factor that influences canopy architecture that directly affects light interception efficiency, photosynthetic efficiency, and, ultimately, yield. The short petiole trait is potentially useful for improving the per unit yield by improving planting density and altering the canopy profile in soybean [2][3][4]. For instance, stable GmMYB14-overexpressing transgenic soybean plants displayed a compact plant architecture with short petioles that can be cultivated under higher density, thereby showing increased yields [5].
University, Jiangsu Province, China. For genetic analysis and gene mapping, the spwp mutant, as the male parent, was crossed with NN1138-2 and KF1 to construct different genetically segregating populations. 180 F 2 wild-type individuals of a KF1 × spwp cross were randomly harvested to obtain a F 2:3 line population for progeny tests and gene finemapping in the following year. NN1138-2, KF1, spwp, and their F 1 , F 2 , and F 2:3 (derived from the F 2 of a KF1 × spwp cross) populations were grown in the normal growing season at the Jiangpu Agricultural Experiment Station of Nanjing Agricultural University. From a line with a separation ratio of 3:1, a pair of near-isogenic lines, viz., the WT (wild type) and MT (mutant type), was constructed through selfing to advance the generations. The nearisogenic lines were planted in greenhouse conditions at Nanjing Agricultural University (Nanjing, China) to study their phenotypic characteristics.

Histological Analysis
The middle part of the pulvini and petioles at the third node, counted from the top to the bottom, were collected from the WT and the spwp mutant plants. All samples were fixed for one week at 4 • C in formalin-aceto-alcohol (FAA) fixing solution (50% alcohol, 5% glacial acetic acid, 5% formalin, and 5% glycerine). The tissues were prepared for paraffin sectioning following the previously described method [24]. Samples were sectioned with a microtome (Leica, RM2245) at a 2 µm thickness and then stained with 1% safranin O/0.5% Fast Green; sections were observed and photographed with a microscope (OLYMPUS BX53). Image J software was used to measure the petiole cell length.

Inheritance Analysis and Molecular Mapping with Simple Sequence Repeat (SSR) Markers
At the R3 stage (a reproductive stage with~5 mm long young pods at one of the four uppermost nodes on the main stem) the phenotypes of individuals in the F 2 and F 2:3 segregating populations were investigated by visual inspection, and the number of mutant-type and wild-type plants was then recorded. A Chi-square (χ 2 ) test was used to analyze the segregation ratio of alleles with the expected ratio at a significance threshold of p-value > 0.05 (χ 2 < 3.84).
Genomic DNA was extracted using a DNAquick Plant System Kit (TIANGEN, DP321) following the standard protocol. Bulk segregant analysis (BSA) was applied to rapidly identify SSR markers potentially linked to the mutated genes [26]. A total of 1015 SSR markers uniformly distributed on the 20 chromosomes were used in the BSA analysis. Among KF1 × spwp F 2 populations, an equal amount of DNA from each of 10 wild-type plants and 10 mutant-type plants formed a wild-type DNA pool and a mutant-type DNA pool, respectively. Primer sequences of SSR markers were obtained from the SoyBase website (http://www.soybase.org, accessed on 10 October 2021) or collected from Song et al. [27]. The PCR amplification program was 95 • C for 3 min followed by 34 cycles of 95 • C for 15 s, 56 • C for 15 s, and 72 • C for 1 min, with final incubation at 72 • C for 5 min before being held at 4 • C. The polymerase chain reaction (PCR) amplification products were separated on 8 % non-denaturing polyacrylamide gels that were stained with 1 g L −1 AgNO 3 for 15 min followed by 16 g L −1 NaOH plus 1.1 mL CH 3 OH for 10 min before being visualized under an LED light box. The primers used in the fine-mapping are listed in Tables S1 and S2.

DNA Library Construction and Whole-Genome Re-Sequencing
Genomic DNA was extracted from 5 plants of the WT and MT (a pair of near-isogenic lines) and then equally pooled to generate a WT DNA pool and an MT DNA pool, respectively. A total of 1.5 µg of genomic DNA for each pool was used to construct their individual sequencing library using a Truseq Nano DNA HT Sample Preparation Kit (Illumina, San Diego, CA, USA), following the manufacturer's recommendations. The final libraries constructed above were sequenced by the Illumina HiSeq4000 platform, and 150 bp paired-end reads were generated with an insert size of around 350 bp. To make sure reads were reliable and without artificial bias in the following analysis, raw reads with ≥10% unidentified nucleotides (N), with > 50% bases having phred quality < 5, or with > 10 nt aligned to the adapter, allowing ≤ 10% mismatches were removed to obtain clean reads. In addition, we trimmed putative PCR duplicates generated by PCR amplification in the library construction process. The clean reads were aligned to the soybean reference genome Williams 82 (Glycine max Wm82.a4. v1) using Burrows-Wheeler Aligner (BWA, v. 0.7.10; settings: mem-t 4-k 32-M-R) [28]. Alignment files were converted into BAM files using SAMtools software [29]. Variant calling was performed by using the Unified Genotyper function in GATK v. 3.8 software [30]. The single-nucleotide polymorphisms (SNPs) were obtained by using the Variant Filtration parameter in GATK (settings: -filterExpression "QD < 4.0 || FS > 60.0 || MQ < 40.0 ", -G_filter "GQ < 20", -cluster WindowSize 4). InDel was filtered via the Variant Filtration parameter (settings: -filter Expression "QD < 4.0 || FS > 200.0 ||Read PosRankSum < −20.0 || Inbreeding Coeff < −0.8 "). SnpEff 4.3t was then applied to annotate all the variants [31]. Theoretically, the genotype of the causal SNPs allele should be completely present in the MT DNA pool, that is to say no reads containing variant target SNP loci should be present in the WT DNA pool. In the present study, an average SNP index of the SNPs in each genomic interval was calculated using a sliding window analysis with a 2 Mb window size and a 50 kb increment. The delta SNP index = SNP index (MT) − SNP index (WT).

RNA Isolation and Quantitative Real-Time RT-PCR Analysis
Total RNAs from leaf petioles and pulvini were extracted using an RNAprep Pure Plant Kit (TIANGEN, Beijing, China) according to the manufacturer's instructions. Approximately 2 µg of RNA was then reverse-transcribed into first-strand cDNA using a HiScript II Q RT SuperMix for qPCR (+gDNA wiper) Kit (Vazyme, Nanjing, China), following a standard protocol. Quantitative real-time RT-PCR (qRT-PCR) was performed using a ChamQ SYBR qPCR Master Mix (Vazyme, China) on a Roche LightCycler480 PCR system, following the manufacturer's instructions. The PCR cycling conditions were 95 • C for 30 s followed by 40 cycles of 95 • C for 10 s and 60 • C for 30 s. Three biological replicates with three technical replicates were used for qRT-PCR assays. The soybean housekeeping gene GmActin11 (Glyma.18G290800) was used as an internal reference to normalize all data. The primers used for qRT-PCR were designed by National Center for Biotechnology Information database (NCBI) Primer BLAST and listed in Table S3.

RNA-Seq analysis
The near-isogenic lines (WT and MT) were grown in a greenhouse for RNA-seq analysis. The stem tips of the WT and MT were individually collected, placed in liquid nitrogen, and stored at −80 • C. Total RNA was isolated using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA integrity was assessed using an RNA Nano 6000 Assay Kit of a Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). The RNA samples were prepared and submitted to the Novogene Biotech Company (Beijing, China) for sequencing. Raw paired-end reads were first filtered with Fastp software to obtain high-quality clean data and then aligned to the Glycine max Wm82.a4. v1 using Hisat2 v2.0.5 software. Gene expression (FPKM, fragments per kilobase of transcript per million fragments mapped) levels were calculated based on the length of the gene and read count mapped to this gene. DEGs were identified using the criteria p-value ≤ 0.05 and |log2 (Fold Change)| ≥ 1. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DEGs was implemented by the clusterProfiler R package.

Phenotypic Characterization of the Spwp Mutant
The phenotypes of the wild-type (WT) and the spwp mutant-type (MT) plants in the greenhouse are shown in Figure 1A. Compared with the WT, the plant type of the spwp mutant is more compact and displays a significant decrease in the length of the petioles. Taking the petiole length from the first to the third leave (from the top to the bottom) as an example, the petiole length was reduced from (10.0 ± 2.0) cm, (13.6 ± 1.3) cm, and (14.3 ± 1.4) cm in the WT to (3.5 ± 0.8) cm, (3.9 ± 1.5) cm, and (3.5 ± 0.8) cm in the MT plants, respectively ( Figure 1B). To elucidate the mechanism of short petiole formation at the cytological level, longitudinal sections of petioles were performed. In the MT plants, the cell length was significantly shorter than in the WT ones ( Figure 1C-E). This result suggests that the decrease in cell elongation is mainly attributed to the shortened petioles in mutant plants.

Phenotypic Characterization of the Spwp Mutant
The phenotypes of the wild-type (WT) and the spwp mutant-type (MT) plants in the greenhouse are shown in Figure 1A. Compared with the WT, the plant type of the spwp mutant is more compact and displays a significant decrease in the length of the petioles. Taking the petiole length from the first to the third leave (from the top to the bottom) as an example, the petiole length was reduced from (10.0 ± 2.0) cm, (13.6 ± 1.3) cm, and (14.3 ± 1.4) cm in the WT to (3.5 ± 0.8) cm, (3.9 ± 1.5) cm, and (3.5 ± 0.8) cm in the MT plants, respectively ( Figure 1B). To elucidate the mechanism of short petiole formation at the cytological level, longitudinal sections of petioles were performed. In the MT plants, the cell length was significantly shorter than in the WT ones ( Figure 1C-E). This result suggests that the decrease in cell elongation is mainly attributed to the shortened petioles in mutant plants. The pulvinus of the MT was significantly smaller than that of the WT, and the MT did not have the knee-like structure at the base of short petioles that was present in the WT (Figure 2A,B). By comparing the anatomical structure of pulvini ( Figure 2C,D) we found that there were fewer and smaller motor cells on both the abaxial side (AB) and adaxial side (AD) in the MT than in the WT ( Figure 2E,F,I,J). However, the xylem area of The pulvinus of the MT was significantly smaller than that of the WT, and the MT did not have the knee-like structure at the base of short petioles that was present in the WT (Figure 2A,B). By comparing the anatomical structure of pulvini ( Figure 2C,D) we found that there were fewer and smaller motor cells on both the abaxial side (AB) and adaxial side (AD) in the MT than in the WT ( Figure 2E,F,I,J). However, the xylem area of the MT plants was larger than that of the WT, both on the AD and AB ( Figure 2G,H,K,L). These observations indicated that cell expansion and proliferation were also affected in the pulvini of the spwp mutant plants compared with the WT.
Genes 2022, 13, x FOR PEER REVIEW 6 of 20 the MT plants was larger than that of the WT, both on the AD and AB ( Figure 2G,H,K,L). These observations indicated that cell expansion and proliferation were also affected in the pulvini of the spwp mutant plants compared with the WT. The leaf petiole angle is mainly controlled by the pulvinus [32,33]. To investigate whether the petiole angle of the MT was affected, we measured the petiole angle from the third leaf position to the fifth leaf position (from the top to the bottom) and found that the petiole angle of the MT was slightly increased in comparison to the WT, but the difference was not significant ( Figure S1). Apart from the differential phenotypes mentioned above, the nyctinastic motions of pulvini and leaflets were significantly diminished in the spwp mutant compared to the wild-type plants.

Genetic Analysis of the Spwp Mutant
To analyze the inheritance pattern of the spwp mutant, we crossed spwp with NN1138-2 and KF1, respectively. All F1 plants of crosses NN1138-2 × spwp and KF1 × spwp showed a similar phenotype to the wild type, indicating that the spwp mutant is recessive. In the F2 population of cross KF1 × spwp, among the 692 plants only 42 showed the spwp mutant phenotype. The ratio of wild-type plants relative to mutant-type plants corresponded to the expected 15:1 segregation ratio for two recessive genes (χ 2 = 0.02, p = 0.87) ( Table 1). However, in the case of the NN1138-2 × spwp cross, the phenotypic segregation The leaf petiole angle is mainly controlled by the pulvinus [32,33]. To investigate whether the petiole angle of the MT was affected, we measured the petiole angle from the third leaf position to the fifth leaf position (from the top to the bottom) and found that the petiole angle of the MT was slightly increased in comparison to the WT, but the difference was not significant ( Figure S1). Apart from the differential phenotypes mentioned above, the nyctinastic motions of pulvini and leaflets were significantly diminished in the spwp mutant compared to the wild-type plants.

Genetic Analysis of the Spwp Mutant
To analyze the inheritance pattern of the spwp mutant, we crossed spwp with NN1138-2 and KF1, respectively. All F 1 plants of crosses NN1138-2 × spwp and KF1 × spwp showed a similar phenotype to the wild type, indicating that the spwp mutant is recessive. In the F 2 population of cross KF1 × spwp, among the 692 plants only 42 showed the spwp mutant phenotype. The ratio of wild-type plants relative to mutant-type plants corresponded to the expected 15:1 segregation ratio for two recessive genes (χ 2 = 0.02, p = 0.87) ( Table 1). However, in the case of the NN1138-2 × spwp cross, the phenotypic segregation ratio of wild-type relative to mutant-type plants was 420:120, which did not fit the expected ratio of 15:1 but did fit a ratio of 3:1 (χ 2 = 2.22, p = 0.14) ( Table 1). We randomly selected 180 F 2 -dominant individuals derived from the KF1 × spwp cross to perform the progeny test. In the F 2:3 lines, the ratio of segregating to non-segregating lines was 93:87, which fitted an expected ratio of 8:7 (χ 2 = 0.2, p = 0.65). These results indicate that the short petioles and weakened pulvini traits of the spwp mutant are controlled by two recessive genes, tentatively designated as spwp1 (short petioles and weakened pulvini 1) and spwp2 (short petioles and weakened pulvini 2). These results also imply that the NN1138-2 genome contains either of those two homozygous recessive genes or alleles.

Fine Mapping of the Spwp1 and Spwp2 Genes
To locate the target genes, a total of 1015 SSR markers were used to screen for polymorphisms between the wild-type and mutant-type DNA pools derived from the F 2 population of the KF1 × spwp cross. We found that three SSR markers, i.e., Satt063, Satt560, and Satt726, on chromosome 14 (linkage group B2), and another three SSR markers, i.e., Satt359, Sat_331, and BE801538, on chromosome 11 (linkage group B1), showed polymorphisms. All 42 F 2 -mutant individuals were used for the further identification of linkage groups and preliminary mapping. New SSR markers from Song et al. [27] were synthesized to determine the mapping regions. After genotyping, the spwp1 gene was mapped between the SSR markers Satt063 and S1465 on chromosome 14 ( Figure 3A); the spwp2 gene was mapped between the SSR markers S1335 and Sat_331 on chromosome 11 ( Figure 4A). ratio of wild-type relative to mutant-type plants was 420:120, which did not fit the expected ratio of 15:1 but did fit a ratio of 3:1 (χ 2 = 2.22, p = 0.14) ( Table 1). We randomly selected 180 F2-dominant individuals derived from the KF1 × spwp cross to perform the progeny test. In the F2:3 lines, the ratio of segregating to non-segregating lines was 93:87, which fitted an expected ratio of 8:7 (χ 2 = 0.2, p = 0.65). These results indicate that the short petioles and weakened pulvini traits of the spwp mutant are controlled by two recessive genes, tentatively designated as spwp1 (short petioles and weakened pulvini 1) and spwp2 (short petioles and weakened pulvini 2). These results also imply that the NN1138-2 genome contains either of those two homozygous recessive genes or alleles.

Fine Mapping of the Spwp1 and Spwp2 Genes
To locate the target genes, a total of 1015 SSR markers were used to screen for polymorphisms between the wild-type and mutant-type DNA pools derived from the F2 population of the KF1 × spwp cross. We found that three SSR markers, i.e., Satt063, Satt560, and Satt726, on chromosome 14 (linkage group B2), and another three SSR markers, i.e., Satt359, Sat_331, and BE801538, on chromosome 11 (linkage group B1), showed polymorphisms. All 42 F2-mutant individuals were used for the further identification of linkage groups and preliminary mapping. New SSR markers from Song et al. [27] were synthesized to determine the mapping regions. After genotyping, the spwp1 gene was mapped between the SSR markers Satt063 and S1465 on chromosome 14 ( Figure 3A); the spwp2 gene was mapped between the SSR markers S1335 and Sat_331 on chromosome 11 ( Figure  4A). : the spwp1 locus was further mapped to a region between SSR markers S1382 and S1465 using theF2 population of NN1138-2 × spwp. (C): the position of the spwp1 locus was finally narrowed down to a genomic region of 183 kb between SSR markers S1413 and S1418 using the  : the spwp1 locus was further mapped to a region between SSR markers S1382 and S1465 using theF 2 population of NN1138-2 × spwp. (C): the position of the spwp1 locus was finally narrowed down to a genomic region of 183 kb between SSR markers S1413 and S1418 using the The genetic analysis of the F 2 population of the NN1138-2 × spwp cross indicated that only one gene controls the phenotype of the spwp mutant under the genetic background of NN1138-2. To test which marker is linked to the mutant phenotype in the NN1138-2 × spwp population, all 111 F 2 -mutant individuals were genotyped using Sat_331 and Satt726 markers. The results showed that the SSR marker Satt726 on chromosome 14 was detected to be linked with the mutant phenotype, and then the spwp1 locus was further delimited to a 1.44 Mb region between SSR markers S1382 and S1465 ( Figure 3B). To isolate the spwp1 and spwp2 genes, 573 F 2:3 -mutant individuals derived from the heterozygous F 2 plants of the KF1 × spwp cross were used for further fine-mapping. Under the polymorphic markers S1413, S1417, and S1418, two, zero, and seven recombination events were detected in the 573 F 2:3 -mutant individuals, respectively, indicating that the spwp1 locus was pinpointed to the segment with a physical distance of approximately 183 kb between SSR markers S1413 and S1418 on chromosome 14 ( Figure 3C). Additionally, under the polymorphic markers S1373, S1381, and S1385, one, zero, and three recombination events were detected, respectively, suggesting that the spwp2 locus was mapped to a physical distance of an approximately 195 kb region on chromosome 11 between SSR markers S1373 and S1385 ( Figure 4B).  The genetic analysis of the F2 population of the NN1138-2 × spwp cross indicated that only one gene controls the phenotype of the spwp mutant under the genetic background of NN1138-2. To test which marker is linked to the mutant phenotype in the NN1138-2 × spwp population, all 111 F2-mutant individuals were genotyped using Sat_331 and Satt726 markers. The results showed that the SSR marker Satt726 on chromosome 14 was detected to be linked with the mutant phenotype, and then the spwp1 locus was further delimited to a 1.44 Mb region between SSR markers S1382 and S1465 ( Figure 3B). To isolate the spwp1 and spwp2 genes, 573 F2:3-mutant individuals derived from the heterozygous F2 plants of the KF1 × spwp cross were used for further fine-mapping. Under the polymorphic markers S1413, S1417, and S1418, two, zero, and seven recombination events were detected in the 573 F2:3-mutant individuals, respectively, indicating that the spwp1 locus was pinpointed to the segment with a physical distance of approximately 183 kb between SSR markers S1413 and S1418 on chromosome 14 ( Figure 3C). Additionally, under the polymorphic markers S1373, S1381, and S1385, one, zero, and three recombination events were detected, respectively, suggesting that the spwp2 locus was mapped to a physical distance of an approximately 195 kb region on chromosome 11 between SSR markers S1373 and S1385 ( Figure 4B).

Whole-Genome Re-Sequencing Analysis
To rapidly identify the causal mutations, next-generation sequencing (NGS) was applied. Two DNA pool libraries were constructed and subsequently subjected to highthroughput whole-genome re-sequencing using the Illumina HiSeq 4000 platform. After filtering 102.49 G of raw data, 101.64 G of clean data was obtained for further analysis. The average Q20 ratio was ∼ 98% and the Q30 ratio was ∼94%, indicating the high quality of the sequencing data (Table S4). The numbers of clean reads were 163,984,765 for the wild type, with an average >31× and genome coverage >90.63%, and 168,202,480 for the mutant type, with an average >31× and genome coverage >90.62% (Table S5). A total of 2,046,994 SNPs and insertion-deletions (InDels) were obtained between the two re-sequenced samples for further analysis.
Because the near isogenic lines were selected from the offspring of a heterozygous single plant with a separation ratio of 3:1 (wild-type plants to mutant-type plants), only one pair of differential genes controlling short petioles and weakened pulvini traits exists

Whole-Genome Re-Sequencing Analysis
To rapidly identify the causal mutations, next-generation sequencing (NGS) was applied. Two DNA pool libraries were constructed and subsequently subjected to highthroughput whole-genome re-sequencing using the Illumina HiSeq 4000 platform. After filtering 102.49 G of raw data, 101.64 G of clean data was obtained for further analysis. The average Q20 ratio was~98% and the Q30 ratio was~94%, indicating the high quality of the sequencing data (Table S4). The numbers of clean reads were 163,984,765 for the wild type, with an average >31× and genome coverage >90.63%, and 168,202,480 for the mutant type, with an average >31× and genome coverage >90.62% (Table S5). A total of 2,046,994 SNPs and insertion-deletions (InDels) were obtained between the two re-sequenced samples for further analysis.
Because the near isogenic lines were selected from the offspring of a heterozygous single plant with a separation ratio of 3:1 (wild-type plants to mutant-type plants), only one pair of differential genes controlling short petioles and weakened pulvini traits exists in them, and the other gene is homozygous, i.e., the genotypes of the near isogenic lines are SPWP1 SPWP1 spwp2 spwp2 vs. spwp1 spwp1 spwp2 spwp2 or spwp1 spwp1 SPWP2 SPWP2 vs. spwp1 spwp1 spwp2 spwp2. To determine the genotypes of near isogenic lines, we used the SNP index method to analyze the association for high-quality SNPs and InDels between the MT and WT DNA pools. The delta (SNP index) between the MT and WT across a 2 Mb window size was measured using a 50 kb step size and plotted for all 20 chromosomes of the soybean genome ( Figure 5). Through analyzing the delta (SNP index) value, only one major peak was identified on chromosome 11 (Gm11), which suggests that the main differential loci of the near isogenic lines were located on chromosome 11, and therefore their genotypes should be spwp1 spwp1 SPWP2 SPWP2 vs. spwp1 spwp1 spwp2 spwp2. The delta (SNP index) plotting, with a threshold of one, revealed that the interval with a physical distance of approximately 1.58 Mb on chromosome 11 between 36.26 Mb and 37.94 Mb was possibly linked to the phenotype of the spwp mutant. This positioning interval is consistent with the results of the linkage analysis, suggesting that the accuracy of the positioning results by the whole-genome re-sequencing analysis.
WT across a 2 Mb window size was measured using a 50 kb step size and plotted for all 20 chromosomes of the soybean genome ( Figure 5). Through analyzing the delta (SNP index) value, only one major peak was identified on chromosome 11 (Gm11), which suggests that the main differential loci of the near isogenic lines were located on chromosome 11, and therefore their genotypes should be spwp1 spwp1 SPWP2 SPWP2 vs. spwp1 spwp1 spwp2 spwp2. The delta (SNP index) plotting, with a threshold of one, revealed that the interval with a physical distance of approximately 1.58 Mb on chromosome 11 between 36.26 Mb and 37.94 Mb was possibly linked to the phenotype of the spwp mutant. This positioning interval is consistent with the results of the linkage analysis, suggesting that the accuracy of the positioning results by the whole-genome re-sequencing analysis.

RNA-Seq Analysis between the Spwp Mutant and WT
As another approach toward identifying the candidate genes and understanding the molecular mechanism underlying the spwp mutant's phenotype, RNA-seq analysis was performed to investigate gene expression changes in the spwp mutant compared with the wild type. After filtering raw reads, an average of 41.53 million and 42.14 million clean reads were obtained from three WT samples and three MT samples, respectively. The percentages of Q30 were all above 93%, and the percentage of uniquely mapped reads ranged from 93.33% to 94.72%, indicating that the sequencing results were of high quality (Table  S6).
Compared with the WT, a total of 1436 differentially expressed genes (DEGs) were detected in the spwp mutant, including 963 up-regulated DEGs and 473 down-regulated DEGs ( Figure 6A; Table S7). Gene Ontology (GO) analysis revealed that the DEGs were highly enriched for biological processes associated with "response to biotic stimulus", "phosphorelay signal transduction system", and "multi-organism process", and were also enriched for molecular function categories related to "ADP binding", "peroxidase activity", and "polysaccharide binding" (Figure 6B). In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment of DEGs showed that the "protein processing in

RNA-Seq Analysis between the Spwp Mutant and WT
As another approach toward identifying the candidate genes and understanding the molecular mechanism underlying the spwp mutant's phenotype, RNA-seq analysis was performed to investigate gene expression changes in the spwp mutant compared with the wild type. After filtering raw reads, an average of 41.53 million and 42.14 million clean reads were obtained from three WT samples and three MT samples, respectively. The percentages of Q30 were all above 93%, and the percentage of uniquely mapped reads ranged from 93.33% to 94.72%, indicating that the sequencing results were of high quality (Table S6).
Compared with the WT, a total of 1436 differentially expressed genes (DEGs) were detected in the spwp mutant, including 963 up-regulated DEGs and 473 down-regulated DEGs ( Figure 6A; Table S7). Gene Ontology (GO) analysis revealed that the DEGs were highly enriched for biological processes associated with "response to biotic stimulus", "phosphorelay signal transduction system", and "multi-organism process", and were also enriched for molecular function categories related to "ADP binding", "peroxidase activity", and "polysaccharide binding" (Figure 6B). In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment of DEGs showed that the "protein processing in endoplasmic reticulum", "zeatin biosynthesis", "brassinosteroid (BR) biosynthesis", and "isoflavonoid biosynthesis" pathways were highly enriched ( Figure S2).

Candidate Genes Analysis by Integrating Whole-Genome Re-Sequencing and RNA-Seq Data
Based on the annotations of the Glycine max Wm82.a4.v1 genome database from SoyBase (http://www.soybase.org, accessed on 10 October 2021), 15/18 putative genes with annotated functions are predicted to reside in the spwp1/spwp2 locus, respectively ( Figure 3D, Table 2; Figure 4C, Table 3). The 18 genes located in the spwp2 locus were screened for sequence variations in the coding DNA sequence between the WT and MT pool genome using the re-sequencing data. Compared with the WT pool genome, twelve SNPs were identified in the coding region of six candidate genes. Of these SNPs, ten formed missense substitutions in the amino acid sequence of five candidate genes (Table 4). These genes included: CHX15 (Glyma.11G229400), a member of the sodium/hydrogen exchanger family; Glyma.11g230000, a member of the phosphoinositide-specific phospholipase C family; a transcription factor, NF-X1 (Glyma.11G230200); and two homologous genes (Glyma.11G230300 and Glyma.11G230600) that encode protein-kinase-domain-containing proteins.  Based on the synteny information obtained from SoyBase, the mapping region of spwp1 and spwp2 belongs to two nonhomologous genomic regions. Furthermore, according to the protein homolog information on Phytozome (https://phytozome-next.jgi.doe.gov/, accessed on 10 October 2021), no duplicated gene pairs were found between the spwp1 locus and spwp2 locus. Hence, spwp1 and spwp2 are two nonhomologous genes that likely have similar functions in regulating the development of petioles and pulvini. To further screen the candidate genes we analyzed the expression levels of all 33 predicted genes (15 genes from the spwp1 locus and 18 genes from the spwp2 locus) using the RNA-seq data. Only two DEGs (Glyma.11g230300 and Glyma.11g230600) in the spwp2 locus were significantly up-regulated in the spwp mutant ( Figure 7A,B). It is interesting to note that Glyma.11g230300 and Glyma.11g230600 are homologous genes that encode protein-kinase-domain-containing proteins and that the functional annotation of the Glyma.14g206000 gene located in the spwp1 locus was also a "protein kinase domain". To validate the RNA-seq analysis results we selected six genes (five genes having the SNPs mentioned above plus Glyma.14g206000) and performed a qRT-PCR analysis. The results of the qRT-PCR analysis showed the same expression trends as the RNA-seq data ( Figure 7C). From the re-sequencing data of the MT pool, three mutation points were identified in the Glyma.11g230300 gene, which all caused amino acid substitution. Moreover, one synonymous mutation and two missense point mutations were detected in the Glyma.11g230600 gene (Table 3). To verify the results of the re-sequencing data and identify other sequence variations in the Glyma.11g230300 and Glyma.11g230600 genes, the coding sequences of these two candidate genes in the wild-type and mutant-type plants were amplified and sequenced. In addition to the SNPs identified by the NGS, two synonymous mutations and one missense point mutation were also detected in both genes, which led to a substitution of Ser with Arg in the Glyma.11g230300 gene and a substitution of Ile with Thr in the Glyma.11g230600 gene (Table S8). These results suggested that Glyma.11g230300 and Glyma.11g230600 may be the most important candidate genes for the spwp2 locus, which is responsible for the short petioles and weakened pulvini phenotype of the spwp mutant in soybean.

Enrichment Analysis on DEGs Associated with the Phenotype of the Spwp Mutant
Notably, the GO analysis indicated that the DEGs were highly enriched in the GO term "phosphorelay signal transduction system" ( Figure 6B). Among which, three homologous genes (Glyma.03G130000, Glyma.11G155100, and Glyma.04G137600) encoding response regulator receiver-domain-containing proteins homologous to Arabidopsis response regulator 9 (ARR9), two homologous genes (Glyma.06G187000 and

Enrichment Analysis on DEGs Associated with the Phenotype of the Spwp Mutant
Notably, the GO analysis indicated that the DEGs were highly enriched in the GO term "phosphorelay signal transduction system" ( Figure 6B). Among which, three homologous genes (Glyma.03G130000, Glyma.11G155100, and Glyma.04G137600) encoding response regulator receiver-domain-containing proteins homologous to Arabidopsis response regulator 9 (ARR9), two homologous genes (Glyma.06G187000 and Glyma.04G177900) encoding a homologous protein of ARR6, and Glyma.08G292400, encoding a homologous protein of ARR3, were all upregulated. The ARRs together with histidine protein kinases such as CKI1 (cytokinin-independent 1) or CRE1 (cytokinin response 1) constitute two-component systems that play a central role in cytokinin signal transduction [34]. To et al. found that ARR5 and ARR6 function additively with ARR3 and ARR4 in the regulation of petiole elongation through the study of ARR multiple insertional mutants [35]. Moreover, the KEGG analysis showed that the DEGs were highly enriched in the "zeatin biosynthesis" pathway ( Figure  S2). The Glyma.17G054500, Glyma.09G225400, and Glyma.03G133300 genes were all upregulated in the spwp mutant, which are homologous to cytokinin oxidase 3 (CKX3), cytokinin oxidase 6 (CKX6), and cytokinin oxidase 1 (CKX1) in Arabidopsis, respectively. The CKX genes are important regulators of active cytokinin levels, which catalyze the degradation of cytokinin [36]. These results showed that both the signal transduction and metabolism of cytokinin were altered in the spwp mutant.
The plant hormone auxin controls diverse aspects of plant growth and development by regulating the fundamental cellular processes of expansion, division, and differentiation [37]. In the present study, we identified eighteen DEGs that are possibly involved in auxin signaling from the RNA-seq data (Table S9). Among them, two auxin efflux carrier genes, PIN3A (Glyma.07G217900) and PIN3B (Glyma.20G014300), were downregulated, while three genes (Glyma.01G114000, Glyma.03G063600, and Glyma.03G063900) encoding auxin-transporter-like proteins were upregulated. Except for polar transport, auxin signal transduction was also affected. The Glyma.15G091000 gene encoding auxin response factor 8 (ARF8) was upregulated. Three Aux/IAA transcription factor genes (Glyma.03G158700, Glyma.20G210400, and Glyma.13G356600) were upregulated. Besides, the transcript levels of five SAUR genes were changed in the spwp mutant compared with WT.
In addition, we examined changes in the expression levels of transcription factors (TFs) based on the RNA-seq data. We identified 151 differentially expressed TFs belonging to 16 families, mainly including MYB (35), ZF (27), bHLH (21), AP2/EREBP (19), and WRKY (18) ( Figure S3). Most of these families have previously been involved in plant growth and development. For example, overexpressing the AtMYB96 gene can cause a dwarf phenotype with reduced lateral roots in Arabidopsis [38]. Two paralogous AtMYB124 (FOUR LIPS, FLP) and AtMYB88 proteins are required to limit cell divisions in the stomatal lineage [39]. In our study, we found that 26 genes encoding MYB TFs were up-regulated in the spwp mutant, while only nine genes were down-regulated. There are also 21 bHLH TF genes that are differently expressed in the spwp mutant compared with the WT. bHLH transcription factors play key roles in phytochrome signal transduction, organ development, and BR-responsive gene expression [40][41][42]. Recently, a study showed that AtLP1 (leafrelated protein 1) and AtLP2, two bHLH homologous proteins, regulate longitudinal cell elongation in Arabidopsis [43].

Spwp Is a New Soybean Short Petiole Mutant with Weakened Pulvini
Short petioles are potentially useful for improving the per unit yield by altering the canopy profile and increasing planting density in soybean. Until now, many mutants with short petioles have been identified in soybean. For example, a previous study reported that the short petiole trait of D76-1609 was controlled by a single recessive gene, lps [2]. The short petiole trait in SS98206SP was also controlled by a single recessive gene designated as lps3, which was mapped on chromosome 13 between SSR markers Sat_234 and Sct_033 [4,44]. The short petiole mutant VP5 was identified from a fast neutron-induced mutant popula-tion, and the phenotype of the VP5 mutant was co-segregated with an 837,919 bp deletion fragment on chromosome 17 [45]. Moreover, the derived short petiole (dsp) mutant was controlled by two recessive genes, dsp1 and dsp2, which were mapped on chromosome seven with flanking markers BARCSOYSSR_07_0787 and BARCSOYSSR_07_0808 and on chromosome 11 between markers BARCSOYSSR_11_0037 and BARCSOYSSR_11_0043, respectively [46]. Recently, Wang et al. [47] found that the rolled leaves and short petioles mutant (rlsp1) was controlled by multiple genes, and identified 10 candidate regions on chromosomes three, six, eight, thirteen, and seventeen using sequencing-based bulked segregant analysis. The Gmilpa1 mutant has shorter petioles and smaller pulvini compared to the wild type, and the underlying gene was isolated through map-based cloning approaches. It encodes an anaphase-promoting complex-like protein that appears to function by promoting cell growth and division [24]. Thus, the genetic regulation of petiole length is quite complex, and more efforts are needed to better understand the genetic basis of petiole length using new short-petiole genetic materials. In this study, we found that the spwp mutant was controlled by two recessive genes, spwp1 and spwp2, which were, respectively, mapped in a 183 kb genomic region on chromosome 14 and in a 195 kb genomic region on chromosome 11 by using the traditional map-based cloning method (Figures 3A-D and 4A-C). To our knowledge, no genes related to the length of leaf petioles have been reported in those two regions. Hence, spwp is a novel short-petiole mutant controlled by two genes.
Pulvini show nyctinastic movement through turgor pressure changes in their motor cells [15]. In the present study, we found that the pulvini of the spwp mutant were incomplete and that the nyctinastic motions of pulvini and leaflets were both diminished in the spwp mutant. Furthermore, the paraffin section observation results showed that the number and size of motor cells on the abaxial and adaxial sites were reduced in the pulvini of the spwp mutant ( Figure 2C,D). Therefore, the defects in nyctinastic mobility in the spwp mutant were likely due to the under-developed pulvini. Although these results are similar to the previous characterization of the soybean Gmilpa1 mutant [24], the developmental defects in the pulvini are caused by different genes. Therefore, this study provides a new genetic resource for soybean ideotype breeding and understanding pulvinus development.

Identification Candidate Genes of the Spwp Mutant through Whole-Genome Re-Sequencing and RNA-Seq
The phenotype of the spwp mutant was controlled by two genes, spwp1 and spwp2, which were mapped to two nonhomologous fragments on chromosome 14 and chromosome 11, respectively. Fifteen genes were predicated within the spwp1 locus. Of these genes, five homologous genes contain the NB-ARC domain, and a recent study found that Rsc4-multiple aspects of plant growth and development. For example, ARABIDOPSIS CRINKLY 4 (ACR4) is essential for epidermal cell differentiation in leaves and seed coats [50], cell layer organization during ovule integument and sepal margin development [51], and proper embryogenesis [52]. OsCR4 is required to maintain the interlocking of the palea and lemma [53], and positively regulates culm elongation in rice [54]. Moreover, ACR4 is part of a mechanism controlling formative cell divisions in the Arabidopsis root [55,56]. Some studies have suggested that the ACR4 functions as a coreceptor of other receptors, such as CLV3 INSENSITIVE KINASES (CIKs) and CLAVATA1 (CLV1), to mediate CLE40 signaling, and that the CLE40-ACR4-WOX5 signaling pathway plays a central role in maintaining the stem cell niche and in controlling columella cell development [57][58][59][60]. It has been proven that WOX5 is involved in cytokinin signaling and auxin signaling pathways [61][62][63].
In the present study, our RNA-seq data showed that the expression levels of DEGs related to cytokinin signal transduction and metabolism were altered in the spwp mutant compared with the WT ( Figure 6B; Table S2). Moreover, a subset of DEGs involved in the auxin signaling transduction network were also detected in the spwp mutant (Table S9). The results of the paraffin section show that cell expansion and cell division were affected in the leaf petioles and pulvini of the spwp mutant ( Figures 1C,D and 2C,D). These results are reminiscent of the candidate genes Glyma.11g230300 and Glyma.11g230600, which may be involved in cytokinin and auxin signal transduction pathways to regulate cell division and cell expansion in the leaf petioles and pulvini of the spwp mutant. However, we cannot rule out the possibility of other genes in the positioning interval responsible for the spwp mutant phenotype, and more work needs to be done in the future to validate the candidate genes. It is possible to dissect the reasons of spwp mutant formation using the proteomics and metabolomics strategies as well as integrate transcriptomics with them to elaborate on the regulatory networks of petiole and pulvinus growth and development in different levels.

The Spwp Mutant Also Showed Defective Leaf Movement
In legume species, a common phenomenon is that leaves open during the day and fold at night, which is diurnal leaf movement named nyctinasty. Structurally, nyctinastic movement is driven by the pulvinus, a specialized motor organ located at the base of leaves and leaflets [21]. In this study, we found that the spwp mutant also showed defective leaf movement, that is, leaflets remained in a horizontal (open) position during day and night. Several factors have been reported to influence the leaf movement. For example, the F-box protein MIO1/SLB1, a component of the SKP1/Cullin/F-box (SCF) E3 ubiquitin ligase complex, influences leaf movement through possibly controlling the length of the pulvinus [64]. Aquaporins, anion channels, and K + channels are thought to regulate volume changes of the motor cells and thus control leaf movement [18,65,66]. Interestingly, a recent study showed that the geometry of the compound leaf is important for leaf movement in Medicago truncatula [67]. Moreover, the plant hormone BR homeostasis is critical for nyctinastic leaf movement but has no effect on pulvinus determination [68]. In this study, the KEGG pathway enrichment analysis showed that the DEGs were highly enriched in BR biosynthesis ( Figure S2). Among them, the transcript levels of Glyma.08G193900 and Glyma.14G059900 were up-regulated in the spwp mutant. The Glyma.14G059900 gene is homologous to the Arabidopsis gene AT5G05690, which encodes a cytochrome P450 monooxygenase that converts 6-deoxocathasterone into 6-deoxoteasterone in the late C6 oxidation pathway of brassinolide biosynthesis [69]. By contrast, the transcript level of Glyma.17G118100 was down-regulated, which is homologous to AT2G01190 (PHYB ACTIVATION TAGGED SUPPRESSOR 1, BAS1) in Arabidopsis. BAS1 is thought to be a control point between multiple photoreceptor signal transduction pathways and BR signaling [70,71]. In addition, GmBRI1b (Glyma.04g218300) was down-regulated in the spwp mutant. GmBRI1b, a homolog of AtBRI1, has been shown to function as a BR receptor [72]. These results suggested that not only BR biosynthesis but also BR signal transduction were altered in the spwp mutant, which was likely responsible for the defective leaf movement of the spwp mutant.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/genes13020185/s1, Figure S1: Comparison of the leaf petiole angle between the wild-type plants (WT) and the mutant-type plants (MT); Figure S2: The scatter diagram of KEGG enrichment; Figure S3: The expression analysis of differentially expressed transcription factors (TFs) based on the RNA-seq data; Table S1: The primers used for mapping the spwp1 locus; Table S2: The primers used for mapping the spwp2 locus; Table S3: List of the primers used for the qRT-PCR analysis of six candidate genes; Table S4: Summary of re-sequencing data quality of the wild-type and mutant-type DNA pools; Table S5: Mapping information for the wild-type and mutant-type DNA pools by whole-genome re-sequencing; Table S6: Throughput and quality of RNA-seq of the spwp mutant and the wild type (WT); Table S7: The differentially expressed genes between the spwp mutant and the wild type (WT); Table S8: The SNPs in the candidate genes detected by Sanger sequencing; Table S9: Auxin-related genes differentially expressed between the spwp mutant and the wild type (WT).

Data Availability Statement:
The data presented in this study are available in the article and supplementary material.

Conflicts of Interest:
The authors declare no conflict of interest.