Short communication
Development of microsatellite markers for the mycoheterotrophic species Burmannia nepalensis (Miers) Hook.f. based on RAD sequencing
2021 Volume 96 Issue 6 Pages 293-298
Details
2021 Volume 96 Issue 6 Pages 293-298
Mycoheterotrophic plants can derive carbon from fungi rather than from photosynthesis. Habitat destruction and sensitivity to environmental perturbation may result in the loss of biodiversity including genetic variation of mycoheterotrophic plants. Burmannia nepalensis (Miers) Hook.f. (Burmanniaceae) is a mycoheterotrophic plant with a wide distribution across southern China and southern and eastern Asia. As part of our endeavor to reveal population genetic patterns of mycoheterotrophic plants, fifteen microsatellite loci were developed by RAD (restriction site-associated DNA) sequencing in 89 individuals from four populations of B. nepalensis. A total of 49 alleles were amplified. The number of alleles per locus ranged from two to six with an average of 3.3. The observed and expected heterozygosity per population varied from 0.000 to 1.000 and from 0.000 to 0.722, respectively. A transferability test showed that only one to five loci could be cross-amplified successfully in four other congeneric species of Burmannia. These markers can be used to reveal population genetic diversity in B. nepalensis, and will help to elucidate the evolutionary history and to enhance efforts for conservation of mycoheterotrophic plants.
Mycoheterotrophs represent a fascinating group of plants. They derive carbon from mycorrhizal fungi, which often simultaneously have a mutualistic interaction with surrounding photosynthetic plants and are thus often referred to as cheaters or exploiters (Bidartondo, 2005; Merckx et al., 2009). The structures of the flower, stem and leaf in mycoheterotrophic plants always show varying degrees of simplification or reduction (Leake, 1994), which is of great importance in evolution and has been attracting much attention from evolutionary scientists for many years (Leake, 1994; Merckx, 2013). However, their tiny size, rarity and strong seasonality make them hard to notice (Merckx et al., 2006), with the result that very limited research has so far been conducted on mycoheterotrophs. In addition, the biodiversity of mycoheterotrophic plants has been experiencing threats from human activities including habitat destruction and fragmentation, and climate change (Merckx et al., 2013). To evaluate the severe impacts of these contributing factors on mycoheterotrophic species, it is important to examine their population genetic diversity, structure and demographic history. However, there exists a major gap in understanding the genetic basis of mycoheterotrophic species: studies on their population genetics have been scarce (Bidartondo, 2005; Klooster et al., 2009), comprising one report for Burmannia wallichii (Zhang and Saunders, 2000) and several for species in the Ericaceae (Klooster and Culley, 2010; Shutoh et al., 2017) and Orchidaceae (Chen et al., 2011; Tsai et al., 2014).
Burmannia L. (Burmanniaceae) is a major genus of mycoheterotrophic plants that contains about 60 species distributed in tropical, subtropical and warm temperate regions. It includes autotrophic, partial-mycoheterotrophic and mycoheterotrophic species with a host preference and specialization toward fungal assemblages (Zhang, 1999; Wu et al., 2010; Zhao et al., 2021), which makes it a good resource to investigate evolutionary mechanisms of mycoheterotrophs. Burmannia nepalensis (Miers) Hook.f., a small (~10 cm), white and slender annual herb with numerous tiny seeds, is one of the few mycoheterotrophic species in the genus with a wide distribution across Yunnan, Guangxi, Guangdong and Hunan provinces of China, as well as in regions from South to Southeast Asia. This species occurs in dense and humid forests with piles of fallen leaves (Wu et al., 2010) and often grows with a clustered pattern at a local site (Zeng, 2021). The only published studies of B. nepalensis are two new observation reports (Maxwell, 1998; Nuraliev et al., 2018) and an examination of its mycorrhizal specificity (Ogura-Tsujita et al., 2013); there are no reports of pollinators. Therefore, nothing is known about the reproductive ecology and population genetics of B. nepalensis. Such research, however, is of pressing importance given the ever-growing rarity of suitable habitats resulting from anthropogenic activity and global climate change (Merckx et al., 2013). In the present study, we aimed to develop microsatellite markers for B. nepalensis, which will enable future investigations of genetic structure and phylogeographic relationships of populations across environmental gradients, and inference of the demographic history of the species.
Of the various PCR-based molecular markers, microsatellites or simple sequence repeats (SSRs) have been the most popular choice for use in studies of population genetics due to high rates of polymorphism, ease of application and well-developed statistical methods (Fox et al., 2019). Although next-generation sequencing-based technology has provided more powerful methods, such as restriction site-associated DNA sequencing (RAD-seq), for addressing questions in ecology, conservation biology and evolution, microsatellite markers remain an efficient and relatively cheap approach in non-model species, especially those with low-abundance and poor-quality DNA templates. In this study, we characterized fifteen microsatellite markers for B. nepalensis based on RAD-seq. This is the first report of microsatellite markers in Burmannia. These markers will provide opportunities for future genetic studies on B. nepalensis as well as other congeneric species.
Individuals were collected from four populations (97 individuals in total) in China: Wenshan (XSJ), Yunnan Province; Baise (YC), Guangxi Province; Wuyanling (WYL), Zhejiang Province; and Conghua (CH), Guangdong Province. In each population, 16 to 29 individuals growing at least 1 m from each other were picked and kept in silica gel. Because of the tiny size of the plants, their genomic DNA was extracted using a Plant Genomic DNA kit (TIANGEN) from the whole dry plant. One individual from the CH population was selected for RAD-seq, which was performed by Novogene following the protocols in Baird et al. (2008). Briefly, genomic DNA was digested with EcoRΙ restriction enzyme prior to the ligation of P1 adapters containing an Illumina sequencing primer site, a forward amplification primer site and a nucleotide barcode. The barcoded samples were then pooled and randomly sheared to ~500 bp using a Fisher Scientific sonicator (700 W). After electrophoresis in 1% agarose, fragments of 300–700 bp were isolated using a MinElute Gel Extraction Kit (Qiagen). After polishing the ends of the DNA using the Quick Blunting Kit (NEB), the size-selected DNA was then ligated to P2 adapters. Fragments with both adapters were amplified to generate the RAD library. The prepared library was sequenced on an Illumina HiSeq 2000 with paired-end read length of 150 bp. After quality filtering, contigs longer than 150 bp were assembled by Velvet software (Zerbino et al., 2009). Sequences containing SSRs were identified using an SSR search script developed by Novogene with the following criteria: at least six repeats for dinucleotides to hexanucleotides, a minimum product size of 100 bp and a minimum distance between two different SSRs of 12 bp. Subsequently, SSR primer pairs were designed by Primer3 (Untergasser et al., 2012).
Two hundred markers with at least six consecutive repeat units were chosen randomly for screening, including 140 markers for dinucleotide repeat units and 60 for trinucleotide repeats. Eight individuals from the four populations (two individuals per population) of B. nepalensis were used to test whether these markers could be amplified successfully. PCR was carried out in 10 μl, including 1 μl of genomic DNA, 0.06 μM M13-tailed (5′-TGTAAAACGACGGCCAGT-3′; Schuelke, 2000) forward primer, 0.25 μM reverse primer, 1 μl of fluorescent dye-labeled M13 (FAM, TAMRA, HEX or ROX; TSINGKE) primer, 2 μl of ddH2O and 5 μl of 2×T5 Super PCR Mix (TSINGKE). Touchdown PCRs were performed under the following thermocycling conditions: an initial denaturation of 4 min at 94 ℃; 7 cycles of 94 ℃ for 40 s, 60 ℃ for 40 s with a decrease of 1 ℃ per cycle, 72 ℃ for 80 s; followed by 28 cycles of 94 ℃ for 40 s, 53 ℃ for 40 s, 72 ℃ for 80 s; and a final extension at 72 ℃ for 8 min. PCR products were scanned with an ABI 3730xl DNA Analyzer (TSINGKE) using the internal size standard GeneScan 500 LIZ (Applied Biosystems) and genotyped by GeneMarker V2.2.0 software (SoftGenetics). Finally, microsatellite loci that could be amplified successfully with clear peaks and suitable fragment lengths in all eight individuals were further tested for polymorphism in the other 89 individuals from the above four populations (XSJ, YC, WYL and CH). The PCR systems and thermocycling conditions followed the protocols described above. Scanning of products as well as allele calling and binning were implemented using the procedure described above. Transferability was also tested in four species of Burmannia (B. oblonga, B. coelestis, B. wallichii and B. decurrens) with 16 to 19 individuals per species, using the same DNA extraction and amplification methods as described above.
During the field observations, we did not find any physical link between rhizomes of different individuals, indicating the absence of vegetative propagation in B. nepalensis. However, we cannot exclude the possibility of clonal reproduction, because a previous study has reported the phenomenon of apomixis in Burmannia species (Zhang, 1999). Therefore, individuals with the same multilocus genotypes (MLGs) were first checked in Genclone 2.0 (Arnaud-Haond and Belkhir, 2007). One individual per MLG was selected to generate a subset of 35 individuals for the following statistical analyses. The number of alleles per locus (A), observed heterozygosity (HO) and expected heterozygosity (HE) as well as deviations from Hardy–Weinberg equilibrium were determined in GenAlEx 6.5 (Peakall and Smouse, 2012). Linkage disequilibrium was detected in FSTAT 2.9.4 (Goudet, 1995). MICRO-CHECKER was used to check the presence of null alleles (van Oosterhout et al., 2004).
A total of 5.534 Gbp raw data were obtained after sequencing, and 5.498 Gbp clean data were obtained after quality filtering. Clean reads were assembled into 418,026 contigs. In total, 3,896 SSRs were identified from RAD data, of which 3,727 were used to design primers successfully. There were 353 (9.74%) dinucleotide repeat motifs, 3,199 (85.83%) trinucleotide repeat motifs, 96 (2.58%) tetranucleotide repeat motifs, 34 (0.91%) pentanucleotide repeat motifs and 44 (1.18%) hexanucleotide repeat motifs among these 3,727 SSRs.
Of the 200 randomly selected primer pairs, seventeen microsatellite loci (8.5%) showed good performance with clear peaks in all eight individuals and were selected for further analysis. Finally, fifteen of them showed polymorphism in the four populations (Table 1). For these fifteen microsatellite loci, no linkage disequilibrium was detected between any two loci after sequential Bonferroni correction, indicating that all of the loci were inherited independently. Signs of null alleles were detected in four loci (Bne020, Bne174, Bne145 and Bne150). There were significant deviations from Hardy–Weinberg equilibrium in most of the loci (P < 0.05) except Bne173, Bne076 and Bne199, likely as the consequences of null alleles, the Wahlund effect (Wahlund, 1928), nonrandom mating and an excess of heterozygotes (Table 2). In total, 49 alleles were amplified in four populations. The number of alleles (A) per locus ranged from 2 to 6 with an average of 3.3. HO and HE in each locus ranged from 0 to 1 (mean = 0.437) and from 0 to 0.722 (mean = 0.366), respectively. Across all loci, the four populations of B. nepalensis showed low levels of genetic variation. The average of number of alleles in each population ranged from 2.1 to 2.5; HE varied from 0.335 to 0.396 (Table 2). These values of genetic parameters in B. nepalensis were similar to those reported for other mycoheterotrophic plants, such as Monotropa hypopitys (2.69 for A) (Klooster et al., 2009), Pyrola japonica sensu lato (0.360 for HE) (Shutoh et al., 2017) and Gastrodia flavilabella (0.444 for HE) (Tsai et al., 2014). The low genetic variation might be attributed to genetic drift because of the fragmented and small habitats, as mentioned in the study of M. hypopitys (Klooster and Culley, 2010). However, this seems unlikely in B. nepalensis since higher mean values of HO than of HE in all four populations were found, indicating heterozygote excess, a phenomenon that has been observed in B. wallichii (Zhang and Saunders, 2000). A recent bottleneck caused by habitat destruction and fragmentation may have given rise to the excess of heterozygotes. Additionally, there are no reports about pollinators in Burmannia, except for the description that Burmannia lutescens was visited by Culicidae (mosquitoes) at dawn and during darkness (Kato, 1996), and we failed to observe any pollinators of B. nepalensis in the field. Therefore, limited pollen dispersal may be one contributor to the low level of genetic variation in B. nepalensis. Furthermore, cross-species amplification showed that six loci could be amplified successfully in at least one congeneric species. Five of these loci (Bne110, Bne174, Bne145, Bne150, Bne114) could be amplified in B. wallichii, and only one to three loci in the other three congeneric species (Table 3), indicating a low transferability of these loci in Burmannia. This implies high genetic divergence between B. nepalensis and congeneric species.
| Locus | Primer sequences (5′–3′) | Motif | Allele size range (bp) | Ta (℃) | Fluorescent dye | GenBank accession number |
|---|---|---|---|---|---|---|
| Bne110 | F: GGCTTTTAGAAGGCTTTTCACATA | (CA)7 | 138–144 | Touchdown | HEX | MW489235 |
| R: ACTAAGGGATCCTGAACATCCTC | ||||||
| Bne020 | F: CATACCACTAAAGATGGTGGTCG | (CT)9 | 160–162 | Touchdown | HEX | MW489236 |
| R: TGACTATATAAAGGGGTCCAGCC | ||||||
| Bne076 | F: GATCTGCTTCAAATTCGCTTCT | (AG)6 | 120–122 | Touchdown | FAM | MW489237 |
| R: CCTTTGGAATCCTATAAATGGAGA | ||||||
| Bne197 | F: TGAAAGATGAACGATGAGGAACT | (TCC)7 | 149–164 | Touchdown | ROX | MW489238 |
| R: ATCGAGGAGCTGCCAAAAAT | ||||||
| Bne023 | F: TGGAGATATGGATTGAACTTTCG | (GT)8 | 153–155 | Touchdown | HEX | MW489240 |
| R: GACATTAAAGTTATTCGCGGTCC | ||||||
| Bne199 | F: CGAAGTAGAAACAAAATCGAAAAAC | (CGG)7 | 116–122 | Touchdown | ROX | MW489241 |
| R: TTTTCTCAAACCCTAACTCCTCC | ||||||
| Bne174 | F: ACTATCTTGGAATCACAAGGCAA | (TTG)7 | 131–140 | Touchdown | FAM | MW489242 |
| R: AGACCTTCACCAAAACAACAAAA | ||||||
| Bne113 | F: AGCAAGCAAGCAAGATCATAGAC | (AG)7 | 102–126 | Touchdown | HEX | MW489243 |
| R: TACAAGGGTCAACCATTTACCAG | ||||||
| Bne014 | F: GGCAGAAGAGAAACAAGTGAAAA | (AG)10 | 114–130 | Touchdown | HEX | MW489245 |
| R: CTTCACTCCTTGTTCATGACCTT | ||||||
| Bne145 | F: AGTGTTACCAAAGGGAAGAGTCC | (AGA)10 | 166–175 | Touchdown | ROX | MW489246 |
| R: CGACAACCAAACTTTTATTCAGC | ||||||
| Bne150 | F: ATAGAGCGGAATGGAGAGATAGG | (TTG)9 | 148–154 | Touchdown | ROX | MW489247 |
| R: GCCTTGGTTTCAGATTTCCTAGT | ||||||
| Bne114 | F: GTATGGGACTTATGAATCAAGCG | (GA)7 | 121–125 | Touchdown | HEX | MW489248 |
| R: CATGCTCTCTCTCTGGTTTGAAT | ||||||
| Bne019 | F: CGTCACTACAGGTATTCCCCATA | (CT)9 | 114–130 | Touchdown | HEX | MW489249 |
| R: AACAGACAAAAGCTACCCCAACT | ||||||
| Bne173 | F: AAACCCTAGCTCTGAATCGAAGT | (CGC)7 | 211–232 | Touchdown | FAM | MW489250 |
| R: TTATCTTTAGGGGAAGAAGGGGT | ||||||
| Bne084 | F: CAATGCACTATCTTCAGCCTCTT | (AG)6 | 97–99 | Touchdown | FAM | MW489251 |
| R: AAAAGTTTTTGCTCTCCTCCAAC |
Ta: annealing temperature.
| Locus | CH (n = 7) | XSJ (n = 12) | YC (n = 6) | WYL (n = 10) | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| A | HO | HE | A | HO | HE | A | HO | HE | A | HO | HE | |
| Bne110 | 3 | 1.000 | 0.561 | 3 | 1.000 | 0.594* | 3 | 1.000 | 0.611 | 3 | 1.000 | 0.620* |
| Bne020 | 2 | 0.000 | 0.408* | 1 | 0.000 | 0.000 | 1 | 0.000 | 0.000 | 1 | 0.000 | 0.000 |
| Bne076 | 2 | 0.286 | 0.245 | 2 | 0.083 | 0.080 | 1 | 0.000 | 0.000 | 1 | 0.000 | 0.000 |
| Bne197 | 3 | 0.714 | 0.500 | 2 | 0.167 | 0.153 | 2 | 1.000 | 0.500* | 2 | 0.700 | 0.455 |
| Bne023 | 1 | 0.000 | 0.000 | 1 | 0.000 | 0.000 | 2 | 0.000 | 0.444* | 1 | 0.000 | 0.000 |
| Bne199 | 2 | 0.286 | 0.245 | 2 | 0.250 | 0.219 | 1 | 0.000 | 0.000 | 2 | 0.300 | 0.255 |
| Bne174 | 2 | 0.000 | 0.245* | 3 | 0.000 | 0.542* | 2 | 0.000 | 0.444* | 2 | 0.000 | 0.180* |
| Bne113 | 4 | 0.714 | 0.704* | 5 | 0.917 | 0.684* | 4 | 1.000 | 0.722* | 4 | 0.400 | 0.530* |
| Bne014 | 4 | 0.857 | 0.694* | 4 | 0.667 | 0.708* | 4 | 0.667 | 0.722 | 3 | 1.000 | 0.580* |
| Bne145 | 1 | 0.000 | 0.000 | 3 | 0.083 | 0.594* | 3 | 0.000 | 0.611* | 3 | 0.100 | 0.485* |
| Bne150 | 1 | 0.000 | 0.000 | 2 | 0.083 | 0.413* | 1 | 0.000 | 0.000 | 2 | 0.000 | 0.420* |
| Bne114 | 2 | 0.857 | 0.490* | 2 | 1.000 | 0.500* | 2 | 1.000 | 0.500* | 2 | 1.000 | 0.500* |
| Bne019 | 2 | 0.429 | 0.337 | 3 | 1.000 | 0.594* | 3 | 1.000 | 0.611 | 2 | 1.000 | 0.500* |
| Bne173 | 3 | 0.286 | 0.357 | 2 | 0.167 | 0.153 | 2 | 0.333 | 0.278 | 1 | 0.000 | 0.000 |
| Bne084 | 2 | 0.857 | 0.490* | 2 | 1.000 | 0.500* | 2 | 1.000 | 0.500* | 2 | 1.000 | 0.500* |
| Mean | 2.3 | 0.419 | 0.352 | 2.5 | 0.428 | 0.382 | 2.2 | 0.467 | 0.396 | 2.1 | 0.433 | 0.335 |
A: number of alleles; HO: observed heterozygosity; HE: expected heterozygosity; n: number of individuals genotyped; * signifcant deviation from Hardy–Weinberg equilibrium (P < 0.05).
| Locus | B. coelestis (n = 19) | B. wallichii (n = 17) | B. decurrens (n = 17) | B. oblonga (n = 16) |
|---|---|---|---|---|
| Bne110 | — | 114–140 | — | — |
| Bne020 | — | — | — | — |
| Bne076 | — | — | — | — |
| Bne197 | — | — | — | — |
| Bne023 | — | — | — | — |
| Bne199 | — | — | — | — |
| Bne174 | 256 | 137–146 | 134–140 | — |
| Bne113 | — | — | — | — |
| Bne014 | 215–217 | — | — | — |
| Bne145 | — | 223 | — | 191 |
| Bne150 | 220–223 | 238 | — | — |
| Bne114 | — | 121–127 | 121–127 | — |
| Bne019 | — | — | — | — |
| Bne173 | — | — | — | — |
| Bne084 | — | — | — | — |
— : unsuccessful amplification.
In conclusion, we have developed and characterized fifteen microsatellite markers based on RAD-seq for B. nepalensis, a widely distributed mycoheterotrophic species. This is the first report of microsatellite loci in Burmannia. These markers will be useful tools for studies of population genetic structure, as well as population history dynamics, in B. nepalensis. They should provide opportunities to investigate the potential influence of ecological and demographic factors on genetic patterns of mycoheterotrophic species.
We thank Lianxuan Zhou, Zhongtao Zhao and Xiaojuan Li from South China Botanical Garden, Chinese Academy of Sciences, for plant sampling. We also thank Zhiming Chen of Nanling National Nature Reserve, Guangdong Province and Yajin Luo of Yachang Orchidaceae National Nature Reserve, Guangxi Province for field assistance. This work was supported by the National Natural Science Foundation of China (No. 31970206) and the Biological Resources Program, Chinese Academy of Sciences (Grant No. KFJ-BRP-017-35).
