1. Introduction

Soybean is a cornerstone of global agriculture, serving as a principal source of plant-derived oil and protein, and playing a key role in nitrogen fixation in agroecosystems (Du et al., 2023; Zhang et al., 2022). Given that the global population is projected to reach 9.7 billion by 2050, the growing demand for protein-rich food and feed resources underscores the urgency to achieve sustainable soybean production (Huang et al., 2024). Developing diverse and genetically rich gene pools is critical to expanding allelic diversity for breeding, enabling the development of high-yielding, climate-resilient soybean genotypes that address future food security challenges. Such initiatives enhance soybean adaptability to changing environments and provide a robust genetic foundation for meeting complex global protein supply demands in the future.

Reliable genomes serve as valuable representatives of entire species, enabling inferences of key conservation genomic parameters (e.g., heterozygosity levels, inbreeding coefficients, and demographic history). The genomic era of soybean began with the landmark assembly of the cultivar Williams 82 (Schmutz et al., 2010), followed by high-quality genomes of other elite cultivars including Zhonghuang 13 (Shen et al., 2018), W05 (Xie et al., 2019), Nongdadou2 (Zhang et al., 2024a), and Tianlong 1 (Sheng et al., 2025). Recent breakthroughs have yielded several T2T assemblies of soybean genomes, such as Jack (Huang et al., 2024), Zhonghuang 13 (Zhang et al., 2023), Williams 82 (Wang et al., 2023), Yesheng71, Yundou1 (Jia et al., 2024), and YSD56 (Lian et al., 2025). These genomic resources have been invaluable for functional genomics and molecular breeding studies, facilitating gene discovery and trait analysis. However, the genetic basis of many complex agronomic traits, particularly those present in unique germplasms like photoperiod-sensitive genic male sterility (PGMS) lines, remains underexplored due to the lack of high-quality reference genomes for such materials. This gap mirrors early-stage challenges observed in other crops, such as watermelon, where multi-variety analyses proved essential for understanding crop evolution (Zhang et al., 2024b). Male sterility systems are pivotal for commercial hybrid seed production in crops, as they enable the efficient utilization of heterosis while eliminating the labor-intensive process of manual emasculation (Vasupalli et al., 2025). In soybean, the development of PGMS lines offers a flexible 'two-line' breeding system, which overcomes the limitations of traditional 'three-line' systems and is crucial for enhancing yield potential (Khan et al., 2023; Vasupalli et al., 2025).

Here, we report a T2T reference genome of landrace 88-428BY, a germplasm released in 1988 in Taiyuan, Shanxi Province, China (Fig. 1a). This line represents China's first reported PGMS resource (Wei, 1991). It was discovered as a natural mutant in 1988 during field trials screening for apomictic resources, identified within a local landrace known as 'Local Tumeidou'. Unlike lines derived from radiation or chemical mutagenesis, 88-428BY arose from spontaneous natural mutation. Although its precise genetic pedigree remains unclear, its core breeding value is well-documented: it enables a highly flexible 'two-line' breeding system. Under long-day conditions (14.5–15.0 h; LD), the plant exhibits normal pollen fertility and produces viable pods; conversely, under short-day conditions (11.5–13.0 h; SD), pollen development is arrested, leading to complete sterility and the formation of fleshy, non-seed-bearing pods instead of normal fertile ones (Yang et al., 2025; Fig. S1). This 'dual-purpose' feature significantly expands the freedom of parental selection for hybrid combinations. Decades of continuous photoperiod and environmental factor trials have confirmed that the fertility of 88-428BY is strictly regulated by photoperiodic changes, showing no response to temperature fluctuations within the tested ranges.

To construct a comprehensive reference genome, we employed a hybrid assembly strategy using PacBio high-fidelity (HiFi), and Oxford Nanopore Technologies (ONT) ultra-long reads to generate backbone contigs. These contigs were subsequently scaffolded into a chromosome-level assembly with the assistance of high-through chromosome conformation capture (Hi-C) datasets. Gap filling resolved remaining sequence gaps, yielding a complete T2T assembly for 88-428BY. This genomic resource provides a valuable foundation for advancing soybean functional genomics, evolutionary studies, and molecular breeding, particularly for leveraging its unique photoperiod-sensitive sterility traits.

