1. Introduction

Woody bamboos belong to the subfamily Bambusoideae within the grass family (Poaceae), comprise approximately 135 genera and over 1700 species (Soreng et al., 2022). They are naturally distributed across all continents except Europe and Antarctica, with the highest species diversity occurring in the tropical and subtropical regions of Asia and Latin America (Gallaher et al., 2022). Dendrocalamus is a paleotropical woody bamboo genus within the Bambusa-Dendrocalamus-Gigantochola complex (i.e., BDG complex, tribe Bambuseae) (Zhou et al., 2017; Liu et al., 2020a). This genus contains many large-sized woody bamboos with about 66 species (Vorontsova et al., 2016; Soreng et al., 2022) and is widely distributed in tropical and subtropical regions of Southeast Asia. Among them, there are about 30 species in China (Li et al., 2006; Yi et al., 2008; Yang et al., 2016; Wang and Li, 2019) and 15 species in Myanmar (FAO, 2006). Two subgenera, D. subg. Dendrocalamus and D. subg. Sinocalamus, are recognized based on morphology (Li and Hsueh, 1988; Li et al., 2006). However, Dendrocalamus was suggested to be paraphyletic with two different clades in recent phylogenetic studies (Liu et al., 2020a, 2024). Three giant bamboo species of D. giganteus, D. calostachyus, and D. sinicus, all native to Myanmar and Southwest China, fell into the same subclade (D2-2) (Liu et al., 2020a). They are also very similar to one another morphologically (Gamble, 1896; Li et al., 2006). Together, they are referred to as the Dendrocalamus giganteus complex.

Dendrocalamus giganteus, commonly known as the giant bamboo, has long been documented as the world's largest woody bamboo species (Mabberley, 1997), with culms reaching up to 30 m in height and 30 cm in diameter (Keng and Wang, 1996; Dlamini et al., 2022). It is thought to be naturally distributed in southern Myanmar and northwestern Thailand (Dransfield and Widjaja, 1995; Singhal et al., 2013), and Yunnan Province, China (Li and Hsueh, 1988; Wang et al., 2016), despite of its wider cultivation in Southeast Asia and in the tropics of Africa and America. This iconic tropical bamboo species has long been associated with the daily life of indigenous people in Southeast Asia (Dransfield and Widjaja, 1995; Rao and Rao, 1999; Yang et al., 2004). It has been recognized as economically and ecologically significant (Dlamini et al., 2022) due to its multiple uses, such as in construction, scaffolding, rural housing, water pipes, buckets, boat masts, woven wares in the rural area, and as raw material for furniture and industrial paper pulp, as well as planted to prevent soil erosion (Ramanayake and Yakandawala, 1997; Bonilla et al., 2010; Razvi et al., 2011; Dlamini et al., 2022). It is also used in constructing bamboo rafts for fishing (Nyein and Mathew, 2017), particularly in southern and southwestern Myanmar, where it is extensively grown. Dendrocalamus calostachyus is found only in the border area of China and Myanmar, and the taxonomy of D. giganteus and D. calostachyus has long been regarded as difficult since the time of Gamble (1896). Meanwhile, Dendrocalamus sinicus is endemic to southwestern Yunnan and found at elevations ranging from 600 to 1500 m (Yang et al., 2018a), and it is precisely recorded as the largest bamboo species reaching up to 38 m (Ma et al., 2024; Mao et al., 2025). It is distinguished by possessing the strongest culms among all known bamboos, making it a valuable species for timber production (Sun et al., 2015). Due to its remarkable growth traits and productivity, D. sinicus holds substantial potential for use in local buildings, papermaking, and the handicraft industry (Hui et al., 2006). Recognized as a rare and ecologically significant species, it carries considerable scientific, economic, and cultural significance.

Demand for giant bamboos as sustainable non-timber forest resources is increasing due to their traditional and industrial uses as a wood substitute in Southeast Asia (Dransfield and Widjaja, 1995). Meanwhile, the conservation of bamboo germplasm is necessary due to overexploitation and genetic degradation (Loh et al., 2000). Moreover, anthropogenic activities and global climate change result in disturbed vegetation, causing habitat degradation and loss of germplasm (Ahuja and Jain, 2015). Overexploitation of bamboo's natural habitats poses a significant threat to its potentially natural populations and historically cultivated clumps, and adopting effective conservation strategies, both in situ and ex situ, is critical for effectively managing this declining resource.

Understanding population genetics is essential for the sustainable utilization and conservation of plant species (Frankham, 2005; Ma et al., 2022). Genetic diversity is crucial for the survival of species and reflects the result of long-term evolutionary processes (Reed and Frankham, 2003), including species life-history characteristics, geographic distribution, population isolation, effective population size, as well as local and global environmental conditions (Hamrick et al., 1992; Eckert et al., 2008; De Kort et al., 2021; Salgotra and Chauhan, 2023). Moreover, levels of genetic diversity reflect the availability of genetic resources that are crucial for both short-term ecological adaptations and long-term evolutionary changes (Tikendra et al., 2021a). Furthermore, such diversity helps forestry species adapt to a changing environment over thousands of years and serves as an essential foundation for evolution. Thus, it is important to understand the population genetics (e.g., genetic diversity and population structure) of the large-sized bamboo species in designing and implementing appropriate conservation strategies and utilization for them, particularly in Southeast Asia (Liu et al., 2013; Desai et al., 2015; Yeasmin et al., 2015).

