Pfam ID (PF02234) searches and BLAST using A. thaliana KRP2 protein sequences (NP_190632) as queries were performed in the Glycine max reference genome (a4.v1 version) to identify members of the soybean KRP gene family; all hits were classified as candidate GmKRPs. The exon/intron organizations of the genomic sequences, CDS, and protein sequences for all KRPs were retrieved from Phytozome (V13, https://phytozome.jgi.doe.gov). The gene structure schematic diagram of the KRPs was drawn and visualized with TBtools (Chen et al., 2020).

Sequence alignments of multiple KRPs from Glycine max and other 6 dicots (Phaseolus vulgaris, Vigna unguiculata, Arabidopsis thaliana, Lotus japonicus, Medicago truncatula, and Cicer arietinum) and 3 monocots (Oryza sativa, Sorghum bicolor, and Zea mays) were performed using MUSCLE. The unrooted phylogenetic tree was made by the neighbor-joining method using MEGA 11 (Tamura et al., 2021). The reliability of an inferred tree was confirmed with bootstrap analysis performed with 1000 replications. The non-synonymous (Ka) and synonymous substitution (Ks) rate ratios of the paralog pairs were calculated using the TBtools program (Chen et al., 2020). The divergence time (T) was calculated following T = Ks/(2 x 6.1 x 10⁻⁹) x 10⁻⁶ Mya (Lynch and Conery, 2000).

The protein lengths, molecular weights, and isoelectric points of the GmKRPs were analyzed using the ProtParam online tool (Gasteiger et al., 2005). Protein subcellular localization was predicted using the WoLF PSORT online tool (Horton et al., 2007). All KRP protein sequences were analyzed to identify conserved protein motifs (Motif scan) using the ‘Multiple EM for Motif Elicitation’ (MEME) program (Bailey et al., 2009). The parameters were as follows: zero or one occurrence per sequence, motif width ranges of 6–50 amino acids, and 10 as the maximum number of finding motifs. Each of the sequences had an E-value of less than 10. The PEST domains were searched (http://emboss.bioinformatics.nl/cgi-bin/emboss/epestfind) (threshold score >+5.0) (Rogers et al., 1986). All SNP datasets located in the coding regions were extracted from whole genome re-sequencing datasets using SNPViz (https://soykb.org/SNPViz2/) as described by Langewisch et al. (2014). Non-synonymous SNPs (coding missense variant) were further analyzed and listed.

The full-length coding region of GmKRP2 was amplified from the root tissue of soybean ‘Williams82’ using gene-specific primers (GmKRP2-F: 5′-caccATGGAGATGGCTCAGGTTAAG-3′ and GmKRP2-R: 5′-TGGCTTCAATTTAACCCACTCG-3′) and Phusion high-fidelity DNA polymerase (Cat#: F530S, Thermo Scientific, Waltham, MA, USA). PCR products with the size of 556 bp were cloned into the Gateway pENTR-D-TOPO vector (Cat#: K240020, Invitrogen, Carlsbad, CA, USA). After sequence verification, it was recombined into the destination vector pMDC83 (Curtis and Grossniklaus, 2003) via LR reactions (Cat#: 11791020, Invitrogen, Carlsbad, CA, USA), which contained the 35S promoter of Cauliflower mosaic virus and fused with GFP at the N-terminal.

pMDC83-GmKRP2a::GFP was introduced into the Agrobacterium rhizogenes strain K599. Hairy roots infected with K599 containing the modified empty vector pMDC83 without the chloramphenicol resistance gene served as the control. Soybean transformation was performed using a previously reported protocol (Govindarajulu et al., 2008). Williams82 was used for hairy root transformation. After positive root selection using hygromycin (25 mg/L) in the medium, the transgenic roots were transferred to the new medium to determine the relative root growth length. At least 30 independent transgenic roots were used as biological replicates.

The pMDC83-GmKRP2a::GFP vector was introduced into Agrobacterium tumefaciens strain GV3101 and then transformed into wild-type Columbia-0 (Col-0) using the floral dip method (Clough and Bent, 1998). Putative transgenic lines were selected with 30 mg L⁻¹ hygromycin on a half-strength MS medium that contained 1% sucrose and 0.8% agar. Three single-copy transgenic lines from the second generation were selected according to the segregation ratio. In this study, seedlings of Col-0 and transgenic plants were grown on half-strength MS medium under a long-day photoperiod (16 h light, 8 h dark), and the temperature was 22°C with 50–70% relative humidity. The hypocotyl and root lengths were measured at 5 and 10 days, respectively.

For the visualization of subcellular protein in Arabidopsis roots, 5-day-old transgenic seedlings harboring GFP fusion protein were visualized. The root tip was observed using a Leica SP8 spectral confocal microscope. For the cell size analysis, Arabidopsis roots were stained with 1 µg/ml propidium iodide. For the visualization of subcellular protein in N. benthamiana leaves, the coding sequence of GmKRP2a was fused to the C-terminal under the control of the CaMV 35S promoter in the transient expression vector pHB-mCherry, referred to as GmKRP2a-mCherry. The transient expression construct was introduced into Agrobacterium strain GV3101 and then used to infiltrate N. benthamiana leaves, as described previously (Lim et al., 2014). For GmKRP2a-GFP transgenic seedlings nuclear staining, the seedling was transferred to 4′,6-diamidino-2-phenylindole (DAPI, 5 mg/L) stain solution and incubated for 1 to 2 min at room temperature. For GmKRP2a-mCherry in tobacco leaves, the leaves were transferred to DAPI (5 mg/L) stain solution and incubated for 20 to 30 min at room temperature. After incubation, seedlings or leaves were placed on slides and then imaged with fluorescence microscope Confocal imaging analysis was performed using a Zeiss LSM 880 NLO laser scanning confocal microscope.

Young Arabidopsis rosette leaves were chopped in LB01 nuclear dissociation buffer. The mixtures were filtered with a 40-μm cell filter, then the released nuclei were centrifuged (1000 rcf, 1 min) and resuspended in LB01 buffer. DAPI (final concentration of 4 μg/ml) was used to stain the nuclei before analysis with a flow cytometer (FACS LSRFortessa, BD, USA). For hairy roots, 2–3 cm of each root tip was cut in LB01 buffer to prepare the nuclear dissociation solution for flow cytometric detection and analysis of the content of nuclear DNA. The results were analyzed using Modfit 5.0 software. At least 1x10⁵ nuclei were used for the analysis.

Total RNA was isolated from 100 mg hairy roots with the plant RNA Purification reagent (Cat#:74904, Invitrogen, UAS) following the manufacturer’s protocols. The RNA quantity was checked with NanoDrop2000 (Thermo Scientific, USA). RNA integrity was analyzed using an Agilent 2100 Bioanalyzer (RNA Nano Chip, Agilent, Santa Clara, CA, USA). High-quality RNA (RIN ≥ 9, OD260/280 ≥ 1.8–2.2, OD260/230 ≥ 1.0, 28S:18S ≥ 1.0, > 1 μg)) was used for library construction. The TruSeq RNA sample prep kit (Cat#: RS-122-2101, Illumina, USA) was used to prepare the RNA-seq library, and the library was quantified with QuantiFluor^® dsDNA System. The paired-end RNA-seq sequencing library was sequenced with the Illumina NovaSeq6000 sequencer (2×150 bp read length, Illumina, Inc. San Diego, CA, USA).

The data were analyzed on the online platform of the Majorbio Cloud Platform (www.majorbio.com). Briefly, Fastp (https://github.com/OpenGene/fastp) was used to obtain high-quality clean reads with default parameters. Clean reads were then separately aligned to the reference genome (Glycine max Wm82.a4.v1), with orientation mode using HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtml) software (Kim et al., 2015). The mapped reads of each sample were assembled using StringTie (https://ccb.jhu.edu/software/stringtie/) in a reference-based approach (Pertea et al., 2015). KEGG pathway analysis was carried out using Goatools (https://github.com/tanghaibao/Goatools) and KOBAS (http://kobas.cbi.pku.edu.cn/home.do) (Xie et al., 2011).

Soybean Williams 82 seeds were germinated on a Petri dish lined with moist filter paper and then seedlings were transferred to half-strength MS solution grown in a growth chamber under a 10-h photoperiod at 25°C/22°C (day/night) temperature and 50 % relative humidity. The plants at the vegetative 1 stage were transferred to half MS solution containing 15% PEG6000 or 150 mM NaCl for 24 h and 48 h, respectively. The roots and first trifoliolate leaves from 5 plants were harvested and used for the GmKRP gene expression analysis. Samples were frozen in liquid nitrogen immediately after harvest and kept at -80°C for total RNA isolation. Soybean seeds were collected from the field.

Total RNA was isolated from soybean hairy roots or Arabidopsis seedlings or soybean seeds using an RNApure Plant Kit (CW0588S, CWBIO, Beijing, China), and cDNA was synthesized using HiScript III qRT SuperMix for qPCR (Cat#: R323-01, Vazyme, China). Gene-specific primers for quantitative PCR and the size of PCR products are listed in Supplementary Table 1 . Real-time PCR amplification was carried out with the CFX connect Real-Time PCR detection system (Bio-Rad) using a ChamQ Universal SYBR qPCR Master Mix (Q711-02). Actin11 (Glyma18g290800) was used as the housekeeping gene to normalize the expression level. All qPCR analyses had three biological replicates and two technical replicates. Differences in variables were analyzed by t-test. The differences were considered to be significant at p < 0.05 or p < 0.01.

To identify soybean KRP members, BLAST and PFAM (PF02234: cyclin-dependent kinase inhibitor) searches were performed using Arabidopsis KRP genes as guides. Nine putative GmKRP loci were found by searching the Glycine max reference genome (Wm82.a4.v1). The gene name, gene ID, protein information, and putative subcellular location are listed in Table 1 . Among the nine putative GmKRP proteins, the amino acid sequences of GmKRPs were relatively consistent in size (ranging from 168 to 224 amino acids), but the PI (isoelectric point) showed wide variation ranging from 5.44 to 9.37. The subcellular location of this gene family was analyzed to better understand its biological function. All GmKRP proteins were predicted to localize in the nucleus (Cell-PLoc 2.0 and CELLO), providing GmKRPs a possibility to play roles in the cell cycle. In addition, all GmKRPs were hydrophilic proteins, and their instability coefficient was over 40, indicating that the KRP family of proteins is unstable.

To better assess the diversity of the KRP gene family beyond soybean, an unrooted phylogenetic tree was constructed using the Neighbor-joining method based on multiple KRP protein sequence alignments from dicot and monocot plants ( Figure 1A ). Nine soybean KRP genes were named according to the cluster of the KRP family, and the others were named based on their chromosomal location or previous reports ( Supplementary Table 2 ). Accordingly, the multiple KRP sequences were divided into three major classes based on the phylogenetic tree. Class I contained 26 sequences from both monocotyledonous and dicotyledonous plants, including three soybean KRPs (GmKRP3, GmKRP4 and GmKRP5); Class II contained 16 sequences, only from monocotyledonous plants; Class III contained the remaining 24 sequences, specifically from dicotyledonous plants, including 6 soybean KRP members (GmKRP1a, GmKRP1b, GmKRP2a, GmKRP2b, GmKRP6 and GmKRP7). In addition, the KRP genes in the seven dicotyledons had closer relationships than those of the monocots, which was consistent with the evolutionary relationship between dicots and monocots. Furthermore, the KRPs from legumes (Glycine max, Phaseolus vulgaris, Vigna unguiculata, Lotus japonicus, Medicago truncatula, and Cicer arietinum) had closer relationships than those of genes from Cruciferae (Arabidopsis thaliana). These results prompt a hypothesis of the divergent evolution of KRP members in dicots and monocots.

The exons and introns of the GmKRP genes were further analyzed by aligning the genomic sequence with the corresponding protein coding sequence. As illustrated in Supplementary Figure 1 , three GmKRP genes had three exons (GmKRP3, GmKRP 4, GmKRP 5), while the rest contained four exons. The KRP members in other species contained three or four exons, except OsKRP2, OsKRP6, and ZmKRP8, which contained seven or two exons ( Supplementary Figure 1 ). Moreover, the number and position of exons were conserved in classes I and III, although the UTR region varied. The last exon in all genes was the most conserved in terms of size, sequence, and encoding protein. In general, the gene structure showed sequence relatedness consistent with the phylogenetic tree based on the protein sequences.

Nine soybean KRP family members were distributed on different chromosomes ( Figure 1B ). We detected 12 fragment duplicated gene couples among GmKRP genes in the soybean genome. No tandem duplications were found in the GmKRP gene family. Therefore, segmental duplication was the main reason for the expansion of the GmKRP genes. The Ka/Ks ratio was further calculated to evaluate the selection pressure during evolution. Determining whether positive Darwinian selection was involved in GmKRP divergence following duplication and estimating the date of the duplication pairs, as shown in Supplementary Table 3 , the ratio between different gene pairs varied from 0.178 to 0.413, indicating that the GmKRP genes primarily evolved under the influence of purifying selection. The duplication time among classes I and II varied from 140–257 Mya (million years ago). The duplication time of gene pairs in class I or class II was 68–83 Mya. However, the separation times of GmKRP1a/GmKRP1b, GmKRP6/GmKRP7, and GmKRP3/GmKRP5 pairs varied from 10 to 14 Mya; this period was consistent with the latest twice whole genome duplication of soybean. These results uncovered that GmKRP expansion was derived from whole genome duplication, resulting in conserved domains and motifs.

The alignment of the soybean KRP protein sequence showed that all GmKRP proteins had a highly conserved CID domain ( Supplementary Figure 2A ). However, the sequence composition at the N terminal showed wide variation. Moreover, 10 conserved sequence motifs were further identified using the MEME program, with an E-value less than 0.05 in the GmKRP amino acid sequence ( Supplementary Figure 1 ). Motifs 1 and 2, located at the extremity of the C-terminal part of the proteins, were shared by all examining KRPs, including GmKRPs. The distribution of other motifs follows phylogenetic classification. The GmKRPs in class I contained motifs 3, 4, 5, 6, 8, and 9. The GmKRPs in class II contained motif 9 or 10. The GmKRPs in class III contained motifs 5, 6, 7, and 10. The sequence information for each motif is provided in Supplementary Figure 2B . In addition to the conserved motifs, PEST domains, which were reported to serves as a proteolytic signal, were found in GmKRP1a, 1b, 2b, and 6 (Rogers et al., 1986, Supplementary Figure 2B ).

To further investigate the putative functions of the GmKRP genes, we analyzed the 2 kb promoter region of all GmKRP genes using PlantCARE web tools. As expected, Core cis-elements, TATA- and CAAT-box, were both found in all promoters in the promoter region. In addition, 12 types of cis-acting elements were found in the promoter region of the GmKRP genes ( Figure 2A ). These regulatory elements were mainly involved in abiotic stress responses (low temperature, anaerobic and drought) and hormone responses (methyl jasmonate, auxin, gibberellin, abscisic acid and salicylic acid). Among all the cis-elements, the light-responsive elements correspond to the highest frequency, and the soybean KRP gene is likely to participate in the stress response and light-induced development process in soybean growth. To investigate whether the expression of GmKRP genes was affected by abiotic stress, the expression levels of GmKRP genes were evaluated by qPCR analysis under abiotic stress (PEG and salt stress). Under PEG stress, the expression levels of GmKRP1a, GmKRP2b, GmKRP3, and GmKRP7 were downregulated in leaves, while three genes involving GmKRP1a, GmKRP2a, and GmKRP4 expression levels showed upregulated in root ( Figures 2B, C ). Under salt stress, several GmKRP genes were upregulated both in leaves and roots ( Figures 2D, E ).

A total of 862 co-expressed genes of the GmKRP family were retrieved from the soybean genome database. These co-expressed genes were used in underlying pathways that GmKRP may be involved in through Gene Ontology (GO) enrichment analysis. These co-expressed genes were mainly enriched in protein modification and phosphorylation, as well as Ras-genes (GO:0007265) and ARF, Adenosine diphosphate-ribosylation factor (GO: 0032011) signal transduction in the biological process ( Supplementary Figure 3A ). For molecular function terms, nucleotide binding, kinase activity, and kinase inhibitor activity were enriched in co-expressing genes. These results are consistent with the signal pathways involved in the regulation of the cell cycle by KRP and its molecular function prediction ( Supplementary Figure 3B ).

To gain insight into the tissue expression pattern of the GmKRP members, a soybean RNA-seq database was extracted, including eight tissues (root, root tip, lateral root, stem, leaf, shoot tip, open and unopen flower) (Libault et al., 2010). As illustrated in Figure 3A , GmKRP2b, GmKRP3, GmKRP4, GmKRP5, and GmKRP7 showed relatively uniform expression in all tissues. GmKRP4 showed relatively high expression in all tissues. Obviously, GmKRP1a, GmKRP1b, GmKRP2a, and GmKRP6 were identified as having tissue-specific expression. GmKRP1a and GmKRP1b was highly expressed in the unopened flower. GmKRP2a was highly expressed in stems. GmKRP6 showed predominantly expressed in the stem, shoot and unopened flower. The expression patterns of GmKRPs in five seed-development stages were further investigated using qPCR ( Figure 3B ). We found that GmKRP2a and GmKRP4 were highly expressed in early seed development stage ( Figure 3C ). In addition, the expression level of GmKRP2a was lower than that of other GmKRP genes at the late stage of seed development ( Figure 3D-G ). Overall, the tissue organ-specific expression pattern indicates that GmKRP transcripts are abundantly expressed in unopened flower and relatively weakly expressed in root.

The SNPs located in the GmKRP coding regions that may lead to function variation were retrieved from 775 soybean germplasms containing both wild and cultivated species using SNPViz software. A total of 52 non-synonymous substitution SNPs were observed in nine KRP genes. Non-synonymous substitution SNPs were absent in the protein-encoding region of GmKRP1b in any of the 775 soybean lines. Compared with the reference genome (Wm82.a4.v1), one and two non-synonymous substitution SNPs were found in GmKRP2a and GmKRP3 in one or two wild soybean germplasms, respectively. Three to six non-synonymous substitution SNPs were found in the coding region of GmKRP4, GmKRP5, and GmKRP7 in less than 20 germplasm resources. However, in GmKRP1a, GmKRP2b, and GmKRP6, more non-synonymous substitution SNPs (8–19 SNPs) were found, which come from more than 100 germplasm resources ( Supplementary Table 4 ). These results indicate that GmKRP1b and GmKRP2a may play a highly conserved function during evolution because few non-synonymous substitution SNPs were found in soybean landraces and elite lines.

To further study the function of the GmKRP genes, one of the highly conservative gene, GmKRP2a gene was cloned from Glycine max Wm82. GmKRP2a was overexpressed in the soybean hairy root via Agrobacterium rhizogenes-mediated transient transgenic system to elucidate the role of GmKRP2a during root development. After initially selecting the positive transgenic roots with hygromycin resistance, the same length transgenic roots of the control and overexpressing lines were transferred to the new medium to investigate the growth rate. The presence of the construct and its expression level were checked using PCR and qPCR, respectively ( Supplementary Figure 4 ). As shown in Figures 4A, B , the transgenic soybean hairy roots exhibited a 35–50% length increase relative to the controls after 7 days. To further investigate gene function, the same construct was used to transformed into Arabidopsis to generate stable transgenic lines. As shown in Figures 4C, E , three independent transgenic Arabidopsis lines with single insert copies were used for further phenotype analysis. An average root growth promotion of 29–34% and 24–37% hypocotyl promotion under dark conditions were observed in GmKRP2a overexpression lines ( Figures 4D, F ). These results indicate that GmKRP2a is a positive regulator of soybean root elongation.

In addition, the subcellular localization of GmKRP2a was investigated using the transient expression of 35S::GmKRP2a–mCherry in tobacco leaves. The results showed clear nuclear fluorescence in the tobacco cells ( Figures 5D, F ) compared with empty vector control ( Figures 5A–C ). The nuclei were stained with DAPI, which allows detection of nuclear DNA. Colocalization of mCherry (red) with the DAPI (blue) was shown by red arrows ( Figures 5G, H ). Furthermore, the subcellular localization of GmKRP2a was further validated using transgenic Arabidopsis carrying 35S::GmKRP2a–GFP. The overlap between DAPI and GFP fluorescence confirmed that GmKRP2a was a nuclear protein, which is consistent with the localization prediction and the function in cell cycle ( Supplementary Figure 5 ).

To address the gene function of GmKRP2a in the cell cycle, the ploidy level in rosette leaves of GmKRP2a-Over and Col-0 plants was determined using flow cytometric analysis. The proportion of cells with an increased DNA content in GmKRP2a-Over transgenic plants was higher than that of Col-0 ( Figures 6A–C ). Moreover, the maximal nuclear DNA level reached 16C in the GmKRP2a-overexpressing roots, whereas the cells of the WT roots displayed an essentially 2C and 4C DNA content distribution. The endoreduplication level was increased in these GmKRP2a-overexpressing plants. Furthermore, the distribution of 2C nuclei was shifted proportionally to the corresponding 4C, 8C, and 16C nuclei. GmKRP2a-overexpressing hairy roots showed an overall higher endoreduplication than the control roots ( Figures 6D, E , P < 0.01), which strongly suggests that GmKRP2a inhibits the cell cycle during the G2/M phase.

Since we previously found that the root elongation was promoted, whereas the cell division was inhibited in GmKRP2a-overexpresion plants, the question on whether GmKRP2a influences cell expansion were further addressed. Therefore, the roots from the GmKRP2a-overexpressing transgenic lines and the WT were observed using a fluorescence microscope after propidium iodide staining. The cell size of the cortex cells in the proximal meristem was more promoted than that of Col-0 ( Figures 7A–D ). The cell size was determined using ImageJ software. The overall increase in size was observed in the GmKRP2a-overexpressing Arabidopsis lines, which induced 30–50% cell expansion ( Figure 7E ).

High-quality RNA-seq was performed to investigate the transcriptomes of 35S::GmKRP2a hairy root lines and the empty vector control for the genome-wide discovery of differentially expressed genes (DEGs) ( Supplementary Table 5 ). A total of 1291 DEGs were identified with the criteria |log2FC|≥1 and p<0.05 ( Supplementary Table 6 ). Among these, 989 and 302 DEGs were up- and downregulated, respectively. A subset of 18 DEGs was selected for the validation of the expressions using qPCR, and the results showed that 18 DEGs displayed expression patterns consistent with the expression levels obtained from RNA-seq analysis ( Supplementary Figure 6 ). Interestingly, there were 102 genes encoding eukaryotic small subunit ribosomal RNA among the downregulated DEGs. Moreover, no SSU_rRNA were upregulated in 35S::GmKRP2a hairy root. In addition, three 5.8S ribosomal RNA were significantly downregulated. These results may indicate that GmKRP2a overexpression decreases the number or activity of ribosomes by downregulating the gene expression of ribosomal small subunit RNA, which further results in a decrease in the level of protein translation in cells and a stagnation of the cell cycle in the G2/M phase.

Furthermore, all functionally annotated transcript sequences were selected and assigned to the reference canonical pathways in KEGG. A total of 127 genes were distributed to 14 different pathways in the KEGG database (p<0.05). KEGG annotation analysis revealed that the plant–pathogen interaction pathway was the most enriched KEGG term with the highest level of gene representation (48 transcripts, Supplementary Table 7 ), followed by lipid metabolism, energy metabolism, and signal transduction ( Supplementary Figure 7 ). This result suggests that increased endoreduplication and inhibited cell division may correlate with pathogen resistance in plants.

As a crucial negative regulator of the cell cycle, the KRP genes have been identified, and a considerable amount of research have been studied in several plant species and have been demonstrated to play significant roles on plant development and in response to various stressors (Polyn et al., 2015; Kumar and Larkin, 2017; Vieira and de Almeida Engler, 2017; Xiao et al., 2017; Sahu et al., 2021). In this study, a comprehensive analysis of genes encoding KRP proteins in the soybean genome was performed, resulting in the identification of nine KRP family members. The phylogenetic analysis and evolutionary process of the KRP family laid the foundation for further functional investigation of KRP genes. Among the 10 species analyzed in this paper, soybean contains the largest number of KRP genes. The soybean genome underwent two rounds of whole-genome duplication, which occurred 59 and 13 million years ago (Schmutz et al., 2010). The expansion of GmKRP genes has arisen more recently due to soybean-specific duplication. This result may indicate that GmKRPs play important roles in regulating biological processes, which need to be demonstrated further. As expected, GmKRPs were grouped more closely to legume KRPs compared to Arabidopsis and monocotyledonous plants.

The gene structures were highly conservative in all 10 species. Plant KRPs shared a low sequence identity with non-plant KRPs. The CID domain at the C terminal, which is required for interacting with the CDKA protein and inhibiting its activity in plants, was the most conserved motif among all plant KRP proteins, (Wang et al., 1998; Zhou et al., 2003b; Torres Acosta et al., 2011). Sequence alignment of the KRP family in this study confirmed the presence of the conserved C-terminal domain in soybean, whereas the randomness of the N-terminal sequence was much higher than that of the C-terminal. The N-terminal region increases KRP1 instability, and the central domain is responsible for nuclear localization (Schnittger et al., 2003; Zhou et al., 2003a). It was considered that the diversity of the N-terminal amino acids caused a diversity of different KRP functions.

The identification of SNPs and their potential effects or functional annotations in one gene family is crucial for studying the genetic basis of phenotypic differences (Schacherer et al., 2009). Several gene-based SNPs have been reported for their roles in controlling traits, such as grain filling, plant height, amylose content, and grain structure in rice and maize (Tanaka et al., 1995; Monna et al., 2002; Wang et al., 2008; Godínez-Palma et al., 2017). In this study, only one non-synonymous substitution of GmKRP2a was identified in one Glycine soja line (PI 562534) across 775 soybean germplasm resources, which resulted in amino acid conversion from Ala¹⁵¹ to Thr¹⁵¹. This amino acid conversion may affect the gene function because it is in the conserved cyclin-binding domain. Sequence variants that affect gene functions can be used to design functional molecular markers (Andersen and Lubberstedt, 2003). Therefore, this identified SNP can be further used to develop functional markers.

The initiation and division of root or leaf primordium or the formation and activation of the meristem during plant tissue development is largely determined by the transition between cell division and cell expansion during the cell cycle (Inze and De Veylder, 2006; Inagaki and Umeda, 2011). These processes and transitions are regulated by a protein complex composed of highly conserved regulators of CYC, CDK, and CKI. However, plant growth involves the integration of many environmental and endogenous signals in which hormones act as an internal signal involved in the determination of plant size. GA, auxin, and brassinosteroids affect proportional changes in organ growth rates (Himanen et al., 2002; Verkest et al., 2005b; Achard et al., 2009; González-Garcia et al., 2011; Zhao et al., 2017; Li et al., 2022). In this study, the promoter region of GmKRP2a contained many methyl jasmonate, auxin, gibberellin, and abscisic acid response elements. Overexpression of GmKRP2a appeared to restrain cell division and enhanced the levels of endoreduplication. Therefore, it could be hypothesized that one of the pathways that hormone control the cell cycle by modulating the expression levels of GmKRP2a.

Drought is one of the main unfavorable factors affecting the growth and development, seed yield, and quality of soybean. The annual loss of soybean yield and quality caused by drought is immeasurable. The trait of deeper roots allows the plant to search for additional soil water while maintaining higher water-use efficiency under drought stress conditions (Calleja-Cabrera et al., 2020). Here, overexpression of GmKRP2a in both soybean transgenic hairy root and stably transformed Arabidopsis plants leads to a significant increase in root elongation, which provides an experimental basis for how to use GmKRP2a to promote the development of roots and improve the drought resistance/tolerance of soybean. Understanding the role of cell cycle regulators during plant development will not only clarify the molecular mechanism, but also guide the agricultural practice. In the future, how GmKRP2a connects root development with environmental signals to feed into cell cycle regulation and better targeted in specific cells should be addressed.

However, the role of KRP is very complicated, and its function in different tissues may be different. For example, overexpression of NtKIS1a in Arabidopsis induces modified plant morphology with smaller organs that contain larger cells, similar to Arabidopsis plants overexpressing KRP genes (Wang et al., 1998; De Veylder et al., 2001; Zhou et al., 2002). NtKIS1a and NtKIS2 showed obvious differences in expression patterns during both the cell cycle and plant development (Jasinski et al., 2002; Jasinski et al., 2003). Overexpressed Orysa; KRP1 reduced cell production during leaf development and dramatically reduced seed filling (Barroco et al., 2006). Overexpression of OsiCKI6 reduces the seed-setting rate, partially due to the poorly developed pollen, and results in a reduced cell number, but larger cells in the leaf and stem coincide with the transgenic plant dwarfed phenotype (Yang et al., 2011). In the current study, no obvious phenotype was observed in GmKRP2a-overexpressing Arabidopsis thaliana aboveground tissue (data not shown).

It has been reported that KRP3, KRP5, KRP6, and KRP7 play roles in the root knot nematodes resistance by inhibiting the size of giant cells or affecting nematode development and offspring (Vieira et al., 2012; Vieira et al., 2014; Coelho et al., 2017). Hamdoun et al. (2016) reported that SIM and SMR1 act in innate immunity in Arabidopsis. These results suggest that precise cell cycle control is not only critical for plant development but also important for plant–pathogen response. In this study, overexpressed GmKRP2a induced the expression level of many plant–pathogen-related genes. Among them, BAK1 (Glyma.05G119500) and FLS2 (Glyma.05G128200) were both induced in the GmKRP2a overexpressing lines. The BAK1 and FLS2 complex acts as a receptor for several LRR-type PRRs (Chinchilla et al., 2007; Heese et al., 2007; Wang et al., 2020). Therefore, we propose that GmKRP2a may play a role in the plant defense response. Further investigation and characterization of the GmKRP2a function will provide a better understanding of crosstalk between cell cycle regulation and plant defense.