In this study, we used the young and healthy leaves of D. farinosus, B. changningensis (R-WTS-BC-005-2015) and B. rigida collected from Changning, Sichuan Province, China (28°29’N, 104°58’E). For molecular studies 21 B. changningensis individual plant samples collected from three sites were used, whereas three individual plants each for B. rigida and D. farinosus were used. Further, 15 individuals per species were used for analysing fibre length, fibre width, fibre wall cavity ratio, pulping yield, cellulose content, stem height and stem diameter.

Total genomic DNA was extracted using the CTAB protocol (Doyle, 1991). The total genomic DNA was used as a template for PCR amplification, and the GBSSI gene was amplified and sequenced. The primer sequences of GBSSI partial regions mentioned by Ye (2010) for PCR amplification ( Table S1 ). The PCR amplification was carried out in a 50 μl reaction volume containing 25 μl Taq master mix (2x Specific Taq Master Mix, 250 units, novoprotein, Suzhou, China), 22 μl ddH₂O, 100 ng total DNA and 2 umol/L of each primer. The PCR program followed was initial denaturation of 95 °C for 5 min, followed by 35 cycles of 95 °C for 30 s (denaturation), 65 °C for 30 s (annealing) and 72 °C for 40 s (extension) and a final extension at 72 °C for 10 min. The PCR products were separated in the 1% agarose gels, and the PCR products were extracted from the gel using SanPrep Column DNA Gel Extraction Kit (Simgen, Hangzhou, China) and sequenced by Sanger sequencing.

The total RNA was isolated by Polysaccharides & Polyphenolics-rich Plant Total RNA Kit (Simgen, Hangzhou, China) from the bulk tissues of B. changningensis leaves, stem, root and apical meristem as per the manufacturer’s instructions. The quality and quantity of RNA were detected using the NanoPhotometer spectrophotometer (Implen, CA, USA) and electrophoresis (1% agarose gel). The libraries were synthesized utilizing the NEB Next® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA). These constructed libraries were sequenced using the Illumina HiSeq forum to obtain paired-end reads (150 bp). After quality control (QC), transcriptome data were further handled and assembled following the procedure described in (Grabherr et al., 2011) using Trinity v2.8.4. The Assembled results by Trinity were processed using Corset as described in (Davidson and Oshlack, 2014).

