A total of 193 leaf samples of Phyllostachys edulis were collected from 37 localities in natural distribution of China during several filed campaigns ( Supplementary Table S1 , Supplementary Figure S2 ). In general, more than five individuals per population were sampled, while the following exceptions: populations No. 7, 13, 27, and 29 (three populations each), No. 14 and 19 (four populations each), and No. 31 and 32 with only two individuals, respectively. The minimally sample size for populations No. 31 and 32 from Sichuan Yibin reflects the challenging field conditions in this region, where low natural regeneration rates of moso bamboo, wildfire and steep terrain significantly limited natural populations sample accessibility. The sampled young leaves were directly silica dried and voucher samples were deposited at Nanjing Forestry university.

Total genomic DNA of 193 moso bamboo samples from 37 populations was extracted from dried leaves using the improved CATB method (Doyel et al., 1996). Subsequently, the DNA purity and concentration were checked by the NanoDrop 1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and assessed using agarose gel electrophoresis with a concentration of 1.0%. The GBS library was constructed following the protocol developed by Qi et al. (2018), with a selected size range of 300–600 bp. Briefly, cleaned DNA was digested using 8 U MspI and 8U PstI-HF (NEB), followed by T4-ligase ligation, random shearing and size selection. The recovered DNA fragments were used as templates for PCR amplification. After purification and quality checking, the resulting library was sequenced with paired-end (2×150 bp) reads on the Illumina Novaseq PE150 platform (Illumina, San Diego, USA) at Beijing Kaitai Gene Technology Co., Ltd. (Beijing, China).

All of the raw sequencing data have been submitted to the NCBI Short Reads Archive database with the accession number PRJNA1133081, ensuring that the data are publicly accessible and reproducible.

Raw reads were filtered using fastp (version 0.23.1) to obtain clean reads (Chen et al., 2018). Data quality filtering was performed in three sequential steps: first, reads containing adapters were removed; second, reads containing more than 10% of ambiguous bases (N) were discarded; and third, low-quality reads were eliminated, defined as those in which the number of bases with a quality score (Q) ≤ 10 exceeded 50% of the total read length. Clean reads were then aligned to the previously assembled genome of P. edulis using BWA (version 0.7.8) (Li and Durbin, 2010). Subsequently, based on the comparison results, the Haplotypecaller module of GATK software (version 4.1.8) was employed for SNP and InDel detection (Mckenna et al., 2010). The raw SNP and InDel dataset were further filtered by vcftools software (concrete parameters: ‘–maf 0.01 –max-missing 0.8 –minDP 4 –min-alleles 2 –max-alleles 2’, version 0.1.16) (Danecek et al., 2011). The filtering criteria were applied as follows: reads were retained for sites with a minimum depth of 4; loci with a minor allele frequency (MAF) less than 0.01 were removed; and loci with SNP genotyping deletion rate exceeding 20% were excluded.

Population genetic analysis was conducted based on the filtered high-quality SNP data. Phylogenetic tree among individuals was constructed using the maximum likelihood (ML) algorithm of Evolview (http://www.evolgenius.info/evolview). The population structure was analyzed using Admixture (version 1.3.0) (Porras-Hurtado et al., 2013), with the assumed number of subgroups (K) ranging from 2 to 10. The optimal K value determined based on the minimum cross-validation (CV) error. Principal component analysis (PCA) and kinship analysis among individuals were performed using gcta (version 1.26.0) (Yang et al., 2011).

Genetic diversity indicators, including Hardy Weinberg equilibrium (HW-P), Number of alleles (Na), effective numbers of alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He), polymorphism information content (PIC), nucleotide diversity (π), fixation index (Fst) and inbreeding coefficient (Fis, Fit) were calculated by Genepop (version 4.7) in R (version 4.2.3) (Rousset, 2008). The Reynolds genetic distance between populations (DR) was calculated by Fst, D R = − ln ( 1 − F s t ) .

