The staminate catkins and pollen were collected from outdoor B. platyphylla trees located at the Northeast Forestry University (45.7662°N, 126.6247°E). Staminate catkin samples were collected on Oct. 20 and Dec. 20, 2021, as well as on Mar. 20 and Apr. 20, 2022, respectively. Pollen samples were collected on May 1, 2022. Each sample was collected from three trees, and three biological replicates were used for each time point.

The collected staminate catkins were fixed in FAA (Formalin-Aceto-Alcohol) with 70% ethanol: acetic acid: 37% formaldehyde = 18: 1: 1, for 48 hours, followed by five rounds of vacuuming for 1.5 hours each. The catkins were then dehydrated using a series of ethanol solutions with increasing concentrations: 50%, 70%, 85%, and 100%, with each concentration step lasting for 3.5 hours. Subsequently, the catkins were dehydrated using a series of xylene solutions with increasing concentrations: 50%, 70%, 85%, and 100%, with each concentration step lasting for 4 hours. The dehydrated material was gradually saturated with paraffin wax. The paraffin mass was sectioned into sections of 7 µm by ultramicrotome. The sections were then stained with 1% toluidine blue for visualization.

The staminate catkins were ground into powder to obtain germline cells. Germline cells were stained with 0.1 μg/mL DAPI (4’,6-diamidino-2-phenylindole) for 30 min and washed three times with PBS (phosphate buffered saline, pH = 7.4). The stained germline cells were observed using a DX51 fluorescence microscope (Olympus, Japan) under bright fields and ultraviolet illumination.

Samples were instantly frozen in liquid nitrogen and stored at -80°C. Total RNA was extracted using RNAprep Pure Plant Kit (DP441, Tiangen, China). The quality of the extracted RNA was assessed using a NanoPhotometer (Implen). The KAPA Stranded RNA-seq Kits (Roche) were employed for sequencing library preparation, and the library quality was evaluated using Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Illumina RNA-Seq was performed by Metware Biotechnology Co. (Wuhan, China) using the Illumina Novaseq6000 system.