The SSR motifs from the transcriptomic data were mined using the online software MISA (Beier et al., 2017) with 10, 6, 5, 4, 3, and 3 values for mono-, di-, tri-, tetra-, penta- and hexa- SSR motifs, respectively ( Table S2 ). The SSR primers were designed by Oligo7 (https://www.oligo.net/downloads.html) (Rychlik, 2007) and synthesized by the Tsingke company (Beijing, China). PCR amplification was performed as described above. The PCR products were separated in the 5% agarose gels, and the PCR products were extracted from the gel using SanPrep Column DNA Gel Extraction Kit (Simgen, Hangzhou, China) and sequenced by Sanger sequencing.

The total DNA of B. changningensis was used to develop the whole genome reads through the Illumina NovaSeq PE150 platform (Novogene Bioinformatics Technology Co. Ltd, Beijing, China). The complete chloroplast genome was assembled by mapping the whole genome reads to reference genome B. emeiensis (HQ337797) using MITObim v1.8 (Hahn et al., 2013) and annotated with the PGA (Qu et al., 2019).

The phylogenetic tree for GBSSI partial regions were constructed using PAUP 4.0a software (Swofford, 2008). Twenty-five GBSSI sequences from different bamboo species were downloaded from the NCBI database and used in this study ( Table S3a ). The DNA sequences were aligned using MEGAX (Kumar et al., 2018), and highly aligned sequences were used to generate the phylogenetic tree using the maximal parsification method in PAUP 4.0a software as per the parameters described in (Winkworth and Donoghue, 2004).

The 24 cp genomes of different bamboo species were downloaded from the NCBI database ( Table S3b ) and all the genes were extracted and concatenated to get the final sequences matrix to construct a phylogenetic tree. The sequences were aligned by MAFFT (Katoh and Standley, 2013) and the phylogenetic tree was constructed by MEGA X software using the Neighbor-Joining method (Kumar et al., 2018). A bootstrap value of 1000 replicates was used to assess the statistical significance.

In this study, we analyzed morphological characteristics such as fibre length, fibre width, fibre wall cavity ratio, pulping yield, cellulose content, stem height, stem diameter, sheath and leaves for B. changningensis, D. farinosus and B. rigida ( Table 1 and Figure 1 ). As expected, B. changningensis shows the intermediate characteristics of D. farinosus and B. rigida in biological traits such as the sheath size, leaf width, and cellulose content. The stem height and diameter were more similar to the putative parent D. farinosus. Interestingly, B. changningensis also contains superior morphological characteristics like longer fibre length, fibre width, lower fibre wall cavity ratio and higher pulping yield than putative parents D. farinosus and B. rigida.

We have developed pair-end transcriptomic data of B. changningensis using the Illumina HiSeq platform to obtain genomic information to develop the DNA markers. A total of 23,270,949 reads containing 6,961,574,574 bp were retained after quality trimming. After de novo assembly, a total of 154,983 contigs and 100,364 unigenes were obtained ( Tables S4 , S5 ). The contigs were subjected to SSR mining through MISA software and identified a total of 3,732 SSR loci, including di-, tri-, tetra-, penta- and hexa-nucleotide motifs. The di- and tetra-nucleotide repeats were the most abundant SSR loci detected, accounting for 34.5% and 35.2%, respectively ( Table S2 ). Sixty-four sequences containing SSR loci were selected randomly for PCR amplification to identify polymorphism between the B. changningensis and its putative parents, D. farinosus and B. rigida using bulk DNA ( Table S1 ). Out of 64 primer pairs, four primer pairs, SSR14, SSR28, SSR31 and SSR34 produced polymorphic PCR amplification ( Figure 2 ). The SSR14 primers had an amplification of 139 bp in D. farinosus, whereas 144 and ~135 bp bands in B. rigida. SSR28 had an amplification size of 168 bp in D. farinosus, while B. rigida contains 159 and ~180 bp bands. The SSR31 primers had an amplification of 255 bp in D. farinosus, whereas 200 bp bands in B. rigida. Similarly, The SSR34 primers had an amplification of 160 in D. farinosus, whereas 167 bp bands in B. rigida. These results confirm that the B. changningensis had both putative parental complementary bands, except the ~135 bp bands of B. rigida.

Moreover, we analysed these four primers in 21 induvial plants of B. changningensis, the presence of putative parental bands in all individual plants confirms the heterozygous nature ( Figure 3 ). Further, we sequenced the 139, 168, 255, and 160 bp bands in D. farinosus and B. changningensis, 144, 159, 200, and 167 bp bands in B. changningensis, and B. rigida to identify the difference in the repeats of the SSR motif. The repeat units of SSR14 were CTCTC and contained a single repeat unit difference between the D. farinosus and B. rigida ( Figure S1A ). At the same time, SSR28 had three repeat unit differences in the CTC motif ( Figure S1B ). B. changningensis contains both putative parental SSR motifs. The repeat unit of SSR31 was AGG and contained a single repeat unit difference between the D. farinosus and B. rigida ( Figure S1C ). Interestingly, this SSR also has 56 bp single insertion and 1 bp deletion in D. farinosus. The repeat units of SSR34 were GCCTC and contained a single repeat unit difference between the D. farinosus and B. rigida ( Figure S1D ). These results suggest that B. changningensis might be the hybrid plant of D. farinosus and B. rigida.

The GBSSI gene is about 3 kb in length, and due to the difficulty in PCR amplification, a portion of the sequence (~820 bp) was used for analysis. Using the genomic DNA of B. changningensis, B. rigida and D. farinosus as a template for amplification, the GBSSI gene was amplified. The PCR product was cloned into TA cloning vector PMD18, and 30 clones for each species were sequenced. Sequencing results identified that a total of three different kinds of clones, BchGBSSI-1, BchGBSSI-2 and BchGBSSI-3, were present in B. changningensis. At the same time, B. rigida contains two clones, and D. farinosus has only one type of clone, BriGBSSI-1, BriGBSSI-2 and DfaGBSSI-1, respectively. Further, we aligned these sequences and identified BchGBSSI-1 and BriGBSSI-1, BchGBSSI-2 and BriGBSSI-2, BchGBSSI-3 and DfaGBSSI-1, which have similar kinds of SNPs. Moreover, a Maximum Parsimonious evolutionary tree was constructed using 23 GBSSI sequences downloaded from the NCBI database ( Figure 4 ). Similar to the above results, B. changningensis two GBSSI genes, BchGBSSI-1 and BchGBSSI-2, clustered with BriGBSSI-1 and BriGBSSI-2, respectively, and the 3^rd GBSSI gene, BchGBSSI-3 was clustered together with DfaGBSSI-1. These results confirm that B. changningensis is the hybrid between the parental bamboo species B. rigida and D. farinosus.

The complete chloroplast genome of B. changningensis was sequenced and analyzed in this study. In the current study, a total of 99,346,872 raw data was developed from the total DNA of B. changningensis. After quality trimming, the raw data was reduced to 83416080, and the GC content was 44.98%. Further, the average length of the reads was 150 bp. We mapped the quality reads to the reference cp genome B. emeiensis (HQ337797), and 14.69% of the reads were mapped. The complete chloroplast genome sequence of B. changningensis (GenBank accession: OM065947) generated was 139505 bp ( Figure 5 ). After annotation, a total of 132 genes were found in the cp genome of B. changningensis. These include 84 genes coding for essential chloroplast function, 40 tRNA genes, and eight rRNA genes. Further, the cp genome consists of a pair of identical inverted repeat regions of size 21,794 bp separated by a large single-copy (LSC) region of 83,041 bp in size and a small single-copy (SSC) region of 12,876 bp.

To identify the female parent of the B. changningensis, we aligned the cp genome sequences of B. changningensis, B. rigida (GenBank accession: MT648824), and D. farinosus (GenBank accession: OM177223). The aligned sequences were examined for the presence/absence of Indels/SNPs. Of the 132 cp genes identified in the B. changningensis, 33 contain the Indels/SNPs. Among them, 32 genes of the D. farinosus displayed Indels/SNPs, whereas the rpoc2 gene has 8 SNPs and 1 Indel in D. farinosus, 4 SNPs and 4 Indels in B. rigida and 3 SNPs in B. changningensis ( Table 2 ). Further, the D. farinosus noncoding region displayed 176 SNPs, 47 insertions, and 34 deletions. The largest insertion and deletions are 50 and 31 bp, respectively ( Table S6 ). Further, a single bp deletion was also found in B. rigida. These results indicate that the cp genomes of B. changningensis and B. rigida were almost similar, and D. farinosus contained SNPs and Indels. Thus B. rigida might be the female parent of the hybrid B. changningensis. Moreover, we also constructed a phylogenetic tree based on the extracted gene sequences from 24 complete chloroplast genomes using the Neighbor-Joining method with 1,000 bootstrap replicates ( Figure 6 ). The hybrid B. changningensis was clustered with B. rigida with high support values, confirming that B. rigida is the female parent of B. changningensis and D. farinosus is the male parent.

SSRs have been widely used in identifying hybrids as a common DNA molecular marker due to the advantages of codominance, simple operation and immunity from environmental effects (Momotaz et al., 2004; Sundaram et al., 2008). After the NGS platform’s discovery, SSR makers’ development became much easier than earlier (Subramanian et al., 2003). In the current study, we have developed the transcriptome data of B. changningensis, which can be used for phylogenetic analyses, variety or species identification and molecular plant breeding belonging to the Bambusa genus besides the hybrid identification studies. Earlier, using 15 EST-SSR markers, (Wu et al., 2009) identified four artificial bamboo hybrid species. Further, (Lu et al., 2009) also identified three artificial hybrid bamboos using SSR markers. Similarly, Yuan identified the artificial hybrids derived from Bambusa multiplex×B. chungii and B. multiplex × Dendrocalamus latiflorus using SSR markers (Yuan et al., 2021). Likewise, our study also determined that B. changningensis is the natural hybrid of D. farinosus and B. rigida.

Molecular phylogenetic studies have shown the surprising ability of gene trees to detect hybrids based on their branching behaviour (Sang and Zhong, 2000; Barraclough and Nee, 2001). ITS gene is one of the most widely used genes for phylogenetic analysis at the generic and infrageneric levels in plants (Álvarez and Wendel, 2003). However, some studies have shown that the variation of the ITS gene fragment is negligible in woody bamboo, and the taxa cannot be well distinguished (Guo et al., 2001). By contrast, GBSSI fragments have more variation than other DNA sequences previously used in woody bamboo (Guo and Li, 2004). Moreover, the GBSSI gene exists as a single copy and evolves more quickly in most Gramineae family genomes like rice, wheat, Leymus, Agropyron, and Arundinarieae (Soltis et al., 1998; Agafonov et al., 2021). As a single-copy nuclear gene, the GBSSI gene has been proven to exhibit higher genetic differentiation in introns of this gene than in ITS regions in closely related species (Mason-Gamer et al., 1998; Peralta and Spooner, 2001). Therefore, we chose GBSSI partial regions for hybrid identification and phylogenetic studies.

Our study identified three copies of GBSSI genes in B. changningensis. Among them, two B. changningensis genes clustered with two copies of B. rigida GBSSI genes, whereas another copy of B. changningensis GBSSI gene clustered with D. farinosus with strong support. These observations might be due to the hybridization of B. changningensis between the parents B. rigida and D. farinosus. Similar to our results, using the GBSSI gene, Goh et al. (2011) also identified a natural bamboo hybrid between D. pendulus and Gigantochloa scortechinii in Peninsular Malaysia. Further, Spooner identified allopolyploidy in wild potatoes using the GBSSI gene (Spooner et al., 2008).

Chloroplast DNA is unaffected by genetic recombination, which is of great significance in elucidating the complex phylogenetic relationship of plants (Huang et al., 1994). Further, cp genomes have a maternal inheritance for most angiosperms, and the molecular system reconstruction based on chloroplast genome data is very effective (Gielly and Taberlet, 1994). Therefore, most researchers used cp genes or cp whole genome sequences to identify maternal parents and analyze the genetic relationship. For example, (Nandhini et al., 2013) identified interspecific hybrid Coffea congensis × Coffea canephora and its parents based on the SNPs present in the rbcL and matK. Similarly, (Khew and Chia, 2011) identified the female parent of orchid hybrid Vanda Miss Joaquim using the matK. Further, (Uchoi et al., 2016) analyzed the phylogenetic relationship of 23 Citrus genotypes using rbcL and matK genes. Chen et al. (2019) used psbA-trnH, matK and trnL sequences as DNA barcoding to establish a rapid identification of Chrysanthemum indicum and Ch. Morifolium. Yu et al. (2022) identified Pinus sylvestris as the paternal parent for P. funebris and P. takahasii through the complete cp genome (Pinaceae species cp genome inherits paternally). Park et al. (2021) identified the maternal and paternal parents in cumber two inbred lines and their F1 hybrid based on the SNPs and indels present on the cp whole-genome sequences. In our study, to avoid interference between B. changnengesis and B. ventricosa, we extracted all the genes from 24 cp genomes and concatenate them to get the final sequences matrix to establish the phylogenetic relationship. Our study also identified that B. rigida is the maternal parent of the B. changningensis based on the SNPs and indels present on the cp whole genome sequence. Interestingly, the rpoC2 gene contains SNPs or indels in the B. changningensis and also in both parents. These results might be because we sequenced the wild species of B. changningensis and don’t know the exact individual parents. These minor variations in the SNPs and indels in the cp genome of wild species might be possible (Dally and Second, 2009). The rpoC2 gene can be considered as the fast-evolving region in these bamboo species that could be used for further evolutionary studies. In conclusion, our results identified that B. changningensis is the natural hybrid between the parents B. rigida × D. farinosus. These findings pave a path toward bamboo speciation and evolutionary studies. The SSR markers developed in this study could be practically used for genotyping of Bambusa and Dendrocalamus genus plants.