Mantel tests were performed to evaluate the associations between genetic distance, geographic distance and environmental distance with 9999 permutations in the R package vegan. AMOVA analysis was carried out by Arlequin (version 3.5.2.2) to quantify the genetic variation among and within the genetic groups of moso bamboo inferred by the population structure analyses. DnaSP (version 6.0) was used to calculate the gene flow (Nm) of moso bamboo population (Rozas et al., 2017).

Functional annotation of all candidate genes was performed using EGGNOG-Mapper (http://eggnog-mapper.embl.de/). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed on the annotation results using TBtools (version 1.9) (Chen et al., 2020).

To validate SNPs identified from GO and KEGG enrichment, three of nonsynonymous mutation sites and one of synonymous site in the coding region of functional genes related to temperature stimulation response were selected for PCR amplification and Sanger sequencing. Flanking sequences of selected SNPs were extracted from the reference genome and PCR primers were designed with Premier Primer 5 (Singh et al., 1998) and synthesized by Beijing Tsingke Biotech Co., Ltd (Beijing, China) ( Supplementary Table S2 ). DNA samples of ten individuals from each of the three subpopulations were randomly selected as templates for PCR amplification. PCR reaction system: 1.1×T3 Super PCR Mix 22μL, 10μM Primer F 1μL, 10μM Primer R 1μL, Template (gDNA) 1μL. PCR reaction conditions: 98°C 2 min; 98°C 10 s, 60°C 30 s, 72°C 45 s, 35 cycles; stretched at 72°C for 5 min and stored at 4°C. The amplified PCR products were sequenced by the Sanger method and sequencing results were analyzed with BioEdit (Hall, 1999).

The Tajima’s D and Fu ‘s Fs values of moso bamboo were calculated by Arlequin (version 3.5) (Fu and Li, 1993). The site frequency spectrum (SFS) of each subgroup was calculated by the easySFS.py script (https://github.com/isaacovercast/easySFS). Stairway Plot (version 2.1) was used to infer the demographic history of moso bamboo (Liu and Fu, 2020).

Species data was obtained from the Global Biodiversity Information Facility (GBIF, https://www.gbif.org) and 19 bioclimatic factors ( Supplementary Table S3 ) with a spatial resolution of 2.5 arc-minutes (~ 5 km) were downloaded from the PaleoClim (http://www.paleoclim.org/) database. Using ArcGIS (version 10.5), SDMtoolbox (version 2.10.4) and MaxEnt (version 3.4.2), the geographical distribution of moso bamboo was simulated and analyzed during three geohistorical history periods, the last interglacial (LIG) 130000–116000 years ago, the last glacial maximum (LGM) 21000 years ago and the current period (Ct, 1979–2013) (Steven et al., 2017).

A total of 1,824 million reads, 248.16 GB clean data was obtained from 193 moso bamboo by GBS ( Supplementary Table S4 ). The data size in 193 moso bamboo samples ranged from 0.90 to 2.43 GB. The clean data was aligned to the moso bamboo genome, the mapping rates were raged from 34.93% to 99.09%, with a mean alignment rate of 89.18% and an effective sequencing depth of 11.65×. Using our specified parameter settings, we identified 51,033 raw SNPs and 7,664 InDels across the 193 samples using bcftools ( Figures 1a, b ). The distribution of SNPs on 24 chromosomes of moso bamboo was examined, and the density of SNP markers was calculated. It was found that there was one SNP variant site per 3518 bp, covering the whole genome of moso bamboo. the density of SNP markers at the beginning and end of each chromosome was higher than that in the middle ( Figure 1a ). Through the statistical analysis of the SNPs annotation results, it was found that the SNPs were located in the Intergenic region ( Supplementary Table S5 ).

ADMIXTURE was used to analyze the genetic structure of the 193 samples collected from the whole natural distribution area. The number of subsets (K) was set from 2 to 10, repeat 10 runs and perform ten-fold cross validation. The optimal number of subsets was determined based on the K value corresponding to the minimum CV (Cross Validation) error. It was found that K = 4, with a CV value of 0.270, was the most appropriate number according to the CV error ( Figure 2a ). This suggests the presence of four potential subpopulations in the entire Chinese moso bamboo population. Unsupervised clustering analysis revealed that distinct subgroups such as the Hunan-Jiangxi subclass, the lower Yangtze River subclass, the southeast coast, the Zhujiang River Basin subclass in South China, and isolated samples in Xian ‘an District, Xianning City, Hubei Province. Some individual samples were found alongside others, reflecting the natural genetic structure traits of moso bamboo ( Figure 2b ).

To better understand the population structure of moso bamboo, PCA and phylogenetic tree were performed. Principal components PC1 and PC2 accounted for 52.91% and 26.77%of the variance, respectively ( Supplementary Table S6 ). The results indicated that moso bamboo was divided into two subpopulations, with 186 samples clustered on the left (96.37%), while only 7 samples from Xian ‘an (6) and Chongyang (1) in Hubei Province gathered on the right ( Figure 3a ). The phylogenetic tree of 193 moso bamboo divided into three genetic groups ( Figure 3b ). Considering China ‘s topography, the large mountains significantly affected the natural distribution of moso bamboo. Consequently, the natural distribution of moso bamboo in China was divided into three potential subpopulation structures ( Figure 3c ). The central α-subpopulation, known as the Luoxiao Mountain group, located in the middle reaches of the Yangtze River catchment area between Luoxiao-Xuefeng Mountains at the junction of Jiangxi and Hunan, north to Lushan City, Jiangxi, south to Leiyang and Dong ‘an, Hunan, mainly distributed in Mufu Mountain, Jiuling Mountain, and Wugong Mountain regions. The eastern β-subpopulation, known as the lower reaches of the Yangtze River population, located below Hukou, Jiangxi. The samples include Huoshan, Yixian, Jurong, Anji, Longyou, Longquan, Shangrao, and the population samples are almost clustered in one branch. The southern γ-subpopulation, known as the Zhujiang River Basin, located from the basins along the southeast coast of China to the Zhujiang River Estuary, extending westward to the Sanchahe area of the Xijiang River Basin.

Based on the SNPs set constructed by 37 moso bamboo populations, we conducted an analysis of their genetic diversity. The results indicated that at the species level, the HW-P (Hardy-Weinberg) value was between 0.638 and 1.000, suggesting that all moso bamboo reached the Hardy-Weinberg genetic balance (P > 0.05), with the balanced gene and genotype frequencies in each population ( Supplementary Table S7 ). At the population level, the average observed number of alleles (Na) was 1.270, the effective number of alleles (Ne) was 1.210, the expected heterozygosity (He) was 0.113, the observed heterozygosity (Ho) was 0.182, and the polymorphism information content (PIC) was 0.087. Notably, the observed Ho values in moso bamboo populations in Hubei tended to be relatively high, while those in Guangdong seem to be generally low. The nucleotide diversity (π) of moso bamboo ranged from 0.105 to 0.280, with an average of 0.129.

In order to explore the degree of genetic differentiation among moso bamboo subpopulations, the F-statistics (Fis, Fst and Fit) were calculated and evaluated for the population differentiation index among each subpopulation ( Figure 4a ). According to Robert ‘s theory (Robert and Samuel, 1978), the value of genetic differentiation coefficient distinguished three levels, 0< Fst< 0.05 (low differentiation), 0.05< Fst< 0.25 (moderate differentiation) and 0.25< Fst (high differentiation). The analysis revealed that the Fst values among the three subpopulations were all above 0.05, indicating a moderate level of genetic differentiation among the moso bamboo subpopulations. Specifically, the Fst values of the central α-subpopulation were higher compared to those of the eastern β-subpopulation and southern γ-subpopulation, indicating that the genomes of the eastern and southern populations were closer, suggesting a higher frequency of gene exchange between them, yet distinct from the central population.

The genetic and geographic distance matrix among 37 populations of moso bamboo indicated that genetic distance ranged from -0.136 to 0.798, with an average of 0.042 ( Figure 4b ). Mantel tests showed a non-significant correlation between genetic and geographic distances for the 37 moso bamboo populations (r = 0.1664, p = 0.7511; Figure 4c ). These results indicate varying different degrees of genetic variation among populations, with an overall small genetic variation.

AMOVA analysis demonstrated that 93.86% of the total variation was within populations, while only 6.14% was observed among groups ( Table 1 ). The gene flow (Nm) was less than 1.0, indicating that the genetic variation of moso bamboo population primarily originated from within the groups.

To investigate the selection pressure on moso bamboo during transmission, we used the constructed SNPs and InDels to calculate the nucleic acid diversity (π) ratio and Fst value within100 kb sliding window of the genome between different groups. We identified 693 windows based on geographical differentiation, selecting regions with significantly low or high π values in the top 5% and high Fst values in the top 5% as the windows of interest ( Figures 5a-c ). These selected regions covered a total of 3681 genes as determined by aligning these windows to the new reference genome annotation file of moso bamboo ( Figure 5d ). Each subpopulations experienced distinct selection pressures, with a significant proportion of unique genes in each, showing obvious differences in the selected genes ( Figure 5d ; Supplementary Table S8 ).

In order to explore the genes and their functions corresponding to the selection pressure of moso bamboo, GO and KEGG enrichment analysis revealed distinct gene enrichments among the three subpopulations ( Figure 6 ). In the central α-subpopulation, significant enrichment in processes such as abscisic acid metabolism, carotenoid metabolism, olefin compound metabolism, response to gibberellin, phospholipid metabolism, cell wall macromolecular biosynthesis, transmembrane transport and transmembrane transporter activity, organelles and plasma membrane, macromolecular catabolism, and endoplasmic reticulum related genes ( Figure 6a ; Supplementary Figure S3 ). The significant KEGG enrichment in processes such as genetic information, nucleocytoplasmic transport, nitrogen metabolism, arginine biosynthesis, membrane transport, alanine, aspartate and glutamate metabolism, amino acid metabolism, energy metabolism and mitochondrial biogenesis related genes ( Figure 6b ; Supplementary Figure S4 ). The eastern β-population exhibited significant enrichment in genes related to supramolecular complexes, regulation of protein transport, temperature stimulation response, abscisic acid biosynthesis and regulation, and carotenoid biosynthesis ( Figure 6a ; Supplementary Figure S5 ). The significant KEGG enrichment in genes related to brite Hierarchies, ubiquitin system, glycerophospholipid metabolism, porphyrin and chlorophyll metabolism, mitochondrial biogenesis, and transcription factors ( Figure 6b ; Supplementary Figure S6 ). In the southern γ-subpopulation, enrichment was observed in catalytic activity, RNA polymerase II transcription regulatory complex, secondary cell wall, cell surface receptor signaling pathway, cell response to endogenous stimuli, hydrolase activity genes, stimulus response-related genes, secondary cell wall genes, defense response to bacteria, and chloroplasts and chloroplast envelopes ( Supplementary Figure S3 ; Supplementary Figure S5 ). Significant KEGG enrichment was observed in genetic information, replication and repair, nucleocytoplasmic transport, amino sugar and nucleotide sugar metabolism, lipid biosynthesis proteins, circadian rhythm plant, fatty acid biosynthesis, nitrogen metabolism, arginine biosynthesis alanine, aspartate and glutamate metabolism and energy metabolism ( Supplementary Figure S4 ; Supplementary Figure S6 ).

By screening the SNP difference sites of the functional genes related to temperature stimulation response in the three subpopulations, the SNPs with stable polymorphism differences were selected, and a total of four SNP sites with good typing were screened. Finally, two nonsynonymous mutation sites (SNP: Chr24_67953493, Chr24_67953514) and one synonymous mutation site (SNP: Chr24_67953497) were screened out in PH02Gene00430 gene, and one nonsynonymous mutation site (SNP: Chr1770_920) was screened out in PH02Gene21923 gene ( Supplementary Table S9 ). Ten samples from each of the three subpopulations were randomly selected as DNA templates for PCR amplification, and the amplified products were subjected to Sanger sequencing. The results showed that HBXN and HBCY in α subpopulation and JSJR in β subpopulation were polymorphic ( Supplementary Table S9 , Supplementary Figure S7 ). The primers of three nonsynonymous SNP mutations can be used to quickly detect the ability of bamboo temperature stimulation response, which are valuable SNPs molecular marker.

According to the neutral test results of moso bamboo test materials ( Table 2 ), the Tajima ‘D and Fu ‘s Fs values for the whole population of moso bamboo, as well as for the central α-subpopulation and the eastern β-subpopulation, were negative but not statistically significant (p > 0.05) ( Figures 7a-c ). This indicates the absence of significant population expansion events in the history of Chinese moso bamboo across the population, as well as the α-subpopulation and the β-subpopulation. However, the test D and FS values of the southern γ-subpopulation were negative, with a significant p-value. The mismatch distribution analysis did not support the sudden expansion of the population (r = 0.0004, p > 0.05) ( Figures 7d-f ). These results show that the moso bamboo population and its subgroups have remained relatively stable state through their historical evolution, without experiencing obvious population expansion events. The historical population dynamics of moso bamboo showed that there were population bottlenecks effect about 8–30 ka (thousand years) and 70–120 ka years ago ( Figures 7d-f ). This stability is related to the growth habit of moso bamboo, characterized by its scattered bamboo nature and reliance on underground stems for rhizome growth every year.

In order to explore the adaptive distribution of moso bamboo in China, we meticulously selected distribution points with intervals exceeding 10 km, and finally obtained 342 valid distribution points. We then selected 19 bioclimates from WorldClim ( Supplementary Table S3 ) for MaxEnt modeling. Our approach involved utilizing 75% of the moso bamboo distribution site data for training modeling, conducting 10 repeated runs, and using the remaining 25% of the distribution site data as a validation subset. The model’s accuracy was assessed based on the highest value of the receiver operating characteristic (ROC) in the Area Under ROC Curve (AUC), providing insight into the predicted geographical distribution areas during three historical geological and climatic periods, the last interglacial period (LIG) 130 000-116–000 years ago ( Figure 8a ), the last glacial period (LGM) 21–000 years ago ( Figure 8b ), and the current period (Ct) 1979-2013 ( Figure 8c ). Across these periods, the model exhibited high accuracy, with average AUC values of 0.977 for both LIG and LGM, and 0.974 for the current period, respectively, indicating robust predictive performance.

Based on the contribution of bioclimatic variables during the three periods, the five predominant factors influencing the prediction potential distribution area of moso bamboo were the warmest season precipitation, temperature seasonality, wettest monthly precipitation, the lowest temperature of the coldest month and annual mean temperature ( Figures 8d-f ). Notably, the warmest season precipitation factor had the most substantial impact. These five key factors effectively analyze the modern natural distribution results of moso bamboo. According to the value of LV (Logistic value), the moso bamboo habitat was divided into four grades. The core suitable area with LV 0.65-1.00 ( Figure 9a ), the middle suitable area with LV 0.40-0.65 ( Figure 9b ), the low suitable area with LV 0.20-0.40, and the non-suitable area with LV 0-0.20 ( Figure 9c ). Throughout the three historical periods, the simulated distribution area of moso bamboo globally exhibited a shrinking trend, with a significant reduction in the core suitable area.

Due to the natural barriers such as mountains and rivers, historical human activities have affected the gene flow among moso bamboo populations, resulting in genetic differentiation and the development of distinct subpopulations structures. The origin of moso bamboo has limited studies because of unsuitability for preservation and lack of fossil evidence. In this study, employing PCA and phylogenetic tree analysis, combined with niche distribution, the Chinese population of moso bamboo divided into three subgroups, central α, eastern β and southern γ. Notably, populations of moso bamboo in Yunnan, Guizhou and Sichuan were found to belong to the same subgroup as those in Fujian. The central α-subpopulation exhibited the highest nucleotide diversity, indicating that it serves as a focal point for variation and diversity. According to the theory of the center of origin of cultivated plants (Hummer and Hancock, 2015), the α-subpopulation in central China may be the origin center of moso bamboo, which is consistent with the results of this study. Additionally, research utilizing SSR and SNP markers has further confirmed this result (Zhao et al., 2021; Hou et al., 2024; Jiang et al., 2017).

The genetic diversity indexes of 37 populations in China were subjected to statistical analysis. The results of He, Ho, Na, Ne, Fis and Fit at the population level of moso bamboo showed a low level of genetic variation and excess of heterozygotes within the domestic moso bamboo population. This is similar to the results of genetic diversity and genetic structure analysis of 34 populations of moso bamboo in China by SSR technology (Jiang et al., 2017).

It is well understood that genetic drift, historical bottlenecks, and climate fluctuations (such as glacial and interglacial cycles) can significantly impact population size, leading to reduced genetic diversity. These factors are particularly relevant in long-term asexual reproductive plants. The evolution of both metapopulations and stable populations is primarily driven by genetic drift and genetic bottlenecks, resulting in subpopulations with low genetic diversity (Angst et al., 2022). In our study, the genetic diversity of moso bamboo population is notably low, while heterozygosity is high, affected by genetic drift, genetic bottlenecks (8–30 ka and 70–120 ka years ago), climate change, human factors and the reproduction characteristics of moso bamboo. Moso bamboo undergoes prolonged vegetative reproduction stages and the infrequent flowering intervals, yet its flowering and seed dispersal exhibit distinct characteristics. At the subgroup level, it was observed that the middle α-subpopulation displayed the highest genetic diversity. Furthermore, the genetic differentiation between the eastern β-subpopulation and southern γ-subpopulation was relatively minor, suggesting that more frequent gene exchange between them. This is consistent with the results of resequencing of moso bamboo populations in 15 representative regions of China (Zhao et al., 2021; Hou et al., 2024).

The AMOVA analysis revealed that the genetic variation within populations exceeded that among populations, a common trend in outcrossing and asexually reproducing perennial species. These species typically exhibit high levels of heterozygosity and maintain significant genetic variation within populations (Hamrick, 1992). The differentiation level within moso bamboo population was generally moderate. The Nm value of gene flow less than 1.0 is usually considered to be the threshold for significant population differentiation (Slatkin and Barton, 1989). In this study, the Nm was 3.21, indicating that there was sufficient gene flow within the population to offset the genetic differentiation caused by gene drift (Hamrick, 1992). These genetic variations contribute to the flexibility and survival of moso bamboo populations in the face of changing environmental circumstances (Olson-Manning et al., 2012). As the most important economic bamboo species, moso bamboo is likely to experience human-mediated genotype movement between distribution regions, further influencing its genetic dynamics.

Based on the analysis of selection pressure, the genetic characteristics among the subgroups of moso bamboo were identified at the molecular level. A high proportion of unique genes were identified in the three subgroups, 48.98% in α-subpopulation, 77.14% in β-subpopulation and 81.15% in γ-subpopulation. This indicates that there were different degrees of selection pressure among subgroups, which may be related to the differences in living environment, climate and human activities. Geoffrey analyzed the selection pressure of 971 sorghum plants and mapped three plant height-related genes dw2, dw3 and SbHT9.1. These genes may represent potential variations of their main agro-climatic traits (Morris et al., 2012).

A variety of functional genes were enriched in the candidate genes screened by selection analysis, from α-subpopulation to β-subpopulation. Notably, genes related to the abscisic acid metabolism process, protein transport regulation, and supramolecular complex were enriched. Abscisic acid (ABA) is intricately involved in regulating plant dormancy, maturation and senescence, stomatal movement and various other processes (Mittler and Blumwald, 2015). Additionally, ABA plays an important role in plant responses to stressors such as drought, high salt, low temperature and pests (Hewage et al., 2020). Supramolecular complexes are essential for executing a variety of physiological functions (Croce and Amerongen, 2014). Therefore, the function of these genes may be related to the adaptability of moso bamboo, potentially influencing its response to environmental challenges and stressors.

In the process of regional differentiation from α-subpopulation to γ-subpopulation, genes related to chloroplast and chloroplast envelope were notably enriched. The area where the γ-subpopulation is located has high temperature, abundant rainfall and high humidity throughout the year. The environmental changes in the process of transmission from the central to the southern may lead to some adaptive changes in chloroplasts. Furthermore, genes associated with stimulus response and secondary cell wall were also enriched. The southern coastal areas are frequently affected by typhoons. It is speculated that these genes may regulate the growth and development of moso bamboo, enhancing the wind resistance and stress tolerance of moso bamboo.

In this study, the verification results of nonsynonymous mutation SNP loci of PH02Gene00430 and PH02Gene21923 related to temperature stimulus response in moso bamboo showed that HBXN and HBCY in α subgroup and JSJR in β subgroup were polymorphic. Functional analysis of these two genes revealed that their corresponding homologous genes were XTH15 (xyloglucan endotransglucosylase/hydrolase 15) and CRLK1 (calcium/calmodulin-regulated receptor-like kinase 1). XTHs (Xyloglucan endotransglucosylase/hydrolases) are considered chief enzymes in cell wall remodeling and play a central role in stress responses. Research shows that the content of xyloglucan components recognized by CCRC-M87/103/104/106 antibodies might be negatively related to banana chilling tolerance (Tan et al., 2024). Heat stress differentially influences the XTH expression profiles in the wheat seedlings (Iurlaro et al., 2016). And CRLK1, a calcium/calmodulin-regulated receptor-like kinase, plays an important role in regulating plant cold tolerance (Yang et al., 2010). Such as Hosta ventricosa perceives temperature changes via CRLK, and its cold tolerance is enhanced when subjected to low-temperature stress (Zhuang et al., 2021). Therefore, these two genes may be closely related to the temperature stimulus response of moso bamboo, and are potential candidate genes. The specific molecular mechanism is worthy of further study.

In the transition from the β-subpopulation to the γ-subpopulation, genes related to the defense response against bacteria were enriched in the southern group. These genes play a crucial role in enhancing plant stress resistance. The southern group, located at the southernmost end of the moso bamboo distribution, faces a growth environment characterized by high temperature, abundant rainfall, and high humidity. These conditions create a conducive breeding ground for bacteria and mildew. The genes related to defense responses against bacteria may help moso bamboo adapt to the challenging southern coastal environment.

There are few studies on the global distribution of potential moso bamboo. SDM is widely used to statistically infer the driving factors influencing species distribution under different conservation, ecological and evolutionary scenarios, aiding in solving problems related to these processes (Zimmermann et al., 2010; Neftalí, 2011). A larger sample size and wider sampling area provide more information about the relationship between species and environment, and leading to more accurate species distribution models (Koster et al., 2014; Peterson et al., 2011; Yi et al., 2018; Yan et al., 2016).

In this study, we obtained global distribution data for moso bamboo for the first time. These data not only indicate the suitable habitats for its distribution area, but also minimize the deviation of simulation results due to the lack of sampling. By analyzing the relationship between the existence probability of moso bamboo and bioclimatic variables, we obtained their corresponding contribution rates and importance. The results showed that precipitation and temperature are the main variables affecting the distribution, which is consistent with the previous research highlighting precipitation as the key meteorological constraint factor affecting the moso bamboo distribution (Shi et al., 2020). Precipitation can impact soil moisture level, with excessively dry or wet soil conditions affecting bamboo growth and development (Koster et al., 2014). The suitable average temperature for the growth of moso bamboo is 10-20°C. When temperatures drop to a critical low temperature of -14°C, moso bamboo growth and metabolic activities will stagnate, resulting in physiological drought. With the synergistic interplay of temperature and seasonal precipitation, moso bamboo can grow vigorously.

Currently, the distribution range of moso bamboo in China is smaller than the potential range we predicted, suggesting a possible expansion in the future. Model calculations indicate that the potential distribution areas of moso bamboo globally are predominantly located in Asia, North America and Europe in the northern hemisphere. In history, moso bamboo has been introduced from China to many countries, with Japan being the only hosting contiguous moso bamboo forests. The prediction results for potential distribution area highlighted the southern region of central Japan as the core suitable area for moso bamboo (Isagi et al., 2014). These Japanese areas mainly exhibit a subtropical monsoon climate, closely resembling the climatic conditions of the southern mainland China, making them conducive the natural growth of moso bamboo. This is an important reason why it can become the main economic moso bamboo species in Japan. Consequently, this study can provide valuable insights for the international introduction of moso bamboo.