Development and Deployment of High-Throughput Retrotransposon-Based Markers Reveal Genetic Diversity and Population Structure of Asian Bamboo

Abstract

Bamboo, a non-timber grass species, known for exceptionally fast growth is a commercially viable crop. Long terminal repeat (LTR) retrotransposons, the main class I mobile genetic elements in plant genomes, are highly abundant (46%) in bamboo, contributing to genome diversity. They play significant roles in the regulation of gene expression, chromosome size and structure as well as in genome integrity. Due to their random insertion behavior, interspaces of retrotransposons can vary significantly among bamboo genotypes. Capitalizing this feature, inter-retrotransposon amplified polymorphism (IRAP) is a high-throughput marker system to study the genetic diversity of plant species. To date, there are no transposon based markers reported from the bamboo genome and particularly using IRAP markers on genetic diversity. Phyllostachys genus of Asian bamboo is the largest of the Bambusoideae subfamily, with great economic importance. We report structure-based analysis of bamboo genome for the LTR-retrotransposon superfamilies, Ty3-gypsy and Ty1-copia, which revealed a total of 98,850 retrotransposons with intact LTR sequences at both the ends. Grouped into 64,281 clusters/scaffold using CD-HIT-EST software, only 13 clusters of retroelements were found with more than 30 LTR sequences and with at least one copy having all intact protein domains such as gag and polyprotein. A total of 16 IRAP primers were synthesized, based on the high copy numbers of conserved LTR sequences. A study using these IRAP markers on genetic diversity and population structure of 58 Asian bamboo accessions belonging to the genus Phyllostachys revealed 3340 amplicons with an average of 98% polymorphism. The bamboo accessions were collected from nine different provinces of China, as well as from Italy and America. A three phased approach using hierarchical clustering, principal components and a model based population structure divided the bamboo accessions into four sub-populations, PhSP1, PhSP2, PhSP3 and PhSP4. All the three analyses produced significant sub-population wise consensus. Further, all the sub-populations revealed admixture of alleles. The analysis of molecular variance (AMOVA) among the sub-populations revealed high intra-population genetic variation (75%) than inter-population. The results suggest that Phyllostachys bamboos are not well evolutionarily diversified, although geographic speciation could have occurred at a limited level. This study highlights the usability of IRAP markers in determining the inter-species variability of Asian bamboos.

1. Introduction