2. Material and methods 2.1. Sample preparation

The 88-428BY plants used in this study were collected from the Taiyuan City (37.78°N, 112.57°E), Shanxi Province, China. Upon collection, all samples were immediately frozen in liquid nitrogen and stored at −80 ℃ until subsequent use. High-quality genomic DNA was extracted from fresh leaf tissues using the DNeasy Plant Mini Kit (Qiagen, Hilden, Germany). After extraction, the quality and integrity of the DNA were assessed via 1% agarose gel electrophoresis and a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), with the A260/A280 ratio maintained at 1.8–2.0 and no significant degradation observed.

2.2. DNA sequencing

For Nanopore sequencing, a DNA library was constructed using the SQK-LSK110 ligation kit (Oxford Nanopore Technologies, Oxford, UK), and the purified library was sequenced on a PromethION instrument (Oxford Nanopore Technologies, Oxford, UK). Base calling was performed using Guppy with default parameters, and reads were filtered for mean_qscore_template ≥ 7. NanoPlot v.1.41.0 (De Coster et al., 2018) was then used to filter the Nanopore reads.

For PacBio HiFi sequencing, a long-read library (~15 kb) was constructed from genomic DNA using SMRTbell® Express Template Prep Kit 2.0 (PacBio, USA) and sequenced on a PacBio Revio instrument (PacBio, USA).

For Hi-C sequencing, a DNA library was prepared following the manufacturer's instructions. Then the library was sequenced by BGI Co., Ltd (Shenzhen, China) on a BGISEQ-500 instrument to generate paired-end reads (150 bp). Low-quality sequences were removed from raw data using fastp v.0.24.0 (Chen et al., 2018) to obtain clean data.

2.3. Genome assembly and evaluation

Contigs were assembled using Hifiasm v.0.19.5 (Cheng et al., 2021) from a combination of PacBio HiFi and ONT ultra-long reads. The Hi-C data were then used for chromosome anchoring and assembly into chromosomes using Juicer v.1.6 (Durand et al., 2016) and 3d-dna pipeline v.201008 (Dudchenko et al., 2017), followed by manual correction with Juicebox v.1.11.08 (Durand et al., 2016). To close gaps, PacBio HiFi reads and ONT long reads were used with LR_Gapcloser (Xu et al., 2018) with parameters of '-t 40 -s nanopore -m 1,000,000 -v 10000'. For further refinement of the genome, the T2T assembly was polished using Racon v.1.5.0 (Vaser et al., 2017).

Benchmarking Universal Single-Copy Orthologs v.5.8.0 (BUSCO) analysis was performed by searching against the conserved 1614 Embryophytagene set (Waterhouse et al., 2018). We applied Merqury v.1.3 (Rhie et al., 2020) with a K-mer value of 17-bp to estimate the quality value (QV). Additionally, ONT and PacBio HiFi reads were aligned with Minimap2 v.2.28 (Li, 2018). The mapping results wear visualized using the "karyoplote" R package (Gel and Serra, 2017).

2.4. Genome annotation

Tandem repetitive sequences were identified using Tandem Repeats Finder v.4.10.0 (Benson, 1999). The interspersed repeats were determined using a combination of homology-based and de novo approaches. The RepBase was used to identify transposable elements (TEs) by searching against the 88-428BY genome assembly at the DNA and protein levels using RepeatMasker v.4.0.7 and RepeatProteinMask, respectively. A de novo repeat library was customized using RepeatModeler v.1.0.4 and LTR_FINDER v.1.0.7 (Xu and Wang, 2007) and then imported to RepeatMasker v.4.0.7 to identify repetitive elements.