The Betula platyphylla v1.1 reference genome was obtained from Phytozome (https://phytozome-next.jgi.doe.gov/info/Bplatyphylla_v1_1). Low-quality reads (Q20 ≤ 50%) were filtered out and Illumina adapters were removed using Cutadapt (Martin, 2011). The quality of the resulting clean reads was assessed using FastQC (Andrews, 2010).

Alignment of the clean reads to the reference genome was performed using Bowtie2 (Langmead and Salzberg, 2012). Gene counts were obtained using featureCounts, and fragments per kilobase of exon per million mapped reads (FPKM) values were calculated using the countToFPKM R package (Liao et al., 2014; Alhendi, 2019). Differentially expressed genes (DEGs) were identified using the DESeq2 R package with the following thresholds: adjusted p-value (p.adjust, FDR method) < 0.05 and absolute log₂ fold change (|log₂FC|) > 1 (Love et al., 2014).

Gene Ontology (GO) enrichment analysis was performed using the clusterProfiler R package (Yu et al., 2012). Gene functional annotations were based on the best match of BLAST searches against TAIR (https://www.arabidopsis.org/) and Swissprot protein databases (https://www.sib.swiss/swiss-prot) of genes. Weighted gene co-expression network analysis (WGCNA) was performed with soft thresholding power β = 8 using the WGCNA R package (Langfelder and Horvath, 2008). The similarity between genes (1-TOM) was calculated, and the genes with similar expression profiles were grouped into the same gene modules using hierarchical clustering and dynamic tree-cut, with a minimum size requirement of 30.

cDNA was synthesized with TransScript® One-Step gDNA Removal and cDNA Synthesis SuperMix (AT311, Transgen). The primers were listed in Supplementary Table 1 . qRT-PCR was performed according to the protocol using TransStart® Top Green qPCR SuperMix (AQ131, Transgen), and three biological replicates were used for each experiment. The relative expressions of all tested candidate genes were normalized to the inner reference genes BpTUB (BPChr11G09309). The relative gene expression level was calculated by the 2^-ΔΔCT method.

The staminate catkins of B. platyphylla, which emerged in June 2021 and dispersed pollen in early May 2022, were subjected to a long period of low temperature during development. To investigate the staminate catkins development across winter, we selected four distinct stages. At S1, staminate catkins were subjected to chilling temperature (0-15°C) in the monthly average temperature from 1°C to 13°C. Staminate catkins at S2-S3 were subjected to freezing temperature (< 0°C) ( Figure 1A ). While at S4, the temperature increased to above zero temperature, and the mean length of staminate catkins increased from 34.2 mm at S1 to 59.6 mm at S4. During freezing temperature, staminate catkins at S2 and S3 exhibited a dehydrated state, which enhanced their freezing tolerance. As the temperature rises, staminate catkins at S4 underwent significant swelling in size and the anthers in catkins approached the pollen-releasing stage ( Figure 1B ; Supplementary Table 2 ).

We investigated the male germline development from B. platyphylla by performing sectioning and DAPI staining. The paraffin section and DAPI staining demonstrated a slow development across the four stages. Each anther of B. platyphylla often contains two locules. At S1, the locules contained interconnected microsporocytes. As the meiosis progresses from S1 to S2, uninuclear microspores differentiated and divided from microsporocytes, which significantly increased the cell numbers within the locules. From S2 to S3, uninuclear microspores with three pores were observed, and the phenotype had no significant changes during three months. At S4, the anther dehiscence occurred, and bicellular microspores were observed ( Figure 1C ).

Our findings suggest a microsporogenesis from S1 to S2, resulting in the formation of uninuclear microspores. From S2 to S3, uninuclear microspore development experiences a freeze-induced slowdown. Finally, at S4, the uninuclear microspores have developed into bicellular pollen, which is ready for release.

RNA-seq was performed on pollen and staminate catkins from S1 to S4. The sequencing yielded a total of 1.2 billion high-quality clean reads. All the clean reads were aligned against the B. platyphylla reference genome (Chen et al., 2021), of which more than 91% were mapped and more than 82% were uniquely mapped, respectively ( Supplementary Table 3 ). To obtain gene expression data, we calculated the fragments per kilobase of exon per kilobase million (FPKM) values of genes from the RNA-seq data for further analysis.

Principal component analysis (PCA) using the gene expression data revealed a strong correlation within three replicates of each experimental group along PC1 (41.99%) and PC2 (19.4%). In contrast to the other experimental groups that could be clearly distinguished from each other, the S2 and S3 experimental groups showed significant overlap, suggesting limited differentiation between these two groups ( Figure 2A ).

We identified 6,513, 365, and 10,377 differentially expressed genes (DEGs) between adjacent stages in the comparisons of S2 vs S1, S3 vs S2, and S4 vs S3, respectively. A total of 12,967 DEGs were found to be differentially expressed in at least one of the examined periods. From S1 to S2, more DEGs (3,314) were down-regulated than up-regulated (3,199). In contrast, more DEGs were up-regulated (5,844) than down-regulated (5,151) from S3 to S4. Only 365 DEGs were found in the comparison of S3 vs S2, indicating a relatively smaller degree of regulatory changes in gene expression at freezing temperature ( Figures 2B, C ). In comparison of S4 and pollen, 4,392 DEGs were upregulated and 6,686 were downregulated ( Supplementary Table 4 ). These findings exhibit the profound impact of substantial temperature fluctuations on male reproduction development, resulting in significant gene expression alteration.

GO enrichment analysis of the DEGs revealed the biological processes associated with different stages of male germline development in B. platyphylla. In the DEGs between S1 and S2, biological processes were enriched in “DNA replication” and “mitotic cell cycle”, suggesting activated cell division and growth during this process. On the other hand, the DEGs between S2 and S3 were found to be enriched in biological processes related to “cellular response to oxygen levels” and “respiratory burst”, suggesting a potential involvement of reactive oxygen species (ROS) in the response to prolonged cold temperature. Between S3 and S4, biological processes related to “cell wall biogenesis” and “water channel activity” were significantly enriched, indicating a focus on cellular structure and water transport in preparation for pollen release ( Figure 2D ). These findings suggested dynamic gene expression patterns and associated biological processes during male germline development in B. platyphylla during overwintering.

To identify gene modules associated with male germline development in B. platyphylla from the four stages to mature pollen, we performed WGCNA to explore relationships between trait factors and gene co-expression modules. The trait matrix represented the developmental characterization in the five experiment groups (S1-S4 and pollen) using binary values. We defined four traits: recovery phase (“00011”) indicates the development process of S4 and pollen; microspore development (“01100”) indicates the development process of S2 and S3; microsporocyte development (“10000”) indicates the development process of S1; cold-inhibited development (“00010”) indicates the development process of S4, respectively. A total of 19,165 genes, filtered based on a mean FPKM value > 0.5, were grouped into 28 co-expression modules ( Figure 3A ; Supplementary Figure 1 ). We calculated pearson correlation coefficients and p-values for each module to assess their association with each trait. The modules significantly correlated with the four traits were identified.

For the recovery phase, the blue module exhibited a significant positive correlation (R = 0.81, p < 0.001, Figure 3B ; Supplementary Figure 2 ). GO enrichment on the genes in the blue module revealed that GO terms of “pollen tube” and “GTPase activity” were enriched ( Figure 3C ). As bicellular microspores emerge, pollen-specific genes such as BpPRK3 (BPChr06G30681), BpPRK4 (BPChr07G02858), BpCPK17 (BPChr04G27193), and BpPPME1 (BPChr03G09691) were upregulated at S4 and pollen ( Figure 4 ). These genes play crucial roles in pollen tube growth (Tian et al., 2006; Myers et al., 2009; Takeuchi and Higashiyama, 2016). CALS5 encodes a callose synthase involved in the synthesis and deposition of callose in the pollen exine wall (Dong et al., 2005). The low expression of BpCALS5 (BPChr13G16035) from S1 to S3, followed by its upregulation at S4, suggested that pollen wall development was achieved by S4 ( Figure 4 ).

On the other hand, the brown module exhibited a significant negative correlation with the recovery phase (R = -0.99, p < 0.001, Figure 3B ; Supplementary Figure 2 ). This module showed significantly enriched in GO terms related to “gene expression”, “mRNA processing”, and “RNA metabolic process” ( Figure 3C ). This implies that cold temperature stimulates the upregulation of gene expression during male reproduction development.

For microspore development, the turquoise module exhibited a significant positive correlation under freezing temperature (R = 0.99, p < 0.001, Figure 3B ; Supplementary Figure 2 ). GO enrichment analysis revealed that the terms associated with “cell division”, “DNA replication”, and “cell cycle G2/M phase transition” were significantly enriched within the turquoise module ( Figure 3C ). During overwintering, cell division emerges as a crucial metabolic process (Perry, 1971). The upregulation of BpEMS1 (BPChr08G07661) at S2 and S3 highlights its essential role in microspore development ( Figure 4 ). The ems1 mutant exhibits male sterility and a lack of tapetal cells, indicating its essential function in microspore development (Zhao et al., 2002). The absence of modules significantly negatively correlated with microspore development suggests that cold temperature primarily regulated male germline development in an activating manner.

For microsporocyte development, the red module showed a positive correlation under chilling temperature (R = 0.97, p < 0.001, Figure 3B ; Supplementary Figure 2 ). These genes belong to GO terms, including “cellular carbohydrate metabolic process” and “glyoxylate metabolic process” ( Figures 3C , 4 ). Among them, BpCPGI (BPChr12G25797), a cytosolic phosphoglucose isomerase was highly expressed at S1, whose homolog in Arabidopsis thaliana is essential for microsporogenesis (Liu et al., 2023) ( Figure 4 ). These findings imply that sugar metabolism plays a crucial role in microsporogenesis under chilling temperature, and ensures sugar supply for subsequent microspore development.

For cold-inhibited development, the green module showed a strong positive correlation (R = 0.99, p < 0.001, Figure 3B ; Supplementary Figure 2 ). The enrichment of GO terms within the green module, including “cell wall” and “water channel activity”, suggested their potential involvement in the late development of staminate catkins ( Figure 3C ). ADR (ANTHER DEHISCENCE REPRESSOR) acts as a key regulator of anther dehiscence by controlling ROS accumulation during secondary thickening in the anther cell wall (Dai et al., 2019). During the key stage of cell wall development, BpADR (BPChr06G00242) was upregulated specifically at S4 ( Figure 4 ). When resource availability was limited, B. platyphylla might prioritize resource allocation towards crucial cellular processes necessary for microspore development and viability by suppressing specific genes and modifying expression patterns in response to low temperature.

Our findings suggest an active sugar metabolism under chilling temperature at S1, which meets the developmental requirements of the microsporocyte and subsequent microspore development. Under freezing temperature from S2 to S3, genes involved in cell division are upregulated, facilitating the development and preparation of uninuclear microspores for division. From S3 to S4, the expression of genes associated with cell wall metabolism and water transport is activated, enabling the successful development of uninuclear microspores into bicellular pollen.

Hub genes, which serve as important regulators within the gene co-expression network, were identified based on thresholds with module membership > 0.8 and gene significance > 0.2. We extracted hub TFs, which may act as central regulators within the network. The blue module encompasses 17 hub TFs potentially involved in male germline development. We found that BpDUO1 (BPChr13G10264) and its potential downstream targets, BpDAZ2 (BPChr09G17992), BpTIP5;1 (BPChr11G07258), and BpDAU2 (BPChr08G17023) exhibited high expression at S4 and pollen, while their expression was low from S1 to S3 ( Figures 4 , 5 ). DUO1, a member of the R2R3 MYB family, is reported to promote the cell division of uninuclear microspores (Borg et al., 2011). The expression pattern of BpDUO1 explains the insufficient cell division of uninuclear microspores under the freezing temperature from S2 to S3.

In the brown module, hub TFs including BpSHI (BPChr05G17937) and BpPDF2 (BPChr14G27036) suggest their essential involvement in the long-term development of staminate catkins from S1 to S3 ( Figure 5 ). These findings align with the well-established roles of SHI (SHORT INTERNODES) and PDF2 (PROTODERMAL FACTOR 2) in regulating stamen development (Ståldal et al., 2012; Kamata et al., 2013). We observed increased expression of BpICE1 (BPChr11G12210) under cold temperature ( Figure 5 ). ICE1 (INDUCER OF CBF EXPRESSION 1) serves as an inducer of CBFs (C-repeat binding factors) and activates various processes related to cold tolerance and adaptation (Chinnusamy et al., 2007). ICE1 is required for maintaining the dehydrated state of anthers, thereby enhancing cold tolerance during pollen development (Wei et al., 2018).

The hub TFs are specifically upregulated at S2 and S3 in the turquoise module ( Figure 5 ). BpAMS (BPChr09G16382), whose homolog in Arabidopsis is known to be involved in anther and microspore development (Ferguson et al., 2017), exhibited specific upregulation under freezing temperature ( Figure 5 ). Moreover, we observed an upregulation of BpMYB65 (BPChr13G00912), BpARF8 (BPChr11G26339), and BpARF17 (BPChr08G06267) during uninuclear microspore development from S2 to S3 ( Figure 5 ). These genes, MYB65 (MYB DOMAIN PROTEIN 65) and ARFs (AUXIN RESPONSE FACTORs) are involved in gibberellic acid (GA)- and auxin-mediated anther development, respectively (Millar and Gubler, 2005; Wang et al., 2017; Ghelli et al., 2023). Furthermore, the increased expression of BpCBF1 (BPChr02G23333) may enhance cold tolerance and the survival of B. platyphylla by regulating the response to cold temperature ( Figure 5 ).

In the red module, BpGAI (BPChr03G02378), is upregulated under chilling temperature at S1 (Ito et al., 2018) ( Figure 5 ). GAI (GIBBERELLIC ACID INSENSITIVE) encoding a DELLA protein, acts as a negative regulator of GA signaling and regulates microsporogenesis and anther dehiscence (Liu et al., 2017; Huang et al., 2022). The upregulation of BpGAI suggests a diminished GA signaling response under chilling temperature at S1. The subsequent downregulation of BpGAI from S2 to S4 suggests that the activation of the GA signaling pathway contributes to the adaption of microspore development to freezing temperature.

In the green module, flowering-associated hub TFs, including BpAP1 (BPChr08G11263), BpMYB24 (BPChr11G26891), and BpWUS (BPChr05G04645), showed a low expression pattern from S1 to S3, and a high expression level at S4 ( Figure 5 ). These TFs have been previously characterized as crucial regulators of stamen development and male identity determination (Deyhle et al., 2007; Song et al., 2011; Huang et al., 2014; Wang et al., 2019). NST1 (NAC SECONDARY WALL THICKENING PROMOTING FACTOR1) and MYB108 are involved in anther dehiscence and pollen release (Mitsuda et al., 2005; Xu et al., 2019). BpNST1 (BPChr09G20586) and BpMYB108 (BPChr08G15721) showed significant upregulation at S4 ( Figure 5 ), implying an activated cell wall metabolism occurred at this stage during anther development.

Our findings reveal the dynamic temperature-dependent regulation of microspore and anther development by GA and auxin responses. At chilling temperature, there is a low GA response at S1, which transitions to a high auxin and GA response from S2 to S3 under freezing temperature. Among the hub TFs, BpAMS, BpMYB65, BpARF8, and BpARF17 may regulate microspore and anther development under freezing temperature from S2 to S3. The hub TFs involved in anther dehiscence and uninuclear microspore division show a low expression from S1 to S3 and a high expression at S4.

Vernalization pathway play a crucial role in the inflorescence bud dormancy in addition to triggering the transition from vegetative to reproductive growth (Lu et al., 2022). We examined the expression of genes related to vernalization pathway. The expression level of BpFLC (BPChr08G03183), was induced by freezing temperature from S1 to S2, then gradually decreased from S2 to S4, along with BpFRI (BPChr13G03877). FRI (FRIGIDA) activates the transcription of FLC, and this activation was reported to be progressively repressed by epigenetic modification at the FLC locus (Li et al., 2018; Zhu et al., 2021). VIP4 (VERNALIZATION INDEPENDENCE 4) was reported to positively regulate FLC expression, and its mutation renders the absence of FLC expression (Zhang and Van Nocker, 2002). BpVIP4 (BPChr04G09387) exhibited high expression from S2 to S3. Specifically, BpFT (BPChr05G17496) showed high expression at S4, low expression at S1, and was undetectable at S2 and S3, indicating its repression under freezing temperature, possibly due to the high expression of its repressor, BpFLC ( Figure 6 ). These observations suggest the involvement of vernalization pathway in the overwintering process of staminate catkins in B. platyphylla.

We used qRT-PCR to verify the reliability of the RNA-seq data. A total of eight genes that may be involved in the development of staminate catkins in B. platyphylla during overwintering were selected. These genes included BpADR, BpAMS, BpAP1, BpDUO1, BpFT, BpINV4, BpMYB65 and BpNST1. We found that the expression patterns of these genes obtained from qRT-PCR were consistent with the RNA-seq data ( Figure 7 ; Supplementary Figure 3 ). The results revealed the reliability of the RNA-seq data.