Bamboo, a monocot and a major grass genera, is a group of evergreen flowering plants belonging to the subfamily Bambusoideae of the family Poaceae [1]. Although the proliferation of bamboo occurs predominantly through rhizomes, most bamboos do reproduce through seeds, flowering at least once in a lifetime. Usually, flowering intervals are long and vary between species, ranging from several to hundreds of years [2,3,4]. More than 1642 bamboo species from 75 genera are known (https://www.inbar.int), among which 100 species are commercially cultivated over 30 million hectares worldwide, particularly in Asia. Several members of the bamboo, including Asian bamboo, are recognised as fast-growing plants, growing up to a height of 35–50 m and up to 30 cm in diameter (https://www.inbar.int). Among cultivated bamboo species, the Asian bamboo can grow at a maximum rate of 100 cm a day and produces huge biomass [5].

Most Asian bamboo species are native to China, although some are known to grow in India, Vietnam and Myanmar. These bamboos account for approximately 0.8% of the forest area worldwide. Some species were introduced to Japan a hundred years ago and became naturalised. More recently, a few naturalized species from Australia, Europe and the Americas have been reported [6]. Most Asian bamboos belong to the genus Phyllostachys in the tribe Arundinarieae. They are chiefly temperate woody bamboos and are tetraploids (2n = 4x = 48) with a 2B karyotype pattern.

Bamboo wood is a non-timber natural raw material having notable industrial importance and economic value in South Asia [7]. Asia is the largest producer of bamboo products in the world, with annual international trade amounting to more than 2.5 billion US dollars (https://www.inbar.int). In spite of being an economically important perennial species, commercial bamboo remains mostly confined to natural populations. Further, the genetic diversity of bamboos has not adequately been explored. The major reasons were the difficulty in assessing the phenotypic variability of clones because of their extended growth period, perennial nature, gigantic size, propagation behavior, non-uniformity of age, long flowering cycle and the extensive area of their natural habitat. However, with the advent of molecular marker-based techniques developed in the 1980s, studies on crop genetic diversity have gained momentum. Subsequently, from 1991, a relatively limited number of molecular fingerprinting studies have been carried out to assess the genetic diversity of the Asian bamboo species using restriction fragment length polymorphism (RFLP), [8,9], randomly amplified polymorphic DNA (RAPD), [10,11], amplified fragment length polymorphism (AFLP), [12,13,14,15,16], simple sequence repeats (SSRs), [17,18], expressed sequence tags-SSR (EST-SSR), [19,20], inter-simple sequence repeats (ISSRs), [15,21] and single-nucleotide polymorphisms (SNP) [22].

In this study, we took advantage of the genome wide abundance of transposable elements to assess genetic variability in bamboo. Transposable elements (TEs) are ubiquitous genetic elements in eukaryotic genomes, capable of self-replicative transposition, affecting genome stability [23,24,25,26,27,28]. Two types of TEs have been identified based on their transposition mechanisms, namely class I retrotransposons and class II DNA transposons [28]. Retrotransposons are RNA-based TEs which duplicate themselves and move within the genome in a semi-conservative manner through a ‘copy-and-paste’ mechanism of an RNA intermediate [28,29,30]. DNA transposons, on the other hand, use a conservative style of transposition and move directly by a ‘cut-and-paste’ mechanism [31,32,33]. Retrotransposons are found in abundance, particularly in plant genomes, outnumbering DNA transposons, accounting for a significant part of the genome such as 68% in wheat [34] and 49%–78% in maize [35,36]. While exploring TEs from 44 bamboo species belonging to 38 genera, Zhou et al. [37] identified TEs as widespread, abundant and diverse in the bamboo genome. In moso bamboo, retrotransposons are reported to occupy 39%–46% [38,39,40] of the genome, accounting for about 65% of the total repetitive elements in the genome [38].

Among the two major types of retroelements, long terminal repeat (LTR)-retrotransposons and non-LTR-retrotransposons, LTR-retrotransposons constitute more than 90% of the retrotransposons found in plant genomes [41,42,43]. They have typical structural features, such as LTR sequences at both ends, transcription and reverse transcription processing signals and target site duplications [24,44]. Additionally, they possess a primer-binding site and a polypurine tract, aiding the synthesis of minus- and plus-strand DNA [45,46]. Based on their characteristics, LTR-retrotransposons are primarily divided into two superfamilies: Ty1-copia and Ty3-gypsy [42]. Recent estimates of the bamboo genome show that 63.2% of the genome is occupied by TEs [40] and 45.7% of the repeat regions belongs to Ty3-gypsy and Ty1-copia types, signifying their role in determining genome size [39]. The genome-wide analysis showed that LTR-retrotransposons are transcriptionally active in the bamboo genome and are responsible for generating 30% of small interfering RNAs (siRNAs) [47]. It has been reported that LTR-retroelements get activated by environmental stress [48,49,50]. Therefore, in the course of genome evolution, LTR-retrotransposon activity could accumulate several variations, making them an ideal source for genome-wide molecular markers [51,52,53,54].

The inter-retrotransposon amplified polymorphism (IRAP) technique produces amplified fragments that are characteristic of a dominant marker [55,56]. Amplification of the IRAP fragment between two LTR-retrotransposons is done using outward-facing primers which anneal to LTR sequences. This method needs neither restriction digestion nor ligation enzyme [57]. Due to technical ease, the IRAP method has been utilized in several studies of the genetic diversity of various plant species. Kalendar et al. [51] studied genetic diversity in barley using IRAP markers, proving their usefulness for diversity studies and the method helped to distinguish between Brazilian and Japanese rice genotypes [58]. Furthermore, it has been used in sunflower [59], Pinus [60], Lilium [61], Persian oak (Quercus brantii Lindl.) and in Bletilla striata [62], and in wild diploid wheat [63] for diversity studies, as well as for population structure and phylogenetic analyses. Guo et al. [62] reported that the results obtained using IRAP markers were similar to those obtained by start codon-targeted (SCoT) markers.

Although LTR-retrotransposons occupy a significant part of the bamboo genome, the genetic diversity information attributable to them remains mostly unknown. Furthermore, no TE-based bamboo markers have been reported to date. In the current study, we report, for the first time, the development of several IRAP markers based on the moso bamboo genome and the use of these markers to assess the IRAP-based genetic diversity and population structure of Phyllostachys bamboo.

2. Materials and Methods

2.1. Plant Material

A total of 58 Asian bamboo accessions were used in the study. The accessions included 47 distinct species belonging to the genus Phyllostachys, of which four species had 15 different varieties shared between them. There were nine varieties of Ph. edulis and two varieties each from Ph. nigra, Ph. bambusoides and Ph. sulphurea. These materials were collected from the forests of the main Asian bamboo growing regions of China spread over the provinces of Zhejiang, Anhui, Sichuang, Jiangxi, Guangdong, Hunan, Henan, Jiangsu and Taiwan. Three species, one sourced from Italy (Ph. nidularia) and two obtained from America (Ph. elegans and Ph. glauca), were also included. The details are listed in

Table 1. List of 58 Phyllostachys accessions (Asian bamboo) collected from different geographical regions used for the analysis of genetic diversity using inter-retrotransposon amplified polymorphism (IRAP) markers.

S. No.	Phyllostachys Accessions	Source Location
1.	Ph. aurea Carr.	Zhejiang, China
2.	Ph. acutiligula G. H. Lai	Anhui, China
3.	Ph. varioauriculata	Anhui, China
4.	Ph. heteroclada	Sichuang, China
5.	Ph. hispida	Zhejiang, China
6.	Ph. hirtivagina	Anhui, China
7.	Ph. edulis (Carrière) J. Houz	Jiangsu, China
8.	Ph. edulis cv. Obliquinoda	Zhejiang, China
9.	Ph. edulis cv. Pachyloen	Jiangxi, China
10.	Ph. edulis cf. huamaoyhw (Wen.) Wen	Hunan, China
11.	Ph. edulis cv. Viridisulcata	Zhejiang, China
12.	Ph. edulis (Carr.) Matsumura	Zhejiang, China
13.	Ph. edulis Mitford cv.	Zhejiang, China
14.	Ph. edulis (Carr.) Mitford cv. Gracilis	Jiangsu, China
15.	Ph. edulis cv. Tubaeformis S.Y. Wang	Hunan, China
16.	Ph. incarnata T.H. Wen, Bull. Bot. Res	Zhejiang, China
17.	Ph. nidularia	Italy
18.	Ph. sulphurea (Carr)A. et C. Riv	Anhui, China
19.	Ph. sulphurea var. viridis	Henan, China
20.	Ph. mannii	Taiwan, China
21.	Ph. virella	Zhejiang, China
22.	Ph. yunhoensxian	Zhejiang, China
23.	Ph. rubromarginata	Guangdong, China
24.	Ph. elegans McClure	America
25.	Ph. fimbriligula	Zhejiang, China
26.	Ph. viridiglaucescens	Zhejiang, China
27.	Ph. robustiramea	Zhejiang, China
28.	Ph. flexuosa	Anhui, China
29.	Ph. glauca	America
30.	Ph. zhejiangensis	Zhejiang, China
31.	Ph. angusta	Anhui, China
32.	Ph. meyeri	Henan, China
33.	Ph. vivax	Anhui, China
34.	Ph. propinqua	Anhui, China
35.	Ph. iridescens	Zhejiang, China
36.	Ph. purpureociliata	Anhui, China
37.	Ph. nuda	Jiangsu, China
38.	Ph. primotina	Zhejiang, China
39.	Ph. arcana	Jiangsu, China
40.	Ph. verrucosa	Hunan, China
41.	Ph. bissetii	Sichuang, China
42.	Ph. aureosulcata	Zhejiang, China
43.	Ph. longiciliata	Anhui, China
44.	Ph. rutila	Zhejiang, China
45.	Ph. bambusoides	Henan, China
46.	Ph. bambusoides Sieb. et Zucc	Anhui, China
47.	Ph. nigella	Zhejiang, China
48.	Ph. dulcis	Zhejiang, China
49.	Ph. aurita	Henan, China
50.	Ph. funhuanensis	Zhejiang, China
51.	Ph. platyglossa	Zhejiang, China
52.	Ph. rubicunda	Zhejiang, China
53.	Ph. atrovaginata	Zhejiang, China
54.	Ph. parvifolia	Zhejiang, China
55.	Ph. corrugata	Anhui, China
56.	Ph. nigra	Hunan, China
57.	Ph. nigra var. henonis	Jiangsu, China
58.	Ph. violascens	Zhejiang, China

. The collected plant materials were planted and maintained in red soil of a botanical garden of Fujian province. The conservation site has a subtropical monsoon climate with four distinct seasons in the year, an average rainfall from 1270 and 2030 mm a year, and an annual average temperature of 17.5 °C.

Fresh young leaves of the bamboo clones were randomly collected and surface-cleaned by gently rinsing with 70% ethanol and preserved in a polythene bag containing colour-changing silica gel (Tsingke, China). The leaf bags were stored in a deep freezer at −80 °C for further analysis.

2.2. Isolation of Genomic DNA

A modified cetyltrimethylammonium bromide (CTAB) method [64] was used to extract genomic DNA from the leaf samples. For this, the leaf samples were cut into small pieces of sizes from 3.0 to 5.0 mm. An amount of 200 mg of the cut leaf pieces was ground in a mortar using liquid N, and quickly transferred into a sterile centrifuge tube containing 850 µL preheated (65 °C) 2% CTAB extraction buffer containing 20 mM EDTA, 100 mM Tris of pH 8.0, 1.4 M NaCl, 2% CTAB, 200 mg/mL PVP and 1% β-mercaptoethanol. The tubes’ contents were mixed by gentle inversion and incubated at 65 °C for 30 min. The tubes were gently inverted twice every 10 min. After incubation, the tubes were allowed to cool for 15 min at 25 °C. After cooling, an equal volume of ice-cold phenol:chloroform:isoamyl alcohol (25:24:1 v/v) mixture was added to the tube; the contents were gently mixed and centrifuged at 12,000 rpm for 10 min. The clear supernatant was collected in another tube and an equal volume of ice-cold chloroform:isoamyl alcohol (24:1 v/v) mixture was added. The contents were gently mixed and centrifuged again at 12,000 rpm for 10 min and the supernatant was collected. To the supernatant, an equal volume of ice-cold isopropanol was added, and the tubes were kept at −80 °C for one hour to precipitate the DNA. The mix was centrifuged at 12,000 rpm for 5 min and the supernatant was discarded. To the pellet, 600 µL of 75% ice-cold ethanol was added and the tube was left standing for 10 min. The pellet-ethanol mixture was centrifuged at 12,000 rpm for 5 min. This step was repeated twice and the pellet was air-dried at 40 °C and then dissolved in 40 µL 1× TE (10 mM Tris-Cl pH 8.0 and 1 mM EDTA pH 8.0) buffer. The concentration and purity of DNA were quantified by a Nanodrop-spectrophotometer (ND1000, ThermoScientific, Wilmington, DE, USA).

2.3. Isolation of LTR-Retrotransposons and IRAP-Primer Design

In an earlier study from our lab, a total of 2,004,644 LTR-retrotransposon-related sequences were identified in the moso bamboo genome, accounting for about 40% of the moso bamboo genome [39]. The LTR sequences were identified using LTRharvest and LTR digest software [65], and the terminal repeats were analysed for similarity both at 5′and 3′ LTR regions using CD-HIT software [66]. The LTR sequences were divided into different clusters using an incremental clustering algorithm with 95% similarity criteria. LTR sequence clusters with more than 30 copy numbers were chosen as candidate sequences for IRAP primer designing. The primers were designed using Primer Premier 5.0 software (http://www.premierbiosoft.com/) and synthesized by Bioengineering (Shanghai) Co. Ltd., (Shanghai, China). Following this, IRAP primers were used for IRAP fragment amplification using appropriate polymerase chain reaction (PCR) conditions. Primers that generated a low number of amplicons were subsequently excluded from the analysis. A set of IRAP primers which provided a high proportion of alleles were, thus, finally shortlisted.

2.4. PCR Amplification of IRAP and Electrophoresis

PCR reactions were performed in 20 µL reaction mixture containing 100 ng genomic DNA, 400 nM primer, and 10 µL PCR master mix (Nanjing Nuoweizan Biotechnology Co., Ltd., Nanjing, China) and the final volume was adjusted to 20 µL by adding nuclease-free water. The annealing temperature of each IRAP primer was determined using gradient PCR. The amplification reaction was carried out in a DNA thermal cycler (DNA Engine^® Thermal Cycler—Bio-Rad). The PCR reaction was run at an initial denaturation temperature of 94 °C for 5 min, followed by 35 cycles of 30 s denaturation at 94 °C, 30 s annealing and 1 min extension at 72 °C with a final extension at 72 °C for 7 min. The annealing temperature was readjusted for each IRAP primer. The amplified product was electrophoresed in 1.5% (w/v) agarose gel at 75–80 V for 2.15 h. The separated alleles were visualised by a gel documentation system (Bio-Rad). The alleles were visually scored as 1 = present; 0 = absent using GelQuest software (https://www.sequentix.de/gelquest/).

2.5. Cloning and Sequencing of IRAP Fragments

Since IRAP markers are dominant, it is important to confirm that it is indeed LTR-retrotransposons that are selectively amplified. Alleles showing clear brightness were randomly excised from the agarose gels and were purified using a DNA gel extraction and purification kit (Simgen, Hangzhou, China). The purified fragment was ligated into the pMD18-T vector between the Not I and EcoR V restriction sites (TaKaRa, Shiga, Japan). The ligated product was transformed into Escherichia coli (E. coli) DH5α competent cells. The recombinant E. coli clones were obtained on LB agar plates containing ampicillin (100 μg/μL), X-gal (40 mg) and isopropyl β-d-1-thiogalactopyranoside (IPTG) (160 µg) and kept at 37 °C overnight. White E. coli colonies were selected by the blue-white screening method, and the insertion was verified by PCR using the corresponding IRAP primer, and the amplified fragment was sequenced (Sino Biological Inc., Beijing, China). The sequences were aligned with the original LTR-retrotransposon sequences using MEGA 7 software [67] to confirm correct amplification of the corresponding LTR sequences. Multiple sequence alignment of the PCR product sequences was also carried out with ClustalW, using the Neighbor-joining method with evolutionary distance for construction of phylogenetic tree created [68] in MEGA software with 500 bootstraps.

2.6. Marker Statistics and Genetic Relations

The total number of alleles, range of allele products, percent of polymorphism and the PIC of each IRAP marker was calculated using the binary data of the corresponding PCR amplicon using the formula PIC = 1 − [f² + (1 − f)²], where ‘f’ is the frequency of the marker in the data set. PIC for dominant markers is a maximum of 0.5 for ‘f’ = 0.5 [69]. Marker allele variation and distribution among the bamboo accessions were also worked out.

2.7. Genetic Distance and Diversity of Bamboo Accessions

The distance/similarity matrix of the bamboo accessions was constructed from allele distribution data. The distance was computed based on Jaccard’s similarity coefficient [70] and was subjected to hierarchical clustering using the unweighted pair-group method with arithmetic average (UPGMA). The UPGMA dendrogram was produced using Free Tree V. 0. 9.1. 50 software [71]. To access the dendogram reliability, over 10,000 bootstrapping (resampling) values were set. The phylogram was visualized by TreeViewX V.0.5.0 software [72].

2.8. IRAP-Statistical Fitness Analysis

To validate the dendrogram and the genetic diversity a statistical fitness analyses were performed using binary data. The cophenetic correlation coefficient (CCC) was estimated between the dendrogram and the observed dissimilarity matrix [70]. Further, 58 Asian bamboo species were classified into different groups using three-dimensional principal component analysis (PCA) based on PC1, PC2, and PC3. The analyses were performed using PAST v. 3.24 software [73], and the scatter plot of three-dimensional PCA was obtained by three-dimension PCA tool, using the OmicShare online tools (http://www.omicshare.com/tools). Both the Broken Stick model and the Jolliffe cut-off value were used to interpret the number of significant components for the total variation obtained from the PCA analysis [74,75,76]. Based on the minimum eigenvalue criteria (of more than 1), significant components were used to calculate the accuracy value with respect to the population structure and hierarchical clustering.

2.9. Analysis of Population Structure

An analysis of population structure and gene flow between the 58 Asian bamboo accessions was performed using a model-based clustering approach to divide the species into sub-populations with the help of STRUCTURE v.2.3.4 software [77]. The program uses Bayesian estimates to identify population structure, under assumptions of admixed ancestry and correlated allelic frequencies using unlinked markers [78]. No prior information was ascribed to the IRAP data while estimating sub-populations. The optimal number of sub-populations (K) was determined by running the programme with K values ranging from 1 to 10, with six independent runs for each K value. To determine the most appropriate K value, the length of the burn-in period parameter was configured to 100,000 and the number of Markov Chain Monte Carlo (MCMC) (Bayesian statistics) replications (simulations) after burn-in was set over 500,000 [79]. The optimum K value was found by an ad hoc statistic ΔK based on the percentage of variation in the log probability of the IRAP marker between successive K values using an online tool, Structure Harvester [80].

2.10. Analysis of Molecular Variance

After determining the sub-populations among the accessions tested, the analysis of the molecular variance (AMOVA) between the sub-populations was estimated using GenAlEx 6.5 software [81]. In the parameter set, a nonparametric permutation and standard permute procedure with 999 pairwise-permutations were used. These values were utilised to measure the total molecular variance between and within the populations.

3. Results

3.1. Development of IRAP Makers and Functionality Assay

A total of 98,850 LTR retrotransposons with both ends of intact LTR sequences were identified in the moso bamboo genome in the present study. The incremental clustering divided these LTR sequences into 64,281 clusters with 95% similarity criteria. Only the clusters that contained more than 30 copies of LTR sequences were considered as candidate clusters for IRAP primer development (Supplementary Table S1). Accordingly, 13 clusters with more than 30 LTR copies and at least one copy of all intact protein domains such as gag and polyprotein were shortlisted, which accounted only 0.02% of the identified clusters. These 13 clusters had a total of 696 copy numbers of LTR sequences, with an average of 53.5 copies per cluster. The highest copy number of 121 was identified in cluster number 3 followed by cluster 4, cluster 15 and cluster 22. Based on the number of clusters and LTR sequence size, a total of 90 markers were initially designed. Each marker was represented by a primer that acted both forward and reverse primers (Figure 1). Among these, 26 primers (29%) that showed proper amplification with clear allele patterns and high reproducibility were shortlisted. Finally, only 16 primers (18%) those showed clearly distinguishable polymorphism, and were, therefore, chosen for further analysis.

A functionality check of the amplicons generated by the selected 16 IRAP markers revealed consistently well-resolved and reproducible amplicon patterns among all the 58 Asian bamboo accessions (Figure 2). There was a total of 215 scorable amplicons (alleles) produced of which 214 were polymorphic (99.5%). Polymorphic alleles were produced with an average of 13.3 alleles per marker. Across the test population, a total of 3282 polymorphic amplicons were generated with an average of 56.6 amplicons per accession. The allele numbers produced per marker ranged between 8 (CL54-R) and 16 (CL34-R and CL63-R) with a size variation ranging from 200 bp to 2700 bp (

Table 2. List of inter-retrotransposon amplified polymorphism (IRAP) primers with description of amplicons used for the analysis of the genetic diversity and population structure of different Phyllostachys species (Asian bamboo).

IRAP Marker	Primer Sequence (5′-3′)	IRAP Primer Location in LTR Retrotransposons (bp)	Allele Size Range (kb)	No. of Polymorphic Alleles	Total Number of Amplicons in the Population	Average Allele Frequency	PIC Value	Melting Temperature (Tm) °C
CL3-F	TATAAAGGTAGCTTTCGGGTATG	428338_436487	0.65–2.7	13	231	0.306	0.350	52
CL4-R	CTGGTATATAGCTGTTGAGCGACG	82500_84292	0.6–2.4	13	192	0.255	0.351	57
CL15-R	CTCGTGTATTCTCCCTTTGC	221059_224205	0.43–2.5	12	227	0.326	0.400	53
CL22-F	TGATCAGAGAAGAAAGGGGA	49377_57704	0.42–2.7	13	125	0.166	0.256	52
CL22-R	CACGCAGAGAGATTGACACG	49377_57704	0.35–2.7	15	173	0.198	0.319	56
CL34-F	GAACGATTACCTCACAGACA	4137_7249	0.45–2.7	15	183	0.210	0.262	52
CL34-R	GAGCAATAAAGAGAAGCCCG	4137_7249	0.4–2.2	16	273	0.294	0.334	53
CL37-F	AGATTGTTTGATTCGGGGGG	153045_162321	0.3–2.5	12	183	0.263	0.325	54
CL37-R	AGCGGCGTGGAGGAGTTACC	153045_162321	0.5–2.7	13	234	0.310	0.376	61
CL42-R	GCAACAACAAACCCTAAAAA	18144_27070	0.3–2.7	14	230	0.283	0.369	50
CL54-R	GCAAGAACATAAGAACAGAA	56871_64575	0.3–2.0	8	209	0.450	0.336	55
CL58-F	AAGGAATCGTCAGTCAACAA	28670_35571	0.3–2.5	11	202	0.317	0.369	52
CL59-F	TGTCAGACAGTACAGGTGCT	1169_4536	0.2–2.7	14	204	0.251	0.329	55
CL61-F	ATGACATAGGGCACACCAGA	680968_684115	0.37–2.0	14	260	0.320	0.275	55
CL62-F	TAAATAGGGAACGAGGAGCC	178936_183841	0.5–2.7	15	190	0.218	0.306	56
CL63-R	ACATTGTTTGATTCGGGGGG	248490_260632	0.45–2.7	16	166	0.179	0.273	55
Mean	-	-	-	13.37	208.75	0.272	0.327	-

Note: LTR, long terminal repeat; PIC, polymorphic information content.

).

The marker CL59-F produced the widest range of amplicon sizes, from 200 bp to 2700 bp, while the amplicons generated by the marker CL61-F had the shortest size range, 370–2000 bp. Furthermore, the marker CL34-R produced the highest number of total amplicons (273) in the population followed by CL61-F (260), CL37-R (234), CL3-F (231) and CL42-R (230). All the primers produced 100% polymorphic alleles, except the marker CL22-F, which produced one monomorphic allele and 13 polymorphic alleles showing a polymorphism of 92.8%. The average allele frequency of the IRAP primers ranged between 0.166 (CL22-F) and 0.450 (CL54-R). The polymorphic information content (PIC) values ranged from 0.256 (CL22-F) to 0.400 (CL15-R), having an average of 0.327.

3.2. IRAP Amplicon Variation within and between the Accessions

Of the total of 3340 amplicons obtained, 3282 (98.3%) were polymorphic, and 58 (1.7%) were monomorphic. Among the 47 Phyllostachys species used in this study (Table 1), Ph. edulis generated the highest number of IRAP amplicons, having an average of 60 alleles across its nine varieties. Ph. edulis generated 72 alleles, followed by Ph. edulis cv. viridisulcata with 71 alleles. Nine moso bamboo varieties generated a total of 538 amplicons, of which eight alleles were monomorphic in all the nine moso bamboo varieties. Ninety-four alleles showed no amplification among the Ph. edulis varieties.

3.3. Multiple Sequence Alignment of IRAP-PCR Products

The multiple sequence alignment of IRAP-PCR amplicons using designed IRAP primers proved that they were indeed LTR-retrotransposon sequences. The neighbor-joining phylogenetic tree showed that the IRAP product sequences were genetically different (Figure 3), suggesting that all sixteen IRAP primers amplified unique PCR bands. A total of 14 bootstrap values were obtained, with 9 values ranging from 56% to 100%, which further confirms that the alleles were unique. The above results suggest that LTR-retrotransposon have unique amplicons for genome integrity and sizes.

3.4. Diversity Analysis of Asian Bamboo Accessions

The genetic similarity values of the 58 Asian bamboo accessions are shown in Supplementary Table S2. Jaccard similarity coefficient values ranged from 0.09 to 0.98 with an average similarity value of 0.25 among the 58 accessions. A total of 1653 pair-wise similarity coefficients were obtained, among which only seven pairs showed high similarity that ranged between 0.80 and 0.99. Twenty-eight pair wise similarity coefficients were intermediate, with values ranging from 0.60 to 0.79. A large proportion (1618 or 98% of similarity coefficients) were less than 0.59 (Supplementary Table S2). The lowest similarity coefficient of 0.091 was seen between Ph. longiciliata and Ph. edulis (Carr.) Matsumura followed by that between Ph. arcana and Ph. edulis (Carr.) Matsumura (0.092). The highest similarity of 0.984 was observed between Ph. nigra var. henonis and Ph. edulis cv. Pachyloen, followed by the similarity between Ph. edulis (Carr.) J. Houz and Ph. edulis (Carr.) Mitford cv. Gracilis.

Based on the UPGMA clustering, the bamboo accessions were grouped into five clusters (Figure 4) viz. Phyllostachys clusters, PhC1 to PhC5. The cluster, PhC1 had 19 accessions grouped within, with an average similarity of 0.33. The second cluster, PhC2 which had a mean similarity of 0.34 contained 13 accessions. The cluster PhC3, which had only two accessions had a similarity of 0.26. The highest average similarity of 0.42 was observed in the sub-cluster PhC4, which included 17 accessions. All of the highly similar bamboo accessions were found in this cluster. The remaining sub-cluster, PhC5 encompassed seven accessions, and had average similarity value of 0.33. The clusters themselves showed affinity among themselves, with the first three clusters getting grouped together, as well as the remaining two getting grouped into a separate group. These two groups showed distinct separation with a bootstrap confidence of 100%. In the dendrogram constructed, bootstrap values ranged from 61% to 100% between clusters and the average bootstrap value observed in this study was 81%. For the clusters, PhC1 separated from PhC2 and PhC3 with bootstrap value of 74%, while PhC2 and PhC3 were separate for 91 times out of 100 bootstrap iterations. For the second group, the clusters PhC4 and PhC5 were distinct for 85% of the bootstrap resamplings.

The lowest similarity values observed between several of the accessions studied were due to the highest percentage (98%) of polymorphic bands generated by the IRAP primers. Of the nine varieties of moso bamboo, Ph. edulis, eight were found to be included in PhC4, while the remaining, Ph. edulis cv. tubaeformis, was found to be placed in sub-cluster PhC5. Most Ph. edulis variants were from the Zhejiang, Jiangsu and Jiangxi provinces of China (Table 1). Among 25 accessions collected from the Zhejiang province of China were found placed in different clusters. No specific clusters represented the species collected from the Zhejiang and Anhui province. The species collected from the Jiangsu province of China were grouped in clusters PhC1 and PhC4, and those from Henan provinces in PhC1 and PhC2. Moreover, two varieties of Ph. nigra were found to be distributed between PhC2 and PhC4. Similar was the case with two varieties belonging to the species, Ph. bambusoides. However, varieties of Ph. sulphurea were found to be included in PhC1 and PhC4. The bootstrap value for Ph. nigra var. henonis and Ph. edulis cv. Pachyloen as well as for Ph. edulis cv. Tubaeformis S.Y.Wang and Ph. sulphurea (Carr)A. et C. Riv was 100%. A bootstrap confidence of 98% was seen for Ph. edulis (Carrière) J. Houz and Ph. edulis (Carr.) Mitford cv. Gracilis. The bootstrap values clearly displayed that a good majority (79%) of cluster nodes were well fitted. None of the nodes were found with very low bootstrap value. Based on the resampling method, it was confirmed that the IRAP markers markedly distinguished the Asian bamboo accessions.

3.5. Statistical Fitness Analysis of Clustering Pattern

Two different statistical analyses such as cophenetic correlation coefficient (CCC) and principal component analysis (PCA) were carried out to confirm the grouping pattern of the Asian bamboo species. The CCC value of 0.8848 confirmed that hierarchical clustering was in significant agreement with the similarity matrix obtained from Jaccard similarity coefficients. Based on the IRAP markers the 58 Asian bamboo accessions indicated significant genetic diversity at those loci.

Three-dimensional PCA of the 58 Asian bamboo accessions, based on the variance-covariance matrix, displayed 13%, 7.9%, and 6.5% of the total variance for the first, second, and third component axes, respectively (Figure 5). A total of 57 axes (principal components) were extracted of which five components were retained based on the broken-stick model (Supplementary Figure S1). The selected components accounted for 38.1% of total variation. Three-dimensional plotting of bamboo accessions based on first three components showed four groups (Figure 5). Further resolution of the genotype grouping using hierarchical clustering of Euclidean distances based on the component scores for the five significant PCs, indicated the members under each group. The first group (Group 1) showed grouping of 13 accessions, Group 2 contained 19 accessions, Group 3 had 18 accessions and Group 4 with eight accessions (Figure 6). Comparing the clusters based on the Jaccard’s similarity coefficients, the Group 1 had members drawn from PhC4 and PhC5, while Group 4 was almost entirely was drawn from PhC4, and contained seven out of nine Ph. edulis accessions used in the study. The Group 2 was exactly similar to PhC1 in both number and membership of accessions. The remaining group, Group 3 contained all the accessions that were members of PhC2 and PhC3, together with three accessions drawn from PhC4 and PhC5. The grouping pattern of the accessions revealed that based on the IRAP polymorphism, Group 2 was more robust and isolated from the subsequent groups of Group 3, Group 1 and Group 4. Group 2 had an accuracy value of 100%, followed by Group 3 (83%), Group 1 (54%) and Group 4 (47%) from the dendrogram based clusters.

3.6. Population Structure of Asian Bamboo Accessions

The population structure analysis performed using IRAP markers to understand the genetic relationship among Asian bamboo accessions, revealed the existence of four apparent sub-populations in the test accessions. The Structure Harvester picked the maximum ΔK value when the inferred number of sub-populations (K) was at four (K = 4), having the maximum delta K value of 265.41 (Figure 7) with the lowest standard deviation (2.34) of the parameter, LnP(K). The resolved sub-populations were designated as PhSP1, PhSP2, PhSP3 and PhSP4 (Figure 8). The first sub-population, PhSP1 accommodated 13 accessions, whereas the remaining sub-populations PhSP2, PhSP3 and PhSP4 carried 19, 17 and 8 accessions, respectively. One of the accessions, Ph. flexuosa remained out of the sub-populations being an admixture of all the four sub-populations. Respective proportions of memberships were 23.0% for PhSP1, 35.7% for PhSP2, 27.1% for PhSP3 and 14.2% for PhSP4. The allele frequency divergence between the sub-populations was maximum between PhSP3 and PhSP4 (0.20), followed by between PhSP2 and PhSP4 (0.19). The lowest divergence was observed between PhSP1, PhSP2 and PhSP3. The estimated values of expected heterozygosity that can be construed as the average distance between members within each sub-population were 0.30, 0.27, 0.29 and 0.08 for PhSP1 to PhSP4 respectively. The proportion of the total genetic variance (F_ST) explained by the sub-populations ranged from 0.21 (PhSP1) to 0.80 (PhSP4), with the remaining PhSP2 and PhSP3 having F_ST values of 0.24 and 0.26 respectively. The inferred ancestry coefficients (Q values) for individual accessions carried by the PhSP1 ranged from 0.63 to 0.99, while that of PhSP2 was between 0.73 and 0.99 (Figure 8). Similarly, PhSP3 had a range of 0.50 to 0.99 for the Q values, and the range for PhSP4 was between 0.74 and 1.00. The average Q values for each sub-population were 0.90 for PhSP1, 0.94 for PhSP2, 0.86 for PhSP3 and 0.91c for PhSP4.

Individual members of the sub-populations included Phyllostachys species such as Ph. hispida, Ph. hirtivagina, Ph. virella, Ph. varioauriculata, Ph. heteroclada, and Ph. mannii having placed in PhSP1, that possessed maximum allele frequency of that sub-population. Similarly, fourteen accessions were found to possess high inferred ancestry coefficients for PhSP2, eight accessions for PhSP3 and five for PhSP4. Seven of the moso bamboo varieties (Ph. edulis) were placed in PhSP4, and four of them had the maximum inferred membership coefficient together with one accession, Ph. bambusoides (Figure 8B). All the four sub-populations contained accessions with admixtures of alleles from all, with PhSP3 carrying maximum of admixed accessions, followed by PhSP1 and PhSP4.

We have drawn a consensus of distribution of accessions across the sub-populations by comparing membership of respective groups obtained from hierarchical clustering using Jaccard coefficients and PCA (Supplementary Table S3). It was identified that PCA showed a very close accuracy with respect to the population structure than that obtained by hierarchical clustering. Of the total of 58 accessions, only 46 accessions were commonly shared between the distinct groups obtained among population structure, principal component, and hierarchical clustering analyses (Supplementary Figure S2). Whereas, both population structure and principal component analyses commonly shared a total of 54 accessions from four sub-clusters, with an accuracy value of 93%. The sub-population, PhSP1 had 11 consensual accessions, with an average similarity of 0.28. The remaining sub-populations, PhSP2, PhSP3, and PhSP4 shared 19, 16, and 8 accessions, respectively sharing average similarity values of 0.33, 0.29 and 0.65 (

Table 3. Description of consensual Asian bamboo accessions picked out from the groupings obtained from dendrogram and principal component analysis (PCA) with respect to deduced population structure based on inter-retrotransposon amplified polymorphism (IRAP). There were a total of 54 accessions falling under different sub-populations, PhSP1 to PhSP4.

Sub-Population Name	Phyllostachys Species Name	Membership No.	Average Similarity
PhSP1	Ph. acutiligula G. H. Lai; Ph. heteroclada; Ph. hirtivagina; Ph. varioauriculata; Ph. hispida; Ph. rubromarginata; Ph. aurea Carr.; Ph. sulphurea (Carr)A. et C. Riv; Ph. virella; Ph. mannii; Ph. yunhoensxian.	11	0.28
PhSP2	Ph. bissetii; Ph. aureosulcata; Ph. propinqua; Ph. verrucosa; Ph. vivax; Ph. nuda; Ph. arcana; Ph. meyeri; Ph. primotina; Ph. glauca; Ph. iridescens; Ph. angusta; Ph. sulphurea var. viridis; Ph. rutila; Ph. purpureociliata; Ph. zhejiangensis; Ph. fimbriligula; Ph. elegans McClure; Ph. viridiglaucescens.	19	0.33
PhSP3	Ph. atrovaginata; Ph. parvifolia; Ph. rubicunda; Ph. funhuanensis; Ph. longiciliata; Ph. aurita; Ph. platyglossa; Ph. corrugata; Ph. dulcis; Ph. violascens; Ph. bambusoides; Ph. nigella; Ph. edulis (Carr.) Mitford cv. Gracilis; Ph. edulis cv. Tubaeformis S.Y. Wang; Ph. nigra; Ph. nigra var. henonis.	16	0.29
PhSP4	Ph. edulis cf. huamaoyhw (Wen.) Wen; Ph. edulis cv. Pachyloen; Ph. bambusoides Sieb. et Zucc; Ph. edulis (Carrière) J. Houz; Ph. edulis cv. Viridisulcata; Ph. edulis Mitford cv.; Ph. edulis cv. Obliquinoda; Ph. edulis (Carr.) Matsumura.	8	0.65

).

3.7. Analysis of Molecular Variance

The genetic diversity of the four sub-populations of the Asian bamboo species divulged from the genome wide IRAP polymorphisms indicated significant level of genetic variation within and among populations. The sums of squared deviation (SS) values within populations and among populations were 1518.27 and 492.12, respectively with the variance (MS) values of 28.12 and 164.04. The percentage of total molecular variance within populations was 75%, and, among populations, it was 25%. A molecular variance of 75% indicates that a strong genetic differentiation occurred within populations. These values confirm the size of the four sub-populations and the ratio of the admixture of alleles in the four sub-populations. The genetic diversity value (P) of the Asian bamboo species was highly significant (p < 0.001) at two hierarchical levels (among populations and within populations).

4. Discussions

To support genetic bamboo improvement, it is essential to understand the genome-wide variation and diversity among the breeding material. Since TEs form 63% of the bamboo genome [40], of which more than 46% are LTR-retrotransposons, retrotransposon-based DNA-fingerprinting could be an ideal technique to study the genome wide diversity of closely related species or breeding lines [82,83]. The use of retrotransposon-based molecular markers to study plant diversity is a cheap and rapid technique which can provide potentially useful molecular information augmenting other DNA- based markers [82]. The availability of the complete bamboo genome (BambooGDB, http://www.bamboogdb.org) [84] enabled a detailed exploration of the LTR-retrotransposons and their copy numbers based on their structure and subsequently for the development of IRAP primers. Since the development of IRAP primers is a one-time investment, potential primers can be continually be applied [85,86,87] for studying genetic diversity of bamboo species and their corresponding congeners.

In the current study, we confirmed that LTRs of bamboo species had different copy numbers, and their sequences contained full complement of LTR retrotransposons. The structure analysis revealed that they are transcriptionally active, and could be functional. The IRAP primers developed from the high copy number clusters although could amplify multiple-sites of the genome, their number seems to be low when considering the potential diversity of the Phyllostachys accessions used. The designed IRAP primers produced both polymorphic and monomorphic alleles that enabled the use of IRAP DNA fingerprinting to assess the diversity among Asian bamboo accessions. Based on this genome-wide analysis, we concluded that only 29% of IRAP markers showed polymorphism implying that inter-LTR regions in the studied genomes of bamboo species were significantly conserved. This implies that the bamboo genome is still under evolution and LTRs are not very active in contributing to the genome wide variations. This is, to the best of our knowledge, the first report of IRAP primer development and the first report of genetic diversity in bamboo using IRAP based fingerprinting.

4.1. Genetic Diversity of Asian Bamboo Species

The IRAP primers designed and used in the current study offered valuable genetic information about economically important Asian bamboo species and their evolutionary pattern. From the pattern of IRAP polymorphism, we conclude that in spite of their abundance in the bamboo genome, LTR retrotransposons significantly lack polymorphism among the 58 accessions investigated in this study. Although, these accessions belonged to different Phyllostachys species as recognized from their collection sites, the lack of IRAP polymorphism suggested significant level of conservation within the bamboo genome. Possible reason for this genome integrity could be the low frequency of sexual reproduction in bamboos that is monocarpic and occur at long breeding cycles [88]. Moreover, bamboos propagate predominantly clonally, leading to less variation on the genomes. However, over a longer period of time, geographic isolation can conserve localized speciation in the bamboo clones leading to identification of different bamboo species, with some level of morphologic variations. Possibly, active LTR-retrotransposons could be playing a major role in this type of speciation, as they are capable of bringing in spontaneous changes in the genome at random loci. We have observed a few polymorphic loci in the present study, but they showed high level of polymorphism between different sub-populations. This could be attributed to the changes infused by active retrotransposons across the genomes. In an earlier study, analysis of 78 accessions of the Asian bamboo species using 23 microsatellite markers had revealed an average of 2.78 alleles per primer [89]. In contrast, the IRAP markers could produce as many as 13 alleles per marker locus, indicating their effectiveness in revealing the genetic diversity of bamboo genomes in the current study. It may be emphasized that, 98% of average polymorphism exhibited by the IRAP amplicons is the highest reported for bamboo. IRAP occurs due to random insertion of retroelements on the genome, resulting in a length variation of the interspersed regions flanked by two elements. If such insertions affect functional genes, the resultant variations are further clonally propagated in the population, resulting in a transient speciation until sexual reproduction occurs. Such transient speciation can be either transmitted to further sexual generations or can be lost or can subsequently crate novel variations. Therefore, species level genetic diversity of the Asian bamboo accessions can be considered predominantly intermediary and is in the course of evolution that may last for several thousands of years to come. High level of genome wide polymorphism using IRAP markers was also reported in other crop species such as among octoploid triticale plants, where 85% polymorphism was observed [90]. Similarly, 94%, 79% and 74% polymorphism have been reported in Bletilla striata [62], Lallemantia iberica [91] and Schistosoma japonicum [92], respectively, using IRAP markers. Like earlier instances, we could employ IRAP markers reliably in assessing the genetic diversity in bamboo.

Based on the inter-retrotransposon distances, 58 Asian bamboo accessions could be separated into four sub-populations. Considering the fact that the study included 47 species of Phyllostachys, the four sub-populations could be considered very low if the bamboo populations are in Hardy-Weinberg equilibrium with typical random mating behavior. This further consolidated our inferences that the study accessions used in the study were transient species adapted to local niches and propagated clonally. There was a significant lack of molecular allelic pattern representing different species, except for variations at few loci across the genome that could have been accumulated more recently in the evolutionary time scale due to transposition activity. These observations were similar to previous report on 78 Asian bamboo accessions using SSR markers, which grouped them into three classes [89]. Similarly, two major clusters have been reported in 50 varieties of Bletilla striata using IRAP markers [62]. Further, Zhao et al. [89] have reported that the genetic variation between Ph. nuda and Ph. propinqua was 0.2143 using SSR markers. In the current study, we obtained similar variation (0.371) between Ph. nuda and Ph. propinqua. Both these species were characterized with having or not having bristles on the back of sheath. Ph. vivax and Ph. aureosulcata have purple-green or yellow-green with purple culm sheaths and clustered in PhSP2, similar to Zhao et al. [89] study. The current results are consistent with current bamboo’s taxonomic classification and agreed with the morphological classification [93]. Our study further proved that, model based and PCA based approaches were significantly better for resolving the population structure of bamboo, in the event of having a few polymorphic IRAP markers, that has produced significantly good number for highly polymorphic alleles. The CCC of 0.88 obtained from the dendrogram based on Jaccard coefficient is suggestive of this. In a previous report on 200 tree accessions in the 20 groups of Olea europaea using IRAP markers, CCC value was 0.96, indicating a good fit between the similarity matrix and the dendrogram [94].

4.2. Population Structure of Asian Bamboo

Structure analysis used in our study implements a model-based approach for inferring population structure using unlinked genotype marker data, identifying genetically distinct populations, admixtures of alleles in populations and assign individuals to specific sub-populations [77]. Using this approach, Jiang et al. [95] identified three sub-populations in 803 accessions of moso bamboo using 20 SSRs markers. Further, Nachimuthu et al. [96] have identified two sub-populations in 192 accessions of rice using 61 SSR markers. In the current study too, the initial analyses of genetic diversity using the hierarchical clustering as well as the PCA were suggestive of a low-level population differentiation within the study panel of Asian bamboo accessions. Therefore, we have fixed a sub-population range of 1 to 10 for the model, assuming admixed populations and correlated allele frequencies. Analysis revealed an optimum population structure consisting of four sub-populations.

The inferred ancestry coefficients of Asian bamboo accessions provided the genetic relationship and gene flow pattern between the sub-populations. Few members of all sub-populations had admixture of alleles, while some members were specifically grouped with maximum frequency of sub-population specific alleles. From the molecular diversity pattern of the Asian bamboo accessions, we concluded that the PhSP2 was the sub-population with particularly lower number of admixed clones. Conversely, the PhSP3 had significantly high level of admixture especially from PhSP2. Similarly, the admixture of alleles was found between PhSP1 and PhSP4. However, we could observe that genetic variations among populations was significantly lower than that within populations indicating several subtle allelic variations among the members of each sub-populations. This was similar to previous studies for bamboo species such as Melocanna baccifera and Bashania fangiana (Dwarf bamboo) using ISSR and AFLP markers, respectively [97,98].

Among the sub-populations, lowest within population variation was observed with PhSP4. This sub-population contained most of the Moso bamboo accessions of the species, Ph. edulis. This genetically very close group, however, contained one accession of Ph. bambusoides, which could be suspected as a mistaken nomenclature. The significant similarity within this group, could be well explained because all the members were different commercial varieties of Moso bamboo. This category of bamboo is the most commercially exploited bamboos and are chiefly propagated clonally. Moso bamboos also offer long flowering intervals, and long breeding cycle [99]. Further, Moso bamboo displayed very low admixture of alleles, which suggests that LTR-elements had a little role in defining its genetic structure.

5. Conclusions

The aim of the current study was to explore the genome wide abundance of LTR-retrotransposons in the Asian bamboo accessions and to development IRAP based markers for investigating genetic diversity and population structure. Although the variation among the interspaces between the LTR retrotransposons in bamboo species was low, few loci showed apparently high polymorphism aiding the analyses. Since transposon activity is related to environmental factors, geographic speciation could be one of the reasons for high IRAP based diversity at certain loci. This is the first report of population structure using IRAP markers in the Asian bamboo species. From the observed pattern of genetic diversity, it is reasonable to assume that the ancestors of Asian bamboo could be few in number with limited variability, which on evolution, adaptably speciated into different species, with subtle genetic change compared to other rapidly multiplying cross-pollinated species. Each of the IRAP primer had unique differentiation, and this marker system offered highly efficient and reproducible alleles for studying inter-retrotransposon-based genetic diversity.