Protein sequences from six plant genomes, including two soybean assemblies (ZH13-T2T and Wm82-T2T), Arabidopsis thaliana (Hou et al., 2022), Phaseolus vulgaris, Vigna angularis, and Medicago truncatula, were downloaded for gene prediction. Homology-based prediction was conducted using GeMoMa v.1.9 (Keilwagen et al., 2019), while AUGUSTUS v.3.2.3 (Stanke and Morgenstern, 2005) was applied to identify coding regions within the repeat-masked genome. The RNA-seq data utilized for gene annotation were generated from five distinct tissues (flower, leaf, pod, root, and stem) collected from 88 to 428BY plants grown under normal conditions. Following quality control, the resulting RNA-seq clean reads were mapped to the assembly using Hisat2 v.2.0.1 (Kim et al., 2015). Transcripts were then assembled and candidate coding regions identified using Stringtie v.1.2.2 (Kovaka et al., 2019) and TransDecoder v.5.7.1. Additionally, RNA-seq reads were de novo assembled with Trinity v.2.15.1 using (Grabherr et al., 2011) parameters of '–max_memory 200 G–CPU 40–min_contig_length 200–genome_guided_bam merged_sorted.bam–full_cleanup–min_kmer_cov 4–min_glue 4–bfly_opts '-V 5– edge-thr = 0.1–stderr'–genome_guided_max_intron 10000'. Assembled transcripts were aligned to the genome using the Program to Assemble Spliced Alignments v.2.4.1 (PASA; Haas et al., 2008), yielding the gene structures from valid alignments. Gene models from all three approaches were integrated into a non-redundant dataset using EVidenceModeler v.1.1.1. The resulting gene models were finally refined using PASA to obtain untranslated regions and alternative splicing variation information. Functional annotation of protein-coding genes followed the protocol described in Sorghum T2T genomic studies (Li et al., 2024). The Plant Transcription Factor Database (https://planttfdb.gao-lab.org/) was used to identify transcription factors (TFs; Jin et al., 2017).

The telomeres were identified using QuarTeT v.1.2.2 (Lin et al., 2023) method with the "-c plant" option. BLAST v.2.16.0 (Altschul et al., 1990) was applied to align the Cent-Gm1 and Cent-Gm2 to 88-428BY genome with an E-value of 1e-6, and then BEDtools v.2.31.1 (Quinlan and Hall, 2010) was used to merge the results with the parameter "-d 1,000,000". If a region with lengths longer than 500 kb, it was classified as a centromere region.

The transfer- (tRNAs) genes were identified using tRNAscan-SE v.2.0.11 (Lowe and Eddy, 1997). The ribosomal (rRNAs) fragments were detected through alignment to Arabidopsis and rice template rRNA sequences via BLAST v.2.16.0 (Altschul et al., 1990) with an E-value of 1e-5. The micro- (miRNAs) and nuclear- (snRNAs) genes were determined by searching against the Rfam database using INFERNAL v.1.1.1 (Nawrocki and Eddy, 2013).

2.5. Ortholog identification

To establish high-confidence orthologous relationships between gene models in the 88-428BY and classical Wm82 genomes, we performed a bidirectional DIAMOND v.0.9.30 (Buchfink et al., 2015) search using protein sequences, with orthologs identified as reciprocal best hits meeting stringent criteria including syntenic alignment, consistent strand orientation and conservation of gene architecture. For cases with multiple homologous hits, only the top-scoring gene pair was retained to ensure unambiguous one-to-one orthology.

2.6. Identification of structure variations (SVs) and SV hotspot regions

Five T2T level soybean genomes (Wm82, Jack, ZH13, Yundou1, and Yesheng71) were first aligned to the 88-428BY genome, using MUMmer v.4.0.0rc1 (Marçais et al., 2018) with parameters of '-g 1000 -c 90 -l 40'. After filtering low-confidence alignments with delta-filter program, the SV calling was performed using the SyRI software v.1.7.0 (Goel et al., 2019). The SURVIVOR v.1.0.7 was utilized to merge SVs, using command line "SURVIVOR merge 1000 1 1 1 0 50" (Jeffares et al., 2017), with only SVs ranging from 50 bp to 100 kb in length retained. Based on the genic regions overlapping with SVs, we annotated the identified SVs with ANNOVAR v.2017-06-01.

SV hotspots were defined as genomic regions where the density of structural variations significantly exceeded the genome-wide average level. To systematically define SV hotspot regions, we first calculated the density of SV breakpoints across the genome using a sliding window approach. Specifically, we calculated the distribution of SV breakpoints for each 200 kb window (with a 100 kb step size) along each chromosome. Then, all 200 kb windows were ranked in descending order according to the numbers of SVs within the window. We defined the top 10% of all windows with the highest frequency of SV breakpoints as SV hotspots, then merged all of the continuous hotspot windows as the "hotspot regions".

2.7. Transcriptome analysis

RNA-seq raw reads from anthers at developmental stages 4 (S4), 8 (S8), and 12 (S12) under long-day (LD: 14.5–15.0 h light) or short-day (SD: 12 h light/12 h dark) conditions, which were generated in our previous study (Yang et al., 2025), were retrieved. After quality control using fastp v.0.24.0 (Chen et al., 2018), clean reads were aligned to the 88-428BY reference genome with HISAT2 v.2.0.1 (Kim et al., 2015). Raw read counts per gene per library were generated from the sorted SAM files using htseq-count v.2.0.5 (Anders et al., 2015). Transcript expression levels were quantified as fragments per kilobase of exon per million mapped reads (FPKM; Roberts et al., 2011) method. Additionally, raw read counts were used to analyze differential gene expression via the R package DESeq2 v.1.18.1 (Love et al., 2014). A gene was defined as a differentially expressed gene (DEG) if it met the criteria of adjusted p-value ≤ 0.05 and absolute fold change ≥ 2. GO enrichment analysis of DEGs was performed using ClusterProfiler v.4.8.1 (Wu et al., 2021), while KEGG enrichment analysis was conducted via KOBAS v.2.0 (Xie et al., 2011).

Gene co-expression networks were constructed using the weighted gene co-expression network analysis (WGCNA) package v.1.72 in R software v.4.1.2 (Langfelder and Horvath, 2008). The optimal soft-thresholding power was set to 14, and the minimum gene module size was designated as 30. The hub gene was determined using cytoHubba (Chin et al., 2014) with the degree algorithm.

2.8. Validation of gene expression via quantitative real-time PCR (qRT-PCR)

Total RNA was isolated, after which reverse transcription was performed with the PrimeScript^TM RT Reagent Kit (TaKaRa, Tokyo, Japan). The synthesized first-strand cDNA was then subjencitedcted to amplification reactions in strict adherence to the standard protocols supplied with the SYBR Green-based RealMasterMix kit (SYBR Green, TIANGEN). Relative expression levels of 10 selected DEGs (ANS2, BBI, DFR2, Gm428.02G00900, Gm428.02G08770, Gm428.15G13580, Gm428.20G09320, HIPP09.1, INV1, NAC25) across anthers tissue at different developmental stages and experimental treatments were quantified via the 2^−ΔΔCt algorithm (Livak and Schmittgen, 2001). Actin11 (accession number: LOC100781831) was designated as the reference gene to adjust for discrepancies in initial template concentrations (Sharmin et al., 2020). All assays were carried out with three biologically independent replicates. The primers were designed using Primer 3 software (Table S20).

3. Results 3.1. A high quality T2T genome assembly and annotation for soybean 88-428BY

We generated a T2T reference genome for soybean landrace 88-428BY using combined sequencing data: 72.99 Gb PacBio HiFi reads (~72 × depth, N50 = 18.58 kb), 254.51 Gb ONT ultra-long reads (~252 × depth, N50 = 53.79 kb), and 244.43 Gb Hi-C data (~242 × depth; Table S1). The initial contig assembly (1.03 Gb, N50 = 43.58 Mb) was refined with Hi-C data, gap-closed and polished to a 1.01 Gb genome (N50 = 52.05 Mb) covering 20 chromosomes (Table S2 and Table S3). This assembly is comparable to published T2T soybean genomes (Table 1, Fig. 1b and Table S4). Validation confirmed exceptional quality, with 99.97% of HiFi reads and 100% of ONT reads mapping to the assembly (Table 1). The assembly shows uniform read coverage (Fig. S2), an overall quality value of 57.71 (ranging from 51.51 to 78.44; Table S5), strong Hi–C interaction continuity (Fig. 1c), and a 99.88% BUSCO completeness score (Table 1 and Fig. 1d). Genome annotation identified 637.94 Mb of repetitive sequences, accounting for 63.03% of the genome. This content was higher than that in other reported T2T soybean genomes (Tables 1 and S6), with LTR retrotransposons (51.56%) as the dominant interspersed repeat type (Table S7). We predicted 55,292 protein-coding genes (Table S8), 77.58% of which showed one-to-one orthology with the Wm82 genome (Table S9), and the gene set achieved 99.19% BUSCO completeness for these genes (Table S10). In total, 54,108 (97.86%) protein-coding genes were annotated with functional information (Table S11). Further, the length distribution of genes, CDSs, exons, and introns among the recently published common bean T2T genome (Wang et al., 2025) and four different varieties of soybeans (88-428BY, Wm82, Yesheng71, and ZH13), showed that the exon and intron length distributions were consistent across species, with small fluctuations in the mRNA and CDS length distributions (Fig. S3). Additionally, 3662 TFs were identified, with bHLH, MYB and ERF families being the most abundant (Fig. S4). We annotated 4156 non-coding RNAs in the 88-428BY assemblies, including 230 miRNAs, 1182 tRNAs, 816 rRNAs, and 1928 snRNAs, respectively (Table S12).

We identified 36 telomeric repeats (motif: AAACCCT) across 20 chromosomes, resulting in 80% of chromosomes harboring telomeres at both termini (Fig. 2 and Table S13). This proportion meets and exceeds the criteria for T2T assemblies (Xie et al., 2024). Centromeres ranged from 0.96 to 4.16 Mb in length and were dominated by tandem repeats and LTRs (Tables S14 and S15), consistent with findings in the Jack T2T genome (Huang et al., 2024). A total of 34 genes were annotated in centromeric regions, showing significant enrichment in terms related to DNA integration, male meiotic nuclear division, female meiotic nuclear division, and galactosylgalactosylxylosylprotein 3-beta-glucuronosyltransferase activity (Fig. S5).

3.2. Transcriptomic responses of 88-428BY to photoperiod treatment and co-expression network analysis

To investigate photoperiod sensitivity, we analyzed anther transcriptomes under LD and SD conditions. RNA-seq clean reads showed higher alignment rates to the 88-428BY genome (87.03–90.83%) than to the Wm82 reference (Table S16; ~83%; Yang et al., 2025). Among the annotated genes, 74.46% (41,173) were expressed (FPKM ≥ 1), including 64.71% (22/34) of the centromeric genes (Fig. S6). We identified 691 common DEGs across anther development stages (Fig. S7 and Table S17). These genes were enriched in KEGG pathways such as pentose and glucuronate interconversions, plant–pathogen interaction, galactose metabolism, starch and sucrose metabolism, and flavonoid biosynthesis, as well as in GO terms including enzyme inhibitor activity and carbohydrate metabolic process (Fig. 3a and b). Results from the qRT-PCR assays demonstrated that the change trends of relative expression levels and FPKM values for the 10 DEGs were generally consistent (Fig. S8).

WGCNA identified 16 co-expression modules (Fig. 3c). The darkmagenta module correlated positively with LD conditions and negatively with SD conditions, comprising 143 genes including 20 TFs. Its top 20 core genes included 6 common DEGs, such as anthocyanidin synthase (ANS2) and dihydroflavonol-4-reductase (DFR2; Fig. 3d). In contrast, the brown4 module showed the opposite correlation with photoperiods, containing 101 genes including 14 TFs. One of its top 20 core genes was a common DEG, identified as a bowman-birk type protein inhibitor (BBI; Fig. 3e). The qRT-PCR validation assays verified the expression profiles of ANS2, DFR2, and BBI (Fig. S9).

3.3. Discovery of SVs and SV hotspot regions

We cataloged 35,236 non-redundant SVs, comprising 17,954 insertions and 17,282 deletions. Most SVs resided in intergenic regions (72.02%), while 2.10% were in exons (Table S18). Among the exonic SVs, 599 were predicted to disrupt protein structure, affecting 660 genes that were significantly enriched in terms like DNA integration and zinc ion binding (Fig. S10). Notably, one common DEG, ammonium transporter 1 (AMT1), harbored a frameshift deletion and exhibited higher expression under LD conditions.

SVs were unevenly distributed across the genome (Fig. 4a). We defined "SV hotspot regions" as the top 10% of 200-kb windows ranked by SV breakpoint density, identifying 522 hotspots harboring 11,730 genes (Table S19). Genes in these hotspots exhibited significantly lower expression (Fig. 4b). Intriguingly, 35 genes from the LD-associated darkmagenta module and 24 genes from the SD-associated brown4 module resided in SV hotspots (Fig. 4c–d). These included two key TFs that were also common DEGs: NAC transcription factor 25 (NAC25; darkmagenta module) showed increasing expression under LD, while AP2-EREBP transcription factor (ERF113; brown4 module) increased under SD during anther development. These integrated analyses link specific SVs and genomic regions to transcriptional networks underlying photoperiod response.

4. Discussion

High-quality soybean genome assemblies are poised to revolutionize crop improvement through integrated breeding and functional genomics. Although several T2T reference genomes have recently been published, genomic resources for PGMS germplasms remain scarce, limiting the dissection of photoperiod-sensitive male sterility mechanisms. This study presents the first high-quality T2T genome assembly of Glycine max landrace 88-428BY, a Chinese PGMS variety. While current soybean research often focuses on quantitative traits such as seed size (Liang et al., 2024), oil/protein content (Zhang et al., 2025), and plant height increase (Sun et al., 2025), the 88-428BY germplasm provides a unique model for molecular dissection of photoperiod-sensitive male sterility. This trait is crucial for hybrid seed production, as conditional male sterility enables efficient hybrid seed production (Khan et al., 2023; Vasupalli et al., 2025).

The assembly demonstrates outstanding completeness and continuity. It exhibits an N50 of 52.05 Mb and a BUSCO completeness score of 99.88%. Notably, telomeres were assembled at both ends of 80% of the chromosomes, exceeding common T2T standards. Repetitive sequences constitute 63.03% of the genome, a proportion higher than that reported for other T2T soybean accessions (35.67–57.57%). This elevated repetitive content may reflect the distinctive genomic architecture of 88-428BY and underscores the significant contribution of repetitive elements to soybean genome plasticity.

Our analysis identified 36 telomeres carrying the conserved AAACCCT repeat motif, a feature consistent across five T2T soybean genomes. In addition, 20 centromeric regions were annotated, with an average length of 2.27 Mb, primarily composed of tandem repeats and LTR retrotransposons. Comparative analysis revealed that the centromeric repeat composition of accession 88-428BY is largely conserved compared to the Jack T2T genome (Huang et al., 2024), supporting the hypothesis that centromeric repetitive elements contribute consistently to chromosomal stability during meiosis. Furthermore, GO enrichment analysis of genes located in centromeric regions showed significant enrichment for terms related to both male and female meiotic nuclear division. The enrichment of meiotic division-related genes in centromeric regions, coupled with the known role of centromeric transcripts in pollen development (e.g., OsMRPL15 in rice), suggests that centromere-associated genes may contribute to fertility regulation in soybean (Xie et al., 2023).

The transcriptomic analyses presented herein expand on our prior findings (Yang et al., 2025), revealing both consistency and substantial advancements. In both studies, KEGG enrichment analysis of DEGs identified overlapping pathways critical to the PGMS trait, including flavonoid biosynthesis, plant–pathogen interaction, galactose metabolism, starch and sucrose metabolism, and pentose and glucuronate interconversions. This consistency validates the robustness of the core transcriptional response governing photoperiod sensitivity. Notably, the adoption of the high-quality 88-428BY T2T genome as a self-reference genome conferred distinct advantages. A key improvement was the elevation of the RNA-seq read unique mapping rate from ~83% (using the Wm82 reference in our prior work) to a range of 87.03–90.83%. This enhanced read alignment accuracy facilitated more precise gene quantification and co-expression network construction. Whereas our previous analysis identified 17 modules and underscored the broad regulatory functions of MADS-box and MYB TFs, the current T2T genome-based WGCNA yielded 16 modules with a more refined architecture. By focusing on the intersection of hub genes and shared DEGs, we pinpointed specific functional effectors: ANS2 and DFR2 in the LD-responsive darkmagenta module, and BBI in the SD-responsive brown4 module. Given the established role of flavonoid biosynthesis in pollen wall formation, the tightly coordinated expression of ANS2 and DFR2 likely contributes to maintaining pollen fertility under long-day conditions (Fang et al., 2021). These findings extend our earlier transcriptomic investigation by contextualizing the results within a T2T genome framework, enabling more accurate genomic mapping.

We identified a total of 35,236 non-redundant SVs, of which 2.10% are localized within exonic regions. Among these exonic SVs, 599 were predicted to disrupt protein structure. A notable example is a frameshift deletion identified in AMT1, a commonly differentially expressed gene. AMT1, encoding an ammonium transporter, exhibits higher expression under LD conditions. This suggests that its frameshift deletion may impair nitrogen transport to developing pollen, a process critical for pollen maturation given that nitrogen deficiency can induce male sterility in crop (Khan et al., 2023). Specifically, nitrogen deficiency caused by impaired AMT1 function can result in reduced pollen viability and seed set, contributing to male sterility, as observed under high-temperature stress where nitrogen metabolism is disrupted (Khan et al., 2023).

Compared with single nucleotide polymorphisms (SNPs), SVs generally exert more substantial effects on gene expression (Alonge et al., 2020). We detected 522 SV hotspot regions, which harbored 11,730 genes exhibiting reduced expression. These regions include significant differentially expressed TFs, such as NAC25 and ERF113. In maize, NAC transcription factors are known to regulate leaf senescence and male fertility (Yuan et al., 2023). Additionally, the NAC transcription factors negatively regulate female fertility by promoting ovule senescence in Arabidopsis thaliana (Van Durme et al., 2023). AP2/ERF transcription factors indirectly influence plant fertility by modulating hormone signaling pathways, particularly through brassinosteroid-mediated development, though their primary role is in abiotic stress responses (Ma et al., 2024). Therefore, the differential expression of NAC25 and ERF113 in SV hotspots suggests their potential involvement in hormone-mediated regulation of pollen development under different photoperiods. While SVs offer valuable insights, the genetic basis of male sterility in soybean can also be attributed to SNPs and small insertions/deletions (InDels), which were not the primary focus of this SV-based study. Recent evidence shows that the ms5 male sterility phenotype is caused by a specific 15-bp deletion and a nucleotide substitution within the GmMSH4 gene, leading to aberrant splicing (Nagayama et al., 2025). Similarly, the identification of the mst-M male-sterile gene via SNP-based mapping further underscores the critical role of small-scale variants in determining fertility (Zhao et al., 2019). Therefore, future studies leveraging large-scale population whole-genome sequencing to generate high-density SNP/InDel datasets, integrated with the T2T reference genome constructed in this study, will be critical for fully delineating the PGMS mechanism in 88-428BY.

5. Conclusions

In summary, we present a high-quality, telomere-to-telomere genome assembly for the PGMS soybean landrace 88-428BY. This resource provides a foundational reference for studying photoperiod-sensitive male sterility in soybean. Integrated multi-omics analyses suggest a model where male fertility under long-day conditions may be maintained by the coordinated expression of flavonoid biosynthesis genes and the potential functional impact of an SV in the ammonium transporter AMT1. Conversely, male sterility under short-day conditions may involve the upregulation of a protease inhibitor (BBI) and structural variation-associated dysregulation of key TFs (NAC25 and ERF113). This study not only fills a genomic resource gap for unique soybean germplasms but also identifies candidate genomic regions and genes for further functional validation.

Acknowledgments

This work was supported by Major Science and Technology Special Program of Shanxi Province (No.202201140601025), the Natural Science Foundation of Shanxi Province (No.202303021211085), the China Postdoctoral Science Foundation and Shanxi Agricultural University Outstanding Doctoral Startup Project (2023BQ95). We thank Baoguo Wei for his kind provision of the 88-428BY.

CRediT authorship contribution statement

Yu-Hua Yang: Writing-original draft, Formal analysis, Methodology, Software, Data curation, Funding acquisition. Zhi-Yuan Bai: Data curation, Formal analysis. Yan Chen: Formal analysis. Meng-En Nie: Formal analysis. Hai-Gang Wang: Formal analysis, Writing-review & editing. Hai-Ping Zhang: Formal analysis, Writing-review & editing. Rui-Jun Zhang: Formal analysis, Writing-review & editing. Zhi-Xin Mu: Supervision, Resources, Conceptualization.

Data accessibility statement

The data that support the findings of this study are openly available in the China National GeneBank (CNGB) Nucleotide Sequence Archive (CNSA) under accession number CNP0007030.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2026.04.006.