Over the past two decades, investigations into bamboo's genetic diversity have been performed using numerous molecular markers, such as random amplified polymorphic DNA (RAPD) (Das et al., 2005; Ramanayake et al., 2007a, 2007b), expressed sequence tag-derived simple sequence repeat (EST-SSR) (Bhandawat et al., 2019), simple sequence repeat (SSR) (Sharma et al., 2008; Meena et al., 2019), and inter simple sequence repeat (ISSR) (Tian et al., 2012; Yang et al., 2012; Nilkanta et al., 2017; Tikendraet al., 2021b). However, due to the lack of genome-wide population investigation, understanding of population structure and genetic diversity of bamboos is far from adequate, particularly for large-sized tropical bamboos. Recently, single nucleotide polymorphisms (SNPs) based on next-generation sequencing (NGS) have become the preferred markers for determining population structure due to their abundance, stability in the genome, biallelic nature, high heritability, and ability to be accurately scored (Verma et al., 2015; Zhao et al., 2019; Yang et al., 2020a). Whole-genome resequencing (WGRS) is thus widely used to explore genetic diversity, population structure, and linkage disequilibrium in many forest species (Hou et al., 2020; Jia et al., 2020; Zhu et al., 2020; Sandercock et al., 2022; Han et al., 2024; Ren et al., 2024). It has also helped to reveal population history, particularly very low genetic diversity and historical bottlenecks, and a recent decline in Ne of the economically important Moso bamboo, Phyllostachys edulis (a temperate woody bamboo, tribe Arundinarieae) in China (Zhao et al., 2021a). Therefore, a genome-scale investigation of the genetic diversity, population differentiation, and demographic history of the three closely related giant bamboo species, D. giganteus, D. sinicus and D. calostachyus is crucial to developing and implementing effective conservation strategies for their natural and domesticated biodiversity.

Taking advantage of the chromosome-level genome assembly of Dendrocalamus sinicus (Ma et al., 2024), we used WGRS to obtain genome-wide data from 284 accessions, including population sampling across all three species of the complex from their potentially native geographical distributions in Myanmar and Yunnan Province, China, plus one D. giganteus accession from Thailand. We aimed to: (i) detect the genome-wide SNPs; (ii) evaluate the genetic diversity and population structure; and (iii) infer the population demographic history of the species complex, to devise effective conservation strategies. It is hoped that our work would provide a deeper understanding of the evolutionary history and genetic structure of the giant bamboos, serving as an important foundation for genomic resources in developing suitable conservation strategies for tropical bamboo resources in Southeast Asia.

2. Materials and methods 2.1. Sample collection and whole-genome resequencing

In the present study, we sampled 176 Dendrocalamus giganteus and 8 D. calostachyus accessions from 59 to 6 populations across the entire geographic range of these two species in Myanmar and Yunnan (China), respectively, and 100 D. sinicus accessions from 13 populations in Yunnan (China). This sampling covered the main documented natural distribution areas of these giant bamboo species, including multiple accessions collected from the remote regions in northern and northwestern Myanmar. All three Dendrocalamus species are clump-forming bamboos and we sampled individuals from different clumps at a distance of at least 100 m apart and as far as possible for each site, and at least two samples were collected for each population with a majority of three to five. Young fresh leaves were cleaned and instantly desiccated using silica gel to prevent DNA degradation and avoid any fungus contaminations throughout the fieldwork. We documented the geographical data information of each accession in Fig. S1, and their detailed information was presented in Table S1.

Whole-genome resequencing was carried out using the MGISEQ platform with a pair-end read length of 150 bp (Table S2). We filtered raw data using Fastp (Chen et al., 2018) to remove low-quality reads using default parameter for subsequent analysis.

2.2. Sequencing reads mapping, SNP calling and filtering

The chromosome-level genome assembly of Dendrocalamus sinicus (Ma et al., 2024) was employed as a reference, which has an estimated genome size of 1.4 Gb. The three species of the D. giganteus complex are all derived from the common hexaploidization event with the same chromosome numbers (2n = 70) and similar genome size (Thakur et al., 2015; Ma et al., 2024). The obtained clean sequencing data above were aligned to the reference genome using bwa-mem v.0.7.17-r1188 (Li and Durbin, 2009) with the default parameters. The conversion of mapping results from SAM to BAM file format, as well as the exclusion of unmapped and nonunique reads, was carried out using SAMtools v.1.6 (Li et al., 2009). In addition, Picard v.2.27.4 (http://broadinstitute.github.io/picard/) was used to filter out duplicate reads and the mapping rates were determined by –flagstat in SAMtools v.1.6.

We used the Genome Analysis Toolkit (GATK v.4.1.9.0) for SNP calling (DePristo et al., 2011) with the HaplotypeCaller to produce gVCF (Genomic variant call) files (-ERC GVCF). We selected SNPs using Selective Variants, applied filtering using VariantFiltration, and identified a total of 127,537,522 SNPs (Dataset 0). The parameters used for SNPs filtering were as follows: “QUAL < 30 || QualByDepth (QD) < 2.0 || mapping quality (MQ) < 40.0 || FisherStrand (FS) > 60.0 || StrandOddsRatio (SOR) > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0”. Subsequently, a total of 36,450,313 high-quality SNPs (Dataset 1) was defined by excluding those with missing genotype ≥ 20% (–geno 0.2) and a minor allele frequency (MAF) ≥ 0.05 (–maf 0.05). To eliminate the linkage between sites in Dataset 1, we employed PLINK (Purcell et al., 2007) with the parameter –indep-pairwise 50 10 0.5 and –indep-pairwise 50 10 0.2. Finally, a total of 5,850,302 and 2,141,230 independent SNPs were obtained in Dataset 2 and 3, respectively, for subsequent analyses.

2.3. Phylogenetic analysis, population stratification and PCA analysis

To examine population stratification and genetic relationships among accessions of the Dendrocalamus giganteus complex, we performed phylogenetic analysis, admixture analysis, principal component analysis (PCA), as well as analyses of population genetic diversity and linkage disequilibrium (LD), identity by state (IBS), and population demographic history.

2.3.1. Phylogenetic analysis of whole-genome SNPs

To elucidate the evolutionary relationships among the sampled accessions, a maximum likelihood (ML) tree was constructed using Dataset 2 by IQ-TREE v.1.6.10 (Nguyen et al., 2015). The general time reversible model with an ascertainment bias correction (-m GTR + ASC) was employed as the substitution model. To assess the tree's reliability, a thousand replications of the ultrafast bootstrap approach were performed. Following the phylogenetic results of Liu et al. (2020a), four congeneric species, Dendrocalamus bambusoides, D. semiscandens, D. fugongensis, and D. tibeticus, representing major subclades of the genus were used as outgroup taxa. Then, the iTOL web server was employed to apply color to the obtained ML tree.

2.3.2. Population genetic structure

Population structure analysis was conducted using ADMIXTURE v.1.3.0 (Alexander et al., 2009) employing a block relaxation algorithm. The Dataset 3 of core SNPs was utilized and the genetic ancestry of each accession was estimated in the default set by specifying the number of genetic clusters (K) ranging from 2 to 9 employing 10-fold cross-validation (CV). Then, the best value of K was determined by identifying the K value with the lowest cross-validation (CV) error estimated by ADMIXTURE. To plot the structure results, the POPHELPER v.2.3.1 (Francis, 2017) and ggplot2 packages in R were used.

2.3.3. Principal component analysis (PCA)

For principal component analysis (PCA), the top 10 principal components (PCs) were extracted using PLINK v.1.90b6.24 with the “–pca” parameter (Purcell et al., 2007). Eigenvalues were calculated to assess the percentage of variation explained by each PC. Subsequently, the first two eigenvectors obtained from the analysis were visualized using the ggplot2 package in R (Wickham, 2009). This analysis was conducted specifically on Dataset 3 and allowed us to capture the major sources of genetic variation and visualize the population structure by representing the accessions in a lower dimensional space defined by the PCs.

2.4. Divergence time estimation

The phylogenetic tree obtained above was dated using treePL (Smith and O'Meara, 2012) with two secondary calibration points. The root age was set at 7.442 million years ago (Ma), based on Liu et al. (2024). The second calibration is the divergence of Dendrocalamus semiscandens and D. fugongensis (1.138–5.579 Ma). The best smoothing parameter for the rate variation of the phylogram was determined by cross-validation from 0.00001 to 100,000, yielding an optimal value of 0.0001 (see Table S3).

2.5. Population genetic diversity and linkage disequilibrium (LD) analysis

Based on the phylogenetic tree, we classified all accessions into seven distinct phylogenetic clades (see details in the Results). The key genetic diversity parameters, including nucleotide diversity (π), population differentiation statistics (F_ST), Tajima's D, observed heterozygosity (Ho), expected heterozygosity (He), and inbreeding coefficient (F_IS) were computed using Dataset 1.

To assess nucleotide diversity (π) within these groups, we utilized VCFtools v.0.1.17 (Danecek et al., 2011) with a sliding window size of 100 kb and a step size of 10 kb, as well as to measure the genetic differentiation (F_ST) among them. Similarly, we employed VCFtools v.0.1.17 to compute Tajima's D in non-overlapping 100-kb windows, using the parameter “– Tajima's D 100,000”. Furthermore, we also utilized the same software to calculate observed (Ho) and expected (He) heterozygosity and inbreeding coefficient (F_IS) with the option “–het”.

The decay of linkage disequilibrium (LD) for each population was calculated based on the pairwise correlation coefficient (r²) statistics for genome-wide unpruned pairwise SNPs using PopLDdecay v.3.42 (https://github.com/BGIshenzhen/PopLDdecay) (Zhang et al., 2019) with default parameters. For each SNP distance, the average r² value was calculated and the LD decay patterns were visualized using R for the seven distinct groups of populations.

2.6. Identity by state (IBS) analysis

We computed pairwise genetic distance by Identity by State (IBS) using PLINK v.1.90b6.21 with the parameter “–distance square 1-ibs flat-missing” (Slifer, 2018) based on Dataset 1. The distribution of IBS values was visualized through a histogram, providing insight into the overall pattern of genetic similarity among accessions. Remarkably, low genetic distances (IBS < 0.105) can point to the occurrence of recent dispersal or introduction events. To illustrate these connections, lines were drawn on the map by linking the respective accessions.

2.7. Gene flow analysis

We used TreeMix (Pickrell and Pritchard, 2012) to infer the direction of gene flow among populations based on Dataset 3. The number of migration events (m) was tested from 0 to 10, with 10 iterations performed for each m value. The “-global” and “-se” options were used to estimate the standard errors of migration proportions, and the “-noss” parameter was applied to avoid overcorrection. The OptM v.0.1.9 package in R (Fitak, 2021) was then used to identify the optimal number of migration edges. The final tree and residual plots for the optimal migration model were visualized using the “plotting_funcs” script provided in TreeMix.

2.8. Identification of deleterious mutation

We annotated Dataset 1 using SnpEff v.4.3t (Cingolani et al., 2012). Deleterious missense variants were predicted using Sorting Intolerant from Tolerant 4G (SIFT4G) (Vaser et al., 2016). SIFT assigns a score between 0 and 1 to each amino acid substitution, and missense SNPs with scores < 0.05 were designated as deleterious, while variants with scores ≥ 0.05 were regarded as tolerated.

2.9. Population demographic history

Demographic history was reconstructed using the pairwise sequentially Markovian coalescent PSMC v.0.6.5-r67 (Li and Durbin, 2011) and sequential Markov coalescent SMC++ v.1.15.4 (Terhorst et al., 2017) models. Firstly, the PSMC model was employed to infer historical demographic changes from individuals. For each group, two accessions with a high sequencing depth (> 17×) (Table S4) were selected as representatives for PSMC analyses with 100 bootstrap replicates. PSMC was run using the default atomic time interval parameter: psmc -N25 -t15 -r5 -p “4 + 25 × 2 + 4 + 6”. Scaled estimates of time and effective population size were converted into absolute values using psmc_plot.pl.

To validate the PSMC results, we also used SMC++ (v.1.15.4) to infer the dynamics of effective population size from population-level data. For SMC++, we selected representative accessions from the seven identified phylogenetic clades, specifically those with sequencing depth greater than 15× (Table S4). This analysis was conducted based on unphased SNPs with a minor allele frequency (MAF) greater than 0.05 using Dataset 1. To delineate the mainly uncalled areas, we followed the masking step as suggested (Malaspinas et al., 2016) by the SNPable toolkit (http://lh3lh3.users.sourceforge.net/snpable.shtml). Subsequently, historical population sizes were estimated using the “SMC++ estimate (–spline cubic)” command. To ensure robust (Ne) estimates, 20 independent replicates were performed. We inferred a mutation rate (μ) of 6.5 × 10⁻⁹ mutations per site per year for this species complex, which aligns with the average estimates ranging from 5 × 10⁻⁹ to 7 × 10⁻⁹ for plant nuclear DNA sequences in previous studies (Gaut et al., 1996). The generation time (g) was set to 65 years for the D. giganteus complex (Jee et al., 2020), and the mutation rate (μ) of 6.5 × 10⁻⁹ mutations per site per year, resulting in an estimated mutation rate of 4.23 × 10⁻⁷ mutation per site per generation for these bamboo species.

3. Results 3.1. Whole-genome resequencing and SNPs calling

Totally, we conducted whole-genome resequencing of 284 accessions of the Dendrocalamus giganteus complex, including 176 identified as D. giganteus, 8 as D. calostachyus, and 100 as D. sinicus collected from across different geographical regions of China and Myanmar, including potentially native area in highly isolated human settlement (Figs. 1a and S1). The detailed localities and elevation information of all accessions was given (Fig. S4 and Table S5) and the majority of sampled D. giganteus and D. calostachyus accessions were from Myanmar, especially in the northern (NORTH), northwestern (NWEST), eastern (EAST), central (CENTRAL), southern (SOUTH) and southwestern (SWEST) regions of Myanmar. An average of approximately 24.4 Gb of clean data per sample was obtained. Approximately 98.97% of the clean reads from all accessions was aligned to the D. sinicus reference genome, resulting in an average depth of 16.86× (Fig. S2). The average mapping rates were 98.63% for D. giganteus and D. calostachyus accessions and 99.59% for D. sinicus accessions, respectively. After initial filtering, we identified a total of 127,537,522 single nucleotide polymorphisms (SNPs). Subsequent filtering steps led to the identification of a core set of 36,450,313 SNPs with a minor allele frequency (MAF) ≥ 0.05 and missing genotype ≤ 20% (Fig. S3). Overall, the sequencing produced high-quality data, providing a solid basis for subsequent analyses.

3.2. Morphological characters, phylogenetic relationships, and population structure

Following current taxonomy, we identified the accessions of the Dendrocalamus giganteus complex into three species, i.e., D. giganteus, D. calostachyus, and D. sinicus based on field observations and floristic accounts in Bambuseae of British India (Gamble, 1896) and Flora of China (Li et al., 2006). Phylogenetic analysis classified the 284 accessions into seven distinct clades (Fig. 1b and e). Group Ⅰ to Group Ⅴ comprise of accessions from D. giganteus and D. calostachyus, while all accessions of Groups Ⅵ and Ⅶ (49 and 51 accessions, respectively; see Table S5) were identified as D. sinicus. The two seemingly isolated accessions, MZW048 and MZW181, are morphologically similar to other accessions of Group Ⅰ. Accessions identified as D. calostachyus were mostly clustered in Group Ⅲ, but there are three samples (MZW182, MZW188, and MZW189) mixed with four accessions of D. giganteus in Group Ⅳ. Group Ⅰ mainly comprised accessions from the northern Myanmar (20 from NORTH and three from NWEST), while accessions from the southern, southwestern, and central regions of Myanmar, were mainly found in Group Ⅱ (44 SOUTH, 57 SWEST, seven CENTRAL, one NORTH, and one from Thailand). Furthermore, the eastern accessions were mainly found in Groups Ⅲ (five EAST) and Ⅳ (six EAST and one CENTRAL). In Group Ⅴ, all the 25 accessions collected from China were grouped with the northern and northwestern accessions (four and 10, respectively) of Myanmar. Morphologically, the diagnostic characters are similar between D. giganteus and D. calostachyus, except that the culm of D. calostachyus appears with silver hairs and is more thickened than that of D. giganteus, the culm leaf sheath blade is erect, and the leaf sheath has long ciliate at margins (Fig. 1c). Generally speaking, the morphology of Group Ⅰ (including other two isolated accessions) to Group Ⅴ is very similar, except that some individuals of Groups Ⅲ and Ⅳ have a slightly different culm character, such as slightly shorter internode, culms with silvery hairs, thicker culm wall and erect culm leaf blade (Fig. 1c). However, the morphology of Groups Ⅵ and Ⅶ of D. sinicus is different from D. giganteus by shorter internode, delayed falling culm leaf sheath, and erect sheath blade (Fig. 1c and Table S6). Interestingly, the accessions of Groups Ⅵ and Ⅶ correspond to the two natural types with straight or bending culms in D. sinicus, respectively.

For genetic structure analysis, the number of identified populations (K) varied from 2 to 9 (Figs. 1e and S6; Table S7) with the optimal K value at 7 (Fig. S7), also aligning well with the classification of the sampled accessions into seven groups. More specifically, this analysis revealed that the deepest split start occurred at K = 4 between Group Ⅰ and Group Ⅱ (Fig. S6). These two groups also exhibited relatively pure ancestral contributions. Moreover, clear genetic differentiation was observed between the populations of Groups Ⅵ and Ⅶ up to K = 6 and they showed the same pure ancestral contribution at K = 6 (Fig. 1e). However, it is evident that all accessions from Group Ⅳ exhibited an admixture of genetic components, exhibiting a more complex pattern when the K values were larger (Fig. S6). In addition, a few accessions in Group Ⅴ displayed admixed genetic components. At K = 7, Group Ⅱ displayed increased diversity and showed admixture with two different genetic components (Fig. 1e). Furthermore, the genetic components of Group Ⅱ exhibited some kinds of complexity, particularly at K = 9 (Fig. S6). This suggested a higher genetic diversity of Group Ⅱ, which was consistent with the analysis of nucleotide diversity (π) as detailed below (Fig. 2c and Table 1).

PCA analysis consistently showed a clear separation of Group Ⅰ and Group Ⅱ from the others based on PC1 and PC2 (Figs. 1d and S8; Table S8). There appeared to be a relatively close genetic relationship between Group Ⅲ, Group Ⅳ, and Group Ⅴ. However, Group Ⅵ and Group Ⅶ showed very close genetic relationships based on PC1, PC2, and PC3, which were somewhat different from the structure and phylogenetic analyses (Fig. S8).

The chronogram inferred from whole-genome sequences indicated that the beginning of the divergence of the Dendrocalamus giganteus complex occurred ~4.64 Ma during the early Pliocene (Fig. 1c). The seven groups then diverged sequentially from 4.64 to 0.58 Ma, corresponding to the early Pliocene to Pleistocene period. Notably, all the five groups of D. giganteus and D. calostachyus (Groups Ⅰ–Ⅴ) diverged before 1.78 Ma, whereas D. sinicus (Groups Ⅵ and Ⅶ) originated after this time and diversified around 0.58 Ma during the middle Pleistocene.

3.3. Genetic differentiation of inferred populations

Population parameters including Ho, He, F_IS, π, F_ST, and Tajima's D were calculated for each group to estimate the patterns of genetic diversity in the populations of the D. giganteus complex (Fig. 2a–d and Table 1). The expected heterozygosity (He) ranged from 0.3087 (Group Ⅴ) to 0.3094 (Group Ⅶ), while the observed heterozygosity (Ho) varied from 0.1081 (Group Ⅵ) to 0.4387 (Group Ⅱ). Group Ⅲ exhibited the lowest observed heterozygosity (Ho = 0.3443) among the D. giganteus and D. calostachyus populations. In contrast, there was relatively little variation in expected heterozygosity (He) among all groups, with an average of 0.3090. Notably, all groups exhibited Ho values exceeding their corresponding He values, except for the D. sinicus populations (Groups Ⅵ and Ⅶ) (Fig. 2b).

The inbreeding coefficient (F_IS) ranged from −0.4191 (Group Ⅱ) to 0.6505 (Group Ⅵ). Notably, all Dendrocalamus giganteus and D. calostachyus populations (Groups Ⅰ to Ⅴ) exhibited negative inbreeding coefficient (F_IS) values (Fig. 2a), while the D. sinicus populations showed positive (F_IS) values, indicating heterozygote excess in D. giganteus and D. calostachyus. It is important to note that some degree of variation exists in both Ho and F_IS among these populations (Fig. 2a and b; Table 1). Moreover, the nucleotide diversity (π) was relatively low among all groups, ranging from 1.65 × 10⁻³ to 6.79 × 10⁻³. Group Ⅱ exhibited the highest genetic diversity (π = 6.79 × 10⁻³), followed by Group Ⅳ (π = 6.74 × 10⁻³), Group Ⅴ (π = 6.68 × 10⁻³), Group Ⅰ (π = 5.69 × 10⁻³), and Group Ⅲ (π = 5.47 × 10⁻³), while Groups Ⅵ (π = 1.65 × 10⁻³) and Ⅶ (π = 1.84 × 10⁻³) exhibited the lowest (Fig. 2c and Table 1). Together, these results suggested the highest genetic diversity in Group Ⅱ while the lowest for the D. sinicus populations.

The positive Tajima's D values may signify population bottlenecks and/or balancing selections within a population (Ching et al., 2002; Liu et al., 2020b; Miao et al., 2021; Pandey et al., 2021; Tandoh et al., 2021) and all these values observed here are positive (Fig. 2d and Table 1). Interestingly, the Group Ⅱ had the highest Tajima's D value (4.2167) among these groups, indicating that it may lack rare alleles.

The mean F_ST values among the seven groups ranged from 0.15 to 0.64 (Fig. 2c). Particularly, Group Ⅲ and Group Ⅳ accessions exhibited the lowest F_ST value (0.15) among all the comparisons, suggesting possible frequent gene exchange between them. On the other hand, the F_ST value between Group Ⅲ and Group Ⅵ accessions was the highest (0.64). This may be attributed to factors such as geographical isolation and limited gene flow among them, thus resulting in substantial genetic differences.

3.4. Linkage disequilibrium (LD) analysis

By utilizing whole-genome SNPs, we investigated the decay of LD across different distances (Fig. 2e) and revealed similar decay patterns among Dendrocalamus giganteus, D. calostachyus and D. sinicus. The decay of LD reached half-maximum average r² at 0.06 kb (r² = 0.369) for Group Ⅵ, 0.06 kb (r² = 0.348) for Group Ⅶ, 0.06 kb (r² = 0.325) for Group Ⅱ, 0.08 kb (r² = 0.343) for Group Ⅰ, 0.2 kb (r² = 0.364) for Group Ⅲ, 0.3 kb (r² = 0.294) for Group Ⅴ, and 3.8 kb (r² = 0.293) for Group Ⅳ, respectively (Tables S9 and S10). Group Ⅱ showed the fastest LD decay, which could be correlated with the high level of nucleotide diversity observed, and Group Ⅳ exhibited the slowest LD decay.

3.5. Genetic distance analysis

To gain insights into the dispersal history of Dendrocalamus giganteus, together with D. calostachyus, with a relatively wide range of distribution, we conducted pairwise Identity by State (IBS) genetic distance among their accessions. Our analysis revealed a general correlation between genetic distance and geographical separation. Notably, we observed an extreme case where accessions exhibited very small IBS values (IBS < 0.105) (Fig. 3a and Table S11). These groups consisted of closely related accessions, indicating recent dispersal or introduction events, which were primarily found within China (accession pairs represented by green lines in Fig. 3b). Furthermore, we revealed additional pairs of accessions with an IBS < 0.105 in China and Myanmar, particularly in the NORTH and NWEST regions. This suggests a potential gene flow recently and connectivity between geographically distant populations. Interestingly, the pairs with the largest IBS values (IBS > 0.27, represented by red lines in Fig. 3b) were also found between the NORTH, NWEST, SOUTH, SWEST, and CHINA regions. Some such pairs were also identified between the NORTH and SOUTH/SWEST regions (Fig. 3b). These findings again reflected a high level of genetic differentiation among Groups Ⅰ (NORTH and NWEST) and Ⅴ (CHINA), Groups Ⅰ and Ⅱ (SOUTH and SWEST), and Groups Ⅱ and Ⅴ. This suggests a potential common ancestry or historical gene flow, possibly via an introduction or migration event.

3.6. Gene flow

We further employed TreeMix to investigate historical gene flow and migration among the seven phylogenetic groups. The optimal number of migration edges was identified as two (m = 2), pointing to the occurrence of two migration events (Fig. S9). Using this parameter, we constructed the maximum-likelihood tree and residual heat map, revealing gene flow from Group Ⅰ (NORTH and NWEST) to Group Ⅳ (Eastern population), as well as from Group Ⅲ (Eastern population) to Group Ⅴ (China and northwestern populations) (Fig. 4). These results were largely consistent with the patterns inferred from population structure analysis.

3.7. Identification of deleterious mutations

We evaluated the genetic load of each population to examine its potential correlation with nucleotide diversity, which represents the accumulation of deleterious alleles that can disrupt gene function and reduce overall individual fitness (Robinson et al., 2016). The Dendrocalamus sinicus groups (Groups Ⅵ and Ⅶ) had a higher genetic load than the D. giganteus and D. calostachyus groups (Groups Ⅰ–Ⅴ) (Table 2). Among the D. giganteus and D. calostachyus groups, the highest proportion of deleterious variants was found in Group Ⅴ while the lowest in Group Ⅲ. In addition, nucleotide diversity (π) was negatively correlated with genetic load (R = −0.79, P = 3.6e-02), depending on the absolute number of deleterious mutation sites. These results suggest that selection efficiency is reduced in regions with low genetic diversity.

3.8. Inference of demographic history

We observed a notable lack of genetic diversity and a remarkably positive Tajima's D value in populations of the Dendrocalamus giganteus complex (Table 1). Thus, we conducted the demographic analysis using the combination of pairwise sequentially Markovian coalescent (PSMC) (Li and Durbin, 2011) and sequential Markov coalescent (SMC++) (Terhorst et al., 2017) to reveal the historical population changes. First, we used the pairwise sequentially Markovian coalescent (Li and Durbin, 2011) to reconstruct the demographic history of each lineage and investigate long-term population changes. PSMC analysis revealed similar demographic histories for the populations of the Dendrocalamus giganteus complex, with effective population sizes (Ne) peaking at approximately 190 thousand years ago (Kya) and 150 Kya in the D. giganteus and D. calostachyus groups (Groups Ⅰ–Ⅴ), respectively, and around 270 Kya in the D. sinicus groups (Figs. 5a and S10). This was followed by a pronounced decline in Ne during the Last Glacial Period (LGP; ca. 115–11.7 Kya) for all populations with a notable bottleneck around 27–14.5 Kya during the Last Glacial Maximum (LGM, ca. 19.0–26.5 Kya) for all D. giganteus and D. calostachyus groups (Sun and An, 2005; Clark et al., 2009; Zhao et al., 2021b). Subsequently, Ne of Groups Ⅵ and Ⅶ gradually increased by reaching approximately ~1000 around 17 and 11 Kya, respectively. After that, both groups experienced decline again in Ne until about 10 Kya.

A comparable demographic history was also supported by SMC++ analysis. Similarly, SMC++ analysis also revealed a rapid decline in the Ne of Dendrocalamus giganteus and D. calostachyus groups between approximately ~100 Kya and ~7 Kya, from the LGP to the mid-Holocene (Fig. 5b and Table S12). Subsequently, Group Ⅰ, Group Ⅱ, and Group Ⅴ underwent a bottleneck event of approximately 14.5 Kya, 11.4 Kya, and 14.8 Kya after the LGM. Additionally, Group Ⅲ and Group Ⅳ experienced a more severe population bottleneck around 7 Kya during the mid-Holocene. After that, there was a gradual increase in Ne to a maximum at about 0.2–0.03 Kya for all groups, which may be related to the initiation and increased cultivation of D. giganteus in southern and southwestern Myanmar (Gamble, 1896) for the utilization of tiger mouth nets (Kyarr Phong) at that period (Nyein and Mathew, 2017).

Similarly, the effective population size (Ne) of Dendrocalamus sinicus showed a sharp decline between approximately 110 Kya and 27 Kya during the LGP (Fig. 5b and Table S12). Following this bottleneck, Ne gradually recovered around 7 Kya, after which it declined again until ~0.14 Kya. These population fluctuations were likely caused by historical climate changes after the mid-Holocene period (Hao et al., 2021), as well as by human disturbances (Qin et al., 2023).

4. Discussion 4.1. Phylogeny and population structure of the Dendrocalamus giganteus complex

In this study, we conducted whole-genome resequencing of giant bamboos, including 284 representative accessions sampled from different geographical regions in Yunnan Province of China and Myanmar. The percentage of clean reads mapped to the reference genome of Dendrocalamus sinicus ranged from 92.84% to 99.86% (Table S2), indicating high-quality data generation. However, the choice of reference genome will certainly have an impact on the results despite that it is a regular approach using reference genome from different taxon in studies of non-model species (Akopyan et al., 2025). Here, D. sinicus has a very close phylogenetic relationship with D. giganteus and D. calostachyus, forming a species complex. The average mapping rates were nearly identical among the three species, and clear population genetic patterns were observed. Taken together, these results indicate that the use of D. sinicus as reference genome is acceptable for comparative analysis here.

For the Dendrocalamus giganteus complex, current taxonomy named them as D. giganteus, D. calostachyus, and D. sinicus, all giant bamboos higher than 30 m with D. sinicus as the precisely documented largest grass species. The 284 accessions covering the major known geographical distributions of these species in Myanmar and Yunnan Province could be grouped into seven distinct clusters by phylogenetic analysis. However, this grouping did not consistently correspond to their geographical origin, particularly for Group Ⅰ, Group Ⅱ, and Group Ⅴ. Moreover, accessions of D. calostachyus are mixed together with those of D. giganteus and in fact these two species are morphologically very similar. Taken together, these results suggest that D. calostachyus could be merged into D. giganteus as one species, which is then sister to D. sinicus. The accessions of D. sinicus are clustered into two groups: straight-culm type (Group Ⅵ) and bending-culm type (Group Ⅶ) (Fig. 1c), consistent with previous findings (Yang et al., 2018a). In Group Ⅴ, all D. giganteus samples from China were clustered together with northern and northwestern samples from Myanmar, probably related to anthropogenic activities and adjacent distribution. It is noteworthy that the accessions from northwestern Myanmar were collected in similar elevation range and their morphological characters were almost identical to those of the Chinese accessions (Fig. 1c and Table S5).

Population structure analysis revealed that an optimal number of groupings (K = 7) also resulted in seven major genetic clusters across all accessions of the Dendrocalamus giganteus complex (Fig. 1e). For D. giganteus, the accessions from Group Ⅰ and Group Ⅲ exhibited relatively pure genetic compositions without mixture from other groups. Conversely, accessions from Group Ⅱ, Group Ⅳ, as well as some accessions from Group Ⅴ, displayed an admixture of genetic components. This suggests that northern and northwestern Myanmar are the hotspot of genetic diversity and a possibly center of origin of the giant bamboos. However, further population genomic study with adequate sampling from Thailand and adjacent Laos is needed to verify this hypothesis. Furthermore, all D. sinicus populations displayed distinct and pure genetic compositions. The results of the principal component analysis (PCA) are consistent with those obtained from the phylogenetic tree and ADMIXTURE analysis, except for two groups of the D. sinicus populations (Groups Ⅵ and Ⅶ) (Fig. 1d and e).

4.2. Population genetic diversity and genetic differentiation

The level of genetic diversity reflects a species' ability to adapt to environmental changes (Reed and Frankham, 2003). In this study, Group Ⅱ exhibited the highest genetic diversity, while Groups Ⅵ and Ⅶ exhibited the lowest within the sampled populations. Compared with other woody plants, the mean nucleotide diversity value of Dendrocalamus giganteus and D. calostachyus populations (π = 6.27 × 10⁻³) was higher than that of the Moso bamboo (π = 7.0 × 10⁻⁴) (Zhao et al., 2021a), peach (π = 1.5 × 10⁻³) (Verde et al., 2013), cassava (π = 2.6 × 10⁻³) (Kawuki et al., 2009), and apricot (π = 2.7 × 10⁻³) (Li et al., 2020), but lower than that of a temperate bamboo, Chimonobambusa tumidissinoda (π = 0.116) (Ye et al., 2024) and date palm (π = 9.2 × 10⁻³) (Hazzouri et al., 2015).

Differentiation coefficients of pairwise (F_ST) are commonly used to assess population differentiation, where values between 0.00 and 0.05 indicate minimal divergence and values greater than 0.15 indicate high divergence (Wright, 1943, 1978; Slatkin, 1985; Nassiry et al., 2009; Yang et al., 2020b; Shu et al., 2021). The mean F_ST values among groups ranged from 0.15 to 0.64, indicating moderate to very high population differentiation (Fig. 2c), which is consistent with our population structure results. Overall, our investigation revealed low intrapopulation genetic diversity and high interpopulation differentiation in the Dendrocalamus giganteus complex (Table 1 and Fig. 2c), consistent with an early study that used ISSR markers on D. giganteus populations sampled just from Yunnan, China (Tian et al., 2012). Other studies using RAPD analysis of D. giganteus from Sri Lanka (Ramanayake et al., 2007b) and nSSRs of D. hamiltonii from India also found low genetic diversity (Meena et al., 2019). In addition, a high level of genetic differentiation has been previously reported for the populations of D. sinicus in China, using simple sequence repeat (SSR) markers (Yang et al., 2018a). Long-lived outcrossing species with wide geographical distributions typically maintain high genetic diversity while with lower genetic differentiation (Hamrick and Godt, 1996). However, in the case of woody bamboos, populations may undergo genetic differentiation over time due to various issues, such as flowering and breeding behavior, population size, and habitat fragmentation (George et al., 2009). The observed high level of genetic differentiation among populations may be due to long vegetative phases and widely separated, unsynchronized flowering events throughout their lifespan, observed in some tropical woody bamboo species, including D. giganteus and D. sinicus (Ramanayake and Yakandawala, 1998; Chen et al., 2017). Additionally, D. giganteus exhibiting sparse flowering (Jee et al., 2020) or low synchrony (Ramanayake and Yakandawala, 1998) may represent genetically diverse materials, particularly under cultivation. Furthermore, other factors, such as anthropogenic pressure, environmental changes, habitat fragmentation, limited gene flow, and selection for local ecological adaptation, may also contribute to reduced diversity within populations and indirectly increase F_ST values (Cruickshank and Hahn, 2014).

Consistent negative inbreeding coefficients (F_IS) observed across the D. giganteus complex populations suggest heterozygote excess, reflecting low levels of inbreeding and moderate genetic connectivity. This pattern is likely due to the distinctive reproductive strategy of these bamboos, which is characterized by an exceptionally long–life cycle of infrequent flowering and poor and limited seed production (Jee et al., 2020). In addition, anthropogenic activities, such as the widespread vegetative propagation of genetically diverse clones in cultivated clumps (e.g., via stem cuttings; Dransfield and Widjaja, 1995), may further reinforce this genetic structure.

Reduced genetic diversity is often shaped by multiple interacting factors, including limited distribution range, breeding system, longevity, evolutionary history, and modes of seed dispersal (Hamrick, 1990; Imai et al., 2016; Wang et al., 2025). In general, species exhibiting inbreeding tend to possess reduced levels of genetic variation. Dendrocalamus sinicus populations exhibited positive F_IS values, indicating a predominance of selfing or inbreeding. Previous studies have reported that many D. sinicus flowering clumps exhibit asynchronous flowering behavior, which may reduce pollination efficacy and result in a lack of seed set (Chen et al., 2017). Additionally, seed dispersal is limited to a relatively short distance (Gu et al., 2012). Together with its narrowed distribution, these factors likely act synergistically to reduce the genetic variation in D. sinicus.

Moreover, the accumulation of deleterious mutations increases the genetic load (Agrawal and Whitlock, 2012; Chen et al., 2020), which can elevate the risk of extinction in small populations (Lande, 1994; Lynch et al., 1995; Dussex et al., 2021). Our study revealed a negative correlation between deleterious mutations and nucleotide diversity as expected. In addition, we found low genetic diversity across all populations, suggesting that this may be related to the accumulation of deleterious mutations.

4.3. Linkage disequilibrium and genetic distance

Generally, domesticated clumps can have higher levels of linkage disequilibrium (LD) than wild populations, and species with faster LD decay are thus considered more archaic (Xia et al., 2015; Moyers et al., 2018; Aleksandrov et al., 2021). Here, our analysis revealed that the Group Ⅱ and Group Ⅰ populations displayed faster LD decay among the Dendrocalamus giganteus complex populations (Fig. 2e). This finding supports our hypothesis that the D. giganteus complex might have an origin from the northern and northwestern Myanmar and spread to the southern and southwestern Myanmar, which is also consistent with the results of the Identity by State (IBS) genetic distance analysis (Fig. 3b). Additionally, the Group Ⅱ and Group Ⅴ populations showed relatively high genetic diversity (π value) compared to the remaining groups, indicating that the early introduction or cultivation of the bamboo species may have occurred from northern and northwestern regions of Myanmar into China and spread to southern and southwestern regions of Myanmar. However, due to the anthropogenic activities and existing sampling in our study, we tend to keep this observation open, awaiting more thorough sampling to fully uncovered the origin and evolutionary routines of the species complex.

4.4. Evolutionary and demographic history of the Dendrocalamus giganteus complex

In representing the possible origin center of the D. giganteus complex, divergence time analysis revealed that Group Ⅰ (northern and northwestern populations) diverged earliest among all populations around 4.64 Ma in the Pliocene. This inference is corroborated by the TreeMix analysis, which shows that historical gene flow occurs from Group Ⅰ (northern and northwestern populations) to Group Ⅳ (eastern population). The divergence of the D. giganteus and D. calostachyus populations mainly occurred during the Pleistocene and the split between the Group Ⅵ and Group Ⅶ lineages of D. sinicus happened last around 0.5 Ma (Fig. 1c). The pronounced climatic and geographic fluctuations during the Pliocene–Quaternary transition (An et al., 2001) likely played a critical role in shaping the lineage divergences within the D. giganteus complex. It is well recognized that global climatic fluctuations over millions of years, particularly from the Pliocene to the Quaternary, have been key drivers of population dynamics and lineage divergence in many tree species (Zhu et al., 2020; Feng et al., 2024; Ren et al., 2024). In addition, these climatic shifts have also influenced the demographic history of the D. giganteus complex.

Investigating the population history and dynamics of a species provides valuable insights into its origin, demographic expansion, and the processes driving population differentiation (Yang et al., 2018b; Ma et al., 2021). In this study, we employed PSMC and SMC++ analyses to infer the demographic history of the D. giganteus complex. Both PSMC and SMC++ analyses revealed a consistent demographic history pattern for the D. giganteus complex and identified several historical events that may have influenced the current demography of these species. A sharp decline in Ne for all the D. giganteus populations between approximately ~100 Kya and ~7 Kya (Fig. 5a and b) was observed. This decline may have been caused by the low temperatures and abrupt climate change during Last Glacial Period (Corrick et al., 2020), as well as anthropocenic impact since the early Holocene (Qin et al., 2023). Moreover, two bottlenecks around 110–27 Kya (during the Last Glacial Period) and around ~7 Kya (after the mid-Holocene period), respectively, were revealed in D. sinicus. Recent findings on Moso bamboo in subtropical China revealed a rapid decline in Ne during the last glacial period (Corrick et al., 2020). The large-sized bamboos within the tropical genus Dendrocalamus in southern and southwestern China may share similar evolutionary histories. Furthermore, this significant reduction in population size is likely to be partially linked to human activities. Pollen records also provided evidence of anthropogenic disturbance in the Cangshan Mountains in Western Yunnan, around 6370 years ago, coinciding with the deforestation for shifting agriculture (Shen et al., 2006).

As a result, the low population-level genetic diversity and intense population size bottlenecks in the populations of the Dendrocalamus giganteus complex, which are closely related to climate change and anthropogenic disturbance, highlight the need for effective conservation approaches and promoting the sustainable usage of large-sized bamboo species under semi-natural conditions, mostly maintained in house-hold cultivation, to conserve their genetic diversity.

5. Conclusion

Our population genomic analyses highlight the historical bottleneck events, low genetic diversity, and strong genetic differentiation among populations of the Dendrocalamus giganteus complex, and the recently reduced effective population size of these important bamboos. Therefore, to prevent further damage to the genetic diversity, conservation strategies should emphasize the protection of populations, specifically those with high genetic diversity (e.g., Group Ⅱ and Group Ⅴ). Additionally, the presence of abundant specific germplasm provides a valuable resource for the establishment of natural Dendrocalamus bamboo forests under a changing environment. Moreover, considering the vegetative propagation of the D. giganteus complex, the collection of bamboo genetic resources would have little effect on its original distribution. Therefore, intensified efforts for ex situ conservation and utilization of germplasm from different populations, particularly the hotspot region of northern and northwestern Myanmar are recommended. On the other hand, D. sinicus, bamboo species endemic to southern and southwestern Yunnan, exhibited a relatively low level of genetic diversity than populations of D. giganteus and D. calostachyus, as well as a significantly declining population size from the historical period. These findings underscore the urgent need to develop more effective conservation strategies for these giant bamboos. Overall, this study contributes to a comprehensive evolutionary and genetic framework for bamboo research and provides important insights for implementing appropriate conservation strategies for tropical woody bamboos in generally and large-sized woody species of the Bambusa-Dendrocalamus-Gigantochloa in particularly.

Acknowledgments

We are grateful to the staff of the Southeast Asia Biodiversity Research Institute, Chinese Academy of Sciences, Yezin, Nay Pyi Taw, Myanmar, the Director of the Forest Research Institute for permission to conduct fieldwork, U Aung Zaw Moe for his useful help in plant identification, and the local staff of MONREC for their help during fieldwork. We would also like to thank Ms. Liang-Min Liu and Dr. Sheng-Yuan Qin for their advice on data analysis and all the colleagues in the research group who helped with this study. This study was supported by the National Natural Science Foundation of China (32120103003), Yunnan Fundamental Research Projects (202401AS070082 and 202101AT070175), and a grant for Reserve Talents for Young and Middle-aged Academic and Technological Leaders in Yunnan Province, China (202105AC160022). Computational resources were supported by the Scientific Data Center, Kunming Institute of Botany, CAS.

CRediT authorship contribution statement

May Zin Win: Writing – original draft, Resources, Investigation, Formal analysis, Data curation. Shuang-Xiu Xu: Data curation. Shu-Yang Gao: Investigation. Jing-Xia Liu: Investigation. Zu-Chang Xu: Investigation. Cen Guo: Investigation. Yun-Long Liu: Investigation. Peng-Fei Ma: Writing – review & editing, Supervision, Methodology, Resources, Visualization, Conceptualization. De-Zhu Li: Writing – review & editing, Supervision, Resources, Methodology, Funding acquisition, Visualization, Conceptualization.

Data accessibility statement

All the newly generated bamboo genome resequencing data in this study were deposited to National Genomics Data Center (NGDC) database under the accession PRJCA055394.

Declaration of competing interest

The authors Peng-Fei Ma and De-Zhu Li are editors for Plant Diversity and were not involved in the editorial review or the decision to publish this article. The other 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.2025.12.018.