Using genetic markers to identify the origin of illegally traded agarwood-producing Aquilaria sinensis trees

Global Ecology and Conservation 22 (2020) e00958

Contents lists available at ScienceDirect

Global Ecology and Conservation journal homepage: http://www.elsevier.com/locate/gecco

Original Research Article

Using genetic markers to identify the origin of illegally traded agarwood-producing Aquilaria sinensis trees Zheng-Feng Wang a, b, c, Hong-Lin Cao a, b, c, *, Chu-Xiong Cai d, Zhang-Ming Wang b, c a

Center for Plant Ecology, Core Botanical Gardens, Chinese Academy of Sciences, Guangzhou, Guangdong Province, 510650, China Key Laboratory of Vegetation Restoration and Management of Degraded Ecosystems, South China Botanical Garden, Chinese Academy of Sciences, Guangzhou, Guangdong Province, 510650, China c Guangdong Provincial Key Laboratory of Applied Botany, South China Botanical Garden, Chinese Academy of Sciences, Guangzhou, Guangdong Province, 510650, China d Dongguan Botanical Garden, Dongguan, Guangdong Province, 523086, China b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 21 November 2019 Received in revised form 3 February 2020 Accepted 3 February 2020

Because of high market demand, agarwood-producing trees are frequently illegally traded in Asia and are usually priced high if they are cut from natural forests. In China, some of these traded Aquilaria sinensis trees are declared to be from natural population(s) and are said to produce agarwood “easily and fast”. To distinguish the origins and genetic differences of natural populations from cultivated ones, we used genetic markers, microsatellites and single nucleotide polymorphisms (SNP) to compare them to trees from both cultivated and known natural origin trees. Our microsatellite results revealed that the illegally traded trees were genetically close to cultivated trees, indicating that they were of cultivated origin. By separating SNPs into genic (in genic regions) and nongenic categories representing functional and non-functional SNPs, our results revealed that the genic SNP markers did not detect more genetic differences between the illegally traded A. sinensis trees and cultivated ones than the nongenic SNP markers did, indicating that they are not functionally discernable from the cultivated trees. Our study revealed that sources labelled as natural by poachers might not have natural origins, which is especially true for agarwood-producing species given their limited natural populations and their long and extensive cultivation history. Our results may reduce the public’s desire for natural agarwood from A. sinensis and other agarwood-producing species and benefit the legal agarwood trade. © 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Keywords: Conservation genetics Genetic diversity Microsatellite SNP RAD-Seq

1. Introduction The poaching of animals and plants is a global problem that threatens biodiversity (Muth and Bowe, 1998; Ervin, 2003; Pattanaik et al., 2008; Sekercioglu et al., 2011). In particular, resources with specific traits are generally most attractive to poachers (Jachmann et al., 1995; Cotterill, 2005; Steinmetz et al., 2010; Gong et al., 2017; Marin et al., 2018), and this is also evident for natural resources (Challender and MacMillan, 2014; Liu et al., 2019a). Poaching from natural, uncultivated resources is more prevalent because traditionally, people believe that individuals grown in nature contain more specific

* Corresponding author. Center for Plant Ecology, Core Botanical Gardens, Chinese Academy of Sciences, Guangzhou, Guangdong Province, 510650, China. E-mail address: [email protected] (H.-L. Cao). https://doi.org/10.1016/j.gecco.2020.e00958 2351-9894/© 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/ licenses/by-nc-nd/4.0/).

2

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

properties or compounds than cultivated ones (Espinoza et al., 2014; Liu et al., 2019a), and thus, the former sell for higher prices than the latter. Flaunting wealth by buying natural-origin products is another way of incentivizing poaching from nature. All these reasons encourage more and more poaching and then make the natural resources for the involved species rarer and rarer. The increases in rarity, in turn, result in increases in the prices for the naturally-collected resources (Challender and MacMillan, 2014). Therefore, supplying the market with cultivated products may not alleviate the poaching pressure on the natural world (Liu et al., 2019a). This is particularly true in agarwood-producing species for which naturallabelled agarwood has a higher market demand than cultivated trees (Mohamed and Lee, 2016; Turjaman et al., 2016). Agarwood-producing species, specifically the Aquilaria species in the Thymelaeaceae family, are primarily distributed in the Indomalesia region (Jim, 2015; Lee and Mohamed, 2016). Among the 21 Aquilaria species, 13 are reported to be known agarwood producers (Lee and Mohamed, 2016). They produce agarwood in their trunks and primary branches due to wounding by worms, lightning- or wind-broken branches, natural microbial or fungal infections or infections that are artificially induced by drilling holes, cutting the bark, and injecting chemicals (Jim, 2015; Azren et al., 2019). To harvest agarwood, people fell trees or cut their branches, which results in the slow growth or death of the trees. Therefore, if harvesting is excessive, it could threaten the natural populations of the specific agarwood species. Furthermore, there are challenges involved in producing agarwood. First, in nature, not all agarwood-producing Aquilaria trees or all their trunks or branches produce agarwood (Jim, 2015; Liu et al., 2018; Azren et al., 2019). Second, it may take up to hundred years for an Aquilaria tree to produce agarwood naturally (Xu et al., 2013; Turjaman et al., 2016; Liu et al., 2018; Ye et al., 2018). Third, artificial methods, including biological inoculations and chemical injections to induce agarwood production are still in development, and the agarwood yield and quality from these approaches require further improvement (Jim, 2015; Chen et al., 2017; Ye et al., 2018; Azren et al., 2019). Finally, people are more in favor of using agarwood from natural trees than from cultivated ones (Jim, 2015; Chen et al., 2017; Ye et al., 2018). Therefore, even though agarwoodproducing Aquilaria species have been extensively cultivated using seedlings and grafting in different countries (Lee and Mohamed, 2016; Yin et al., 2016), the yields and quality of the resulting agarwood are still far from meeting the commercial demand, and the poaching of Aquilaria trees in the nature is still actively on-going despite of strict prohibitions (Jim, 2015). In addition to chopping the trunks or branches of Aquilaria trees that may contain agarwood, if poachers find individuals with a generally small stem size that can potentially produce agarwood, they may even dig out entire trees and sell them on the black market, where they can obtain high prices. Consequently, Aquilaria species have been over-exploited, and their number of natural populations are decreasing rapidly. All Aquilaria species are therefore listed in Appendix II of the Convention on the International Trade in Endangered Species (CITES) regulation. Aquilaria crassna, A. khasiana, A. malaccensis and A. rostrata are in the critically endangered category, A. microcarpa is endangered, A. banaensis, A. beccariana, A. cumingiana, A. filaria, A. hirta, A. rugosa, A. sinensis and A. yunnanensis are vulnerable by the International Union for Conservation of Nature (IUCN) Red List of Threatened Species. A. sinensis is endemic to China, and it is the only agarwood-producing species used for crude drug production in Chinese traditional medicine (Mohamed and Lee, 2016). Its agarwood is called “Chenxiang” in Chinese, and it has been in high demand on the domestic market under the rapid economic development of China in recent years (Yin et al., 2016). At present, some illegally traded A. sinensis trees are grown in Dianbai (DB) county, Guangdong province, China (Fig. 1). To prevent these trees from being stolen again, the peasants who bought them protect them with fences (Fig. S1). Local peasants call them “Qi-nanproducing” A. sinensis, in which “Qi-nan” represents the type of high-quality agarwood. However, at least three different Aquilaria species, such as A. sinensis, A. crassna, and A. malaccensis (synonym of A. agallocha) in China, are claimed to be “Qinan-producing” species (Zhang et al., 2011; Li et al., 2016; Yang et al., 2016). Therefore, calling these illegally traded trees “Qinan-producing” is only a gimmick to obtain high prices on the market. These situations raise two questions for these A. sinensis trees: are they from natural populations, and are they genetically different from cultivated ones? To answer these questions, we located and sampled a total of 16 illegally traded agarwood-producing A. sinensis trees growing in Dianbai county for our DB-ILLEGAL population, and we chose and sampled three cultivated and four natural populations for comparison (Table 1, Fig. 1). The three cultivated populations are large leaf (LL) and small leaf (SL) cultivars in Dongguan city (DG-LL, DG-SL) and one village population in Zhongshan city (ZS-VILLAGE). The four natural populations are the Dinghushan (DHS) population in Zhaoqin city, the Luofushan (LFS) population in Huizhou city, the Kuichong (KC) population in Shenzhen city and the Wuguishan population, also in Zhongshan city (ZS-WGS). According to the buyers of illegally traded A. sinensis trees, poachers might have stolen the trees from Shenzhen city decades ago. However, because Shenzhen is one of the fast-developing cities under Chinese reform and opening policies, its fast urbanization has largely destroyed local A. sinensis populations. Therefore, we could only find and sample one A. sinensis population from the Kuichong reserve in Shenzhen city. We chose the Dongguan and Zhongshan populations to use as cultivated examples for comparison because Dongguan and Zhongshan are traditional A. sinensis plantation areas (Mei, 2016) and are renowned for their high-quality agarwood. Historically, the agarwood produced in Dongguan and Zhongshan was used as a tribute for hundreds of years, ever since the Tang Dynasty (A.D. 618) (Huang et al., 2014; Mei, 2016). In Dongguan, people have historically cultivated two varieties for A. sinensis, large leaf and small leaf, and these varieties are planted in the Dongguan Botanical Garden in Dongguan city. In Zhongshan city, A. sinensis trees that were planted in the past are now only sporadically distributed in villages (Mei, 2016). However, Wuguishan, in Zhongshan city, is a reserve where the A. sinensis is believed to be natural. We did not sample the A. sinensis trees from the other areas of Dianbai county, because these areas had no natural A. sinensis populations. We did not sample cultivated A. sinensis populations in Dianbai county, since Huang et al. (2014) reported that the A. sinensis plantations in Dianbai county and Dongguan city were genetically similar.

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

3

Fig. 1. Map showing the sampling sites for Aquilaria sinensis populations in Guangdong province, China. The blue circle and arrow indicate the possible source population of the DB-ILLEGAL population. Population name abbreviations corresponding to sites and/or cultivars are shown in Table 1, and the same applies to the other figures. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

For our analytical methods, we used microsatellite and genome-wide single nucleotide polymorphism (SNP) markers to conduct a population genetics study. To obtain SNP markers, we first assembled a draft genome of A. sinensis and then used it as a reference and applied the RAD-seq (Restriction-site Associated DNA Sequence) method.

2. Materials and methods 2.1. Sample collection In 2017 and 2018, a total of 207 A. sinensis samples were collected from eight populations in six locations, with DB-ILLEGAL from Dianbai county, DHS from Zhaoqin city, LFS from Huizhou city, KC from Shenzhen city, ZS-WGS and ZS-VILLAGE from Zhongshan city, and DG-LL and DG-SL from Dongguan city (Table 1, Fig. 1). The sample sizes were 16e30 individuals from each population. All the samples are identified as A. sinensis by collectors and also by our authors, H-L Cao and C-X Cai, two taxonomists who have studied A. sinensis for many years and H-L Cao has contributed several A. sinensis specimens to the Herbarium of South China Botanical Garden. We included all the populations (207 individuals) for microsatellite analysis. After examining the microsatellite results, we only included the DB-ILLEGAL, DG-LL, DG-SL and KC populations and reduced the examination to 12, 10, 10 and 7 individuals each, respectively, for SNP analyses due to budget limitations (Table 1). For the microsatellite analysis, we collected two or three leaves from each tree and placed them inside plastic bags containing silica gel. For the SNP analysis, we collected two or three leaves and immediately placed them in liquid nitrogen until DNA extraction.

2.2. Microsatellite genotyping and analysis Eighteen microsatellite loci (Wang et al., 2018b) were used to genotype all the A. sinensis individuals (Table 1). For each population, the genetic diversity parameters of the observed and unbiased expected heterozygosity (HO, HE) were estimated using GenAlEx 6.501 (Peakall and Smouse, 2012). The genetic relationships among the populations were visualized by principal component analyses (PCA) in Adegenet 2.1.1 (Jombart, 2008) in the R package. Furthermore, both a model-based method in STRUCTURE 2.3.4 (Pritchard et al., 2000) and a non-model-based method by discriminant analysis of principal components (DAPC) (Jombart et al., 2010) were used to investigate the genetic structures among populations. For the STRUCTURE analysis, assuming the use of both the admixture and correlated allele frequency models, 20 independent runs were performed for each possible population cluster (K, from 1 to 8) using 5
105 iterations after a burn-in period of 5
105 on the total multi-loci genotypes without prior concern regarding the origin of their populations. The choice of the probable K value was made as recommended in the STRUCTURE user’s manual and by the DK method (Evanno et al., 2005). For the inferred K, we used CLUMPP 1.1.2 (Jakobsson and Rosenberg, 2007) to calculate the mean membership

4

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

Table 1 Sampling information and genetic diversities of eight Aquilaria sinensis populations in Guangdong Province, China. Locations

Populations

Abbreviations for populations Population status N (for microsatellite analysis) N (for SNP analysis) HO

HE

Dianbai county Zhaoqin city Huizhou city Shenzhen city Zhongshan city Zhongshan city Dongguan city Dongguan city

Illegally traded Dinghushan Luofushan Kuichong Wuguishan Village Large leaf cultivar Small leaf cultivar

DB-ILLEGAL DHS LFS KC ZS-WGS ZS-VILLAGE DG-LL DG-SL

0.509 0.579 0.401 0.475 0.519 0.588 0.510 0.524

e Natural Natural Natural Natural Cultivated Cultivated Cultivated

Overall*

16 27 27 25 30 30 26 26

12 e e 7 e e 10 10

0.458 0.638 0.419 0.489 0.548 0.566 0.436 0.491

207

39

0.506 0.513

N: sample size; HO: observed heterozygosity; HE: unbiased expected heterozygosity. Note: The HO and HE were calculated for the microsatellite loci.

coefficient of 20 runs per individual. All the STRUCTURE analyses were performed using a combination of functions from StrataG 2.1 (Archer et al., 2017) in the R package; DAPC was performed using Adegenet 2.1.1 in the R package. 2.3. Reference genome assembly We used one DB-ILLEGAL tree to assemble a draft genome for A. sinensis by both Illumina- and PacBio-based shotgun sequencing. The Illumina shotgun sequencing employed a 150-bp paired-end library generated on the Illumina HiSeq X Ten platform and PacBio shotgun sequencing used the PacBio Sequel system. We obtained a total of 462,090,393 pairs (approximately 139G) in Illumina reads (NCBI SRA database accession number SRR8892930) and approximately 16G in PacBio reads (NCBI SRA database accession number SRR8892931). Both the Illumina and PacBio library preparation and sequencing were conducted by Annoroad Gene Technology Co. in Beijing, China. For the quality trimming of the Illumina shotgun paired-end reads, Sickle v. 1.33 (Joshi and Fass, 2011) was used, retaining base quality values larger than 30 and removing reads with lengths of less than 80 bp. RECKONER v. 1.1 was applied (Długosz and Deorowicz, 2017) for read error correction. Based on the trimmed and error-corrected shotgun reads, GenomeScope (Vurture et al., 2017) was used to estimate the genome size and level heterozygosity rate in A. sinensis with the best k-mer estimated from Kmergenie v. 1.7044 (Chikhi and Medvedev, 2014). The quality-trimmed and error-corrected Illumina shotgun paired end reads were also used to correct the errors in the PacBio long reads, using HALC v. 1.1 (Bao and Lan, 2017). A hybrid de novo assembly strategy was used to construct the A. sinensis reference genome for SNP calling. Prior to assembly, the quality-filtered Illumina shotgun reads were first self-assembled by Platanus v. 1.2.4 (Kajitani et al., 2014). This process included three parts, namely contig assembly, scaffold construction and gap closure. We used bubble crush 0.2 for the parameter maximum difference, branch cutting 0.3 for maximum difference, and default values for the other parameters in the contig assembly; bubble crush 0.2 for the maximum difference and the default values for the other parameters in scaffold construction; and the default parameter values for gap closure. After the gap closure, the scaffolds with undefined bases “Ns” (gaps) were split at the gaps. Then, using all the contigs without “Ns” as backbones, DBG2OLC (Ye et al., 2016) was used to overlap the PacBio reads (hybrid assembly). At this step, we used a k-mer length of 22, an adaptive k-mer threshold for each contig of 0.01, a k-mer-matching threshold for each contig of 4, minimum matching k-mers for each two reads of 100, and chimera removal option 1. After the overlapping, a consensus step by DBG2OLC was used to finally call the hybrid-assembled contigs. We used ra (https://github.com/rvaser/ra), a fast assembler for long PacBio reads, to perform the self-assembly for the PacBio shotgun reads. Then, we applied Quickmerge (Chakraborty et al., 2016) to merge the hybrid-assembled contigs with the PacBio self-assembled contigs. The merged contigs were then improved by BAUM (Wang et al., 2018a), followed by polishing with Racon v. 1.3.2 (Vaser et al., 2017) and Pilon v. 1.23 (Walker et al., 2014). For the polishing step, Racon and Pilon were each run three times. The assembled genome was further improved by Fsfix (Teh et al., 2017) to fix some insertion or deletion errors. Finally, we reduced the duplications in the genome due to heterozygosity using purge_haplotigs (Roach et al., 2018). The final assembled genome (NCBI accession number SMDT00000000) was evaluated by Benchmarking Universal ~o et al., 2015). Single-Copy Orthologs (BUSCO) v. 3.1.0, using the plant dataset eudicotyledons_odb10 (Sima 2.4. SNP calling To obtain the SNPs, we used the RAD-seq method (150 bp paired-end) (Wang et al., 2018b) to sequence 39 individuals from four populations (Table 1). Subsequently, we obtained 5G-12G reads (NCBI SRA database accession numbers SRR8885786SRR8885824) for each individual; RAD-seq library preparation and sequencing were conducted by Annoroad Gene Technology Co. in Beijing, China. Using the above assembled genome as a reference, we applied dDocent v. 2.7.6 (a bash pipeline for RAD sequencing) (Puritz et al., 2014) to call and filter the SNPs in our RAD sequencing data. In dDocent pipeline, we used yara v. 1.0.2 (Siragusa et al., 2013) to map the reads. After dDocent SNP calling, we used Plink v. 1.9 (Purcell et al., 2007) to prune the SNPs that

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

5

showed deviations from Hardy-Weinberg equilibrium (HWE) and high-linkage disequilibrium (LD). For LD pruning, we used the threshold value of 0.5, as calculated with the squared correlation value between phased alleles for all pairwise SNP pairs by using the “inter-chr” command in Plink. The pruned SNP datasets were used for subsequent analyses.

2.5. Genic SNP annotation To determine whether the DB-ILLEGAL population was genetically different from the cultivated one, we classified the SNPs into genic- or nongenic-based ones if they were in the gene sequences. To this end, we first masked the repetitive elements for our assembled genome using RED v. 2.0 (Girgis, 2015). Then, we downloaded the A. sinensis RNA-seq reads (NCBI SRA for database accession numbers SRR3043094, SRR3043096, SRR3043098, SRR3095921, SRR3303973, SRR3303974 and SRR3303975) from the GenBank website and quality-trimmed the reads with fastp v. 0.19.8 (Chen et al., 2018b). Subsequently, by integrating the RNA-seq data, the gene sequences were predicted in the A. sinensis repetition masked genome by Augustus v. 3.3.2 (Stanke et al., 2006) with a training set from Arabidopsis thaliana. Finally, we identified the genic SNPs if they were in the exon regions of the genes.

2.6. SNP analysis For both genic and nongenic SNPs, we performed four analyses to examine the genetic relationships among individuals. First, we visualized them by PCA with SNPRelate v. 1.17.1 (Zheng et al., 2012). Second, we used assignPOP v. 1.1.5 (Chen et al., 2018a) for sample discrimination by measuring the individual membership probabilities among the populations based on a machine-learning framework. This approach essentially included three steps. Step one uses PCA for dimensionality reduction. Step two uses Monte-Carlo cross-validation to estimate the mean and variance of the assignment accuracy. For this step, we used the following parameters: train of individual proportion: 0.5, 0.7, and 0.9; train of loci proportion: 0.1, 0.25, 0.5, and 1; loci sample method: random; iterations: 100; and model: support vector machine (svm). Step three uses K-fold crossvalidation to estimate the membership probability. For this step, we used the following parameters: k-fold: 4e6; train loci proportion: 1; loci sample method: random; and model: svm. After the analysis, each tested individual was automatically assigned to the most likely source population by assignPOP according to the highest membership probabilities. Third, we analyzed the SNP dataset using STRUCTURE 2.3.4 with the same parameter values as in the microsatellite analysis above to investigate the genetic structures among populations. Considering only the four populations analyzed here, we did not determine the probable optimal K value but rather illustrated all the clusters from two to four instead. Finally, the relationships among individuals were investigated using the SNP phylogenetic tree generated via SNPhylo v. 20160204 (Lee et al., 2014), with 1000 bootstrap replicates.

3. Results 3.1. Microsatellite data Among the populations, DHS had the highest HO, while ZS-VILLAGE had the highest HE and LFS had the lowest HO and HE values (Table 1).

Fig. 2. Principal components analysis of the microsatellite loci data for Aquilaria sinensis individuals from different populations in Guangdong Province, China. Each plot shows two principle components: PC1 vs. PC2, PC1 vs.PC3 and PC2 vs. PC3, with the percentages of variance explained in brackets. The same type of symbols indicate the individuals from the same populations.

6

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

In PCA, the first three PCs accounted for 10.155%, 7.845%, and 6.417% of the genetic variance among the A. sinensis individuals (Fig. 2). Based on these findings, the DB-ILLEGAL population could be distinguished from the other populations, except the two DG populations. In the DAPC analysis, based on 50 retained principle components, the Bayesian Information Criterion (BIC) values (Fig. S2) indicated that K ¼ 7 was the optimal number of clusters. Using K ¼ 7, the DB-ILLEGAL population was essentially assigned to four genetic clusters, but none of them was distinct or contained only DB-ILLEGAL individuals (Fig. 3). In the STRUCTURE analysis, the mean log-likelihood values indicated that K ¼ 7 was the most probable number of genetic clusters (Fig. S3), while DK showed two obvious peaks, with the highest at K ¼ 2 and the second highest at K ¼ 4. Thus, we plotted the individuals at these three K values (Fig. 3B). At K ¼ 2, the DB-ILLEGAL population was generally grouped with LFS, KC and two DG populations but separated from DHS and two ZS populations. With increasing K values, the DB-ILLEGAL population was clearly divided into different genetic clusters, and some of its individuals showed a genetic admixture of different clusters. 3.2. Reference genome assembly The K-mer analysis, which used Kmergenie on the Illumina shotgun paired-end reads of A. sinensis, indicated that the best k-mer value was 121. GenomeScope used this k-mer size to estimate that the genome size of A. sinensis was 733,623,729 bp, with a level heterozygosity rate of 0.303%. A hybrid de novo genome assembly resulted in 3368 scaffolds, with a total assembly length of 699,794,123 bp, smaller than the estimated genome size by k-mer analysis. The minimum and maximum lengths of the scaffolds were 5543 bp and 1,777,979 bp, respectively, with an average length of 207,777.44 bp and an N50 of 323,621 bp. The BUSCO evaluation indicated that 95.8% (including 91.0% complete and single-copy and 4.8% complete and duplicated) of the 2121 plant core eukaryotic genes were completely captured by the genome assembly, while 1.3% were captured as fragmented and 2.9% were missing. 3.3. SNP data A total of 60,381 SNPs were obtained by dDocent, but only 6492 of them remained after HWE and LD filtering. After annotation, 1305 were found in the exon regions and were considered genic SNPs, while the remaining 5187 were nongenic.

Fig. 3. Plots of the A) DAPC (discriminant analysis of principal components) for genetic clusters K ¼ 7 and B) STRUCTURE at K ¼ 2, 4, and 7 based on microsatellite loci data for different Aquilaria sinensis populations in Guangdong Province, China. In the plots, each vertical bar represents one individual, and its corresponding population is shown under the vertical bars. The colors of the bars show the assignment probability of the clusters. In panel A), the individuals under the black horizontal bars are used in the SNP analysis. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

7

3.3.1. Nongenic SNPs In the PCA, the first three PCs based on nongenic SNPs accounted for 5.413%, 5.128%, and 4.679% of the total variation, respectively. The PC1 and PC2 clearly separated most of the DB-ILLEGAL individuals from the other three populations, but PC3 indicated substantial divergences among the DB-ILLEGAL individuals (Fig. S4). The assignment accuracy results of the overall and individual populations based on assignPOP (Fig. S5) indicated that the proportion of the samples larger than 0.9 and loci larger than 0.1 displayed high discriminatory power for our populations and could distinguish each population correctly, except for the DG-SL population. Then, using the 5187 nongenic SNP loci with a proportion of 1, three k-fold (4-, 5-, and 6-fold) cross validations showed that, although most of the individuals in the four populations could be assigned back to their own populations, they showed clear mixtures of different cluster components (Fig. S6A). In the STRUCTURE analysis (Fig. S6B) at K ¼ 2, most of the individuals appeared intermixed in two clusters, indicating that K ¼ 2 was not optimal. At K ¼ 3 and K ¼ 4, clear divergences appeared among the DB-ILLEGAL individuals, with some showing high assignment probabilities to their own cluster and others showing substantial admixing with DG components. The SNP phylogenetic tree showed that the DB-ILLEGAL population was separated half-half into two branches (Fig. S7). One branch contained six DB-ILLEGAL, all KC, and two DG-SL individuals, but the KC individuals were in an independent subbranch and were clearly separated from the DB-ILLEGAL ones; the other branch contained the remaining DB-ILLEGAL individuals. 3.3.2. Genic SNP In the PCA for the genic SNPs, the first three PCs accounted for 6.086%, 5.233%, and 4.938% of the total variation, respectively. The PCA results indicated that the DB-ILLEGAL individuals were more separated from the KC than from the two DG populations and displayed clear genetic divergences on PC3 (Fig. 4). In the assignPOP analysis, the 1305 genic SNPs showed high discrimination power for all four populations together (overall) and for DB-ILLEGAL, DG-LL and KC separately at sample proportions larger than 0.9 and loci larger than 0.25, but not for the DG-SL population (Fig. S8). All the DB-ILLEGAL individuals were assigned to their own populations, but with a clear mixture of different cluster components in most of the individuals (Fig. 5A). Two DG-SL individuals could be assigned to the DB-ILLEGAL population, but the results were not consistent under different K-fold cross-validations. In the STRUCTURE analysis (Fig. 5B), at K ¼ 2, the DB-ILLEGAL population was clustered together with the DG populations and separated from KC. At K ¼ 3, some divergences appeared among the DB-ILLEGAL individuals, but many of them still clustered together with the two DG populations. At K ¼ 4, more divergences were pronounced among DB-ILLEGAL individuals, with some of them clearly belonging to two clusters (blue and yellow) and others mixing with different clusters. Finally, the phylogenetic tree revealed that most of the DB-ILLEGAL and two DG-SL individuals were clustered into one branch, and two DB-ILLEGAL individuals were excluded from this branch and grouped with one DG-LL individual (Fig. 6).

4. Discussion Poaching natural plants, animals, and their products can negatively affect population sizes and the survival of species that are protected. If these stolen resources are illegally traded, then pinpointing their origins (sites, populations, herds, breeds, or varieties, etc.) and revealing their genetic background (cultivated vs. natural) (Manel et al., 2002; Espinoza et al., 2014; Huang et al., 2014; Hoelzel, 2015; Pimm et al., 2015; Wasser et al., 2015; Ng et al., 2017; Marin et al., 2018; Liu et al., 2019a) are important for local resource conservation. In this study, we used both microsatellite and SNP markers to study the genetic background of illegally traded A. sinensis trees.

Fig. 4. Principal components analysis of the 1305 genic SNP loci data for Aquilaria sinensis individuals from different populations in Guangdong province, China. Each plot shows two principle components, with PC1 vs. PC2, PC1 vs.PC3 and PC2 vs. PC3, and the percentages of the variance explained by PC1, PC2 and PC3 are shown in brackets. Individuals with the same symbol and color are from the same population. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

8

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

4.1. The origin of illegally traded A. sinensis trees Our microsatellite analysis results indicate that the illegally traded A. sinensis trees are more genetically related to the Dongguan cultivars (DG-LL and DG-SL) than to the others. These findings do not support the declaration by the peasants that the illegally traded A. sinensis trees were from Shenzhen city. Or, if the declaration is true, the possible reasons for this finding might be that we did not find and collect samples from the populations where the trees were from, or these populations have become extinct in Shenzhen city. However, the illegally traded A. sinensis trees that were genetically close to those from Dongguan suggested that they most likely descended from a cultivated ancestor. A substantial genetic admixture within these trees also supports this possible cultivated origin because different genetic sources could be introduced to one common area in ancient plantations (Mei, 2016), which is a common phenomenon in crop cultivation (Poets et al., 2015; Hamadeh et al., 2018; Choi et al., 2019). In fact, A. sinensis has a long and extensive cultivation history (Mei, 2016). A. sinensis has been widely cultivated in China since the Tang Dynasty (A.D. 618). Among its plantation areas, due to specific soil and climatic conditions, the agarwood from Dongguan was believed to have the highest quality and was thus commonly used as tribute to the emperors. Therefore, collecting seeds and poaching A. sinensis trees from the Dongguan area were frequent activities in the past (Mei, 2016). These seeds and stolen trees might then have been planted outside Dongguan in South China. Because there are no measures to prevent the seeds of these planted trees from escaping to nature, away from their sources, when their progenies colonize new habitats and grow up in nature, they are easily wrongly treated as “natural”.

Fig. 5. Membership probability results obtained by A) assignPOP via 4-, 5-, and 6-fold cross-validations and B) STRUCTURE for genetic group K ¼ 2, 3, and 4, based on 1305 genic SNP loci for four Aquilaria sinensis populations from Guangdong Province, China. In the plots, each vertical bar represents one individual and its corresponding population is shown under the vertical bars. The colors of the bars show the assignment probabilities to the genetic clusters. In panel A), the individuals assigned to the other populations are indicated by using the other population name written above their bars. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

9

Fig. 6. Maximum likelihood phylogenetic tree constructed using 1305 genic SNP loci for individuals from four Aquilaria sinensis populations in Guangdong province, China. The bootstrap values are placed next to the nodes.

4.2. Functional divergence between illegally traded and cultivated A. sinensis Here, we used SNP markers to compare the genetic differences between illegally traded A. sinensis trees and cultivated ones. However, because no reference A. sinensis genome was available, to obtain the SNP markers, we assembled its genome. Therefore, we will briefly discuss this genome first. According to the flow cytometry analysis, the 1C-value of A. sinensis was 0.94 pg (Farah et al., 2018), which could be converted to a genome size of approximately 919 Mb with 978 Mb/1 pg (Gregory et al., 2007). The estimated and assembled genome sizes from our k-mer analyses were approximately 734 and 700 Mb, respectively, smaller than 919 Mb, implying that more sequencing data are needed to complete our genome. However, unlike most previous RAD-seq analyses that did not employ a reference genome (Catchen et al., 2013) or that used highly fragmented genome assembled by RAD-seq per se (Puritz et al., 2014), in this study, we assembled a reference genome using both short and long reads. Although this assemblage may be incomplete, a BUSCO evaluation indicated that it covered 95.8% of the core genes in eudicotyledons, which indicates that it is less fragmented. Therefore, using this genome as a reference should have provided accurate SNP calling for our analysis.

10

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

Our results from the genic SNP analyses indicate that the illegally traded A. sinensis trees do not form a distinct group and are not genetically different from the cultivated ones. In fact, the nongenic SNPs revealed more divergence between illegally traded A. sinensis trees and cultivated ones than the genic SNPs did (Figs. 4 and 5B vs. Figs. S4 and S6B). Generally, since nongenic SNPs are selectively neutral, their evolution may be primarily influenced by genetic drift; and because genic SNPs are in genic regions, their evolution may primarily be influenced by selection pressure. Hence, the higher genetic divergence between illegally traded A. sinensis trees and the other populations, as illustrated by nongenic SNPs, should reflect the effects of genetic drift, while less genetic divergence between them as indicated by genic SNPs could reflect some functional similarity among illegally traded A. sinensis trees and the other populations due to selection pressure under cultivation. Genetic drift could be attributed to the limited seed dispersal distance in A. sinensis. A. sinensis is primarily pollinated by moths such as pyralids, noctuids, and geometrids, which could contribute to long-distance gene flow, while its seeds are dispersed by both gravity and hornets (Chen et al., 2016). The seeds of A. sinensis are diamond-shaped and relatively large with a length of 14 mm and a width of 5 mm (Liu et al., 2019b), and most of them are directly dispersed around their mother trees (personal observation). For seeds dispersed by hornets, the average dispersal distance was approximately 80 m, with the longest distance being unknown (Chen et al., 2016). However, the longest dispersal distance of A. malaccensis seeds by hornets was only up to 500 m (Manohara, 2013). Therefore, the seed dispersal distance in A. sinensis may be limited. Because the seed flow is the determinant of the genetic structure, even under long-distance pollen flow (Wang et al., 2016), limited seed flow may have resulted in the genetic drift we detected in nongenic SNPs in the DB-ILLEGAL population. 5. Conservation implications Overall, our study indicates that the illegally traded A. sinensis trees are of cultivated origin. Similar to A. sinensis, mass cultivation has been reported for five other major agarwood-producing Aquilaria species in different countries (Lee and Mohamed, 2016). For these species, large-scale plantations can lead to many “natural” populations descended from the cultivated populations, because there are no measures restricting their gene flow in nature (Singh et al., 2015). Furthermore, because the natural populations of agarwood-producing species have declined significantly, poachers may also turn to cultivated trees (Jim, 2015; Subasinghe and Hettiarachchi, 2016), declare them natural, and sell them at high prices on the market. All these results suggest that without genetic background discrimination, intentionally false labelling is possible and will continue in the agarwood market and will promote illegal trade, especially if there is no law to punish or prohibit such labelling. In this case, to reduce or stop the poaching and illegal trade of agarwood species, government agencies and scientific departments should guide the agarwood market to restrict the use of “natural” or “wild” labelling for protected agarwood species and include the genetic backgrounds of cultivated and natural populations of agarwood species in a certificate system for agarwood trade. They should also educate the public to wean themselves from the desire for “natural” agarwood and not to unknowingly become part of the illegal trading practice. Funding This study was supported by Science and Technology Project of Dongguan city (2015108101002); Conservation of traditional Chinese medicine by Forestry Administration of Guangdong; Special foundation for conservation of wild animals, plants and wetlands by Forestry Administration of Guangdong (2130299); Science and Technology Project of Guangdong Province (2017A030303060); Project for conservation and management of wild animals and plants in Guangdong province; Zhongshan Finance Bureau; Guangzhou Wild Life Conservation and Management Office (SYZFCGe[2017]032); and China National Forestry and Grassland Administration (2130211). Declaration of competing interest The authors declare no conflict of interest. Acknowledgements We thank Chuang-Hui Yuan, Shu-Mei Tan, Lin-Fang Wu, Qian-Cai Jiang, Zong-Jian Tan, Pan-Pan Liu, Hao-bing Liao, XiangXu Huang, Dan Lian, SieRu Lai and the others for their helping in field sampling. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.gecco.2020.e00958. References Archer, F.I., Adams, P.E., Schneiders, B.B., 2017. Stratagy: an r package for manipulating, summarizing and analysing population genetic data. Mol. Ecol. Resour. 17, 5e11. https://doi.org/10.1111/1755-0998.12559.

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

11

Azren, P.D., Lee, S.Y., Emang, D., Mohamed, R., 2019. History and perspectives of induction technology for agarwood production from cultivated Aquilaria in Asia: a review. J. Forestry Res. 30, 1e11. https://doi.org/10.1007/s11676-018-0627-4. Bao, E., Lan, L., 2017. HALC: high throughput algorithm for long read error correction. BMC Bioinf. 18, 204. https://doi.org/10.1186/s12859-017-1610-3. Catchen, J., Hohenlohe, P., Bassham, S., Amores, A., Cresko, W., 2013. Stacks: an analysis tool set for population genomics. Mol. Ecol. 22, 3124e3140. https:// doi.org/10.1111/mec.12354. Chakraborty, M., Baldwin-Brown, J.G., Long, A.D., Emerson, J.J., 2016. Contiguous and accurate de novo assembly of metazoan genomes with modest long read coverage. Nucleic Acids Res. 44, e147. https://doi.org/10.1093/nar/gkw654. Challender, D.W.S., MacMillan, D.C., 2014. Poaching is more than an enforcement problem. Conserv. Lett. 7, 484e494. https://doi.org/10.1111/conl.12082. Chen, G., Liu, C., Sun, W., 2016. Pollination and seed dispersal of Aquilaria sinensis (Lour.) Gilg (Thymelaeaceae): an economic plant species with extremely small populations in China. Plant Diversity 38, 227e232. https://doi.org/10.1016/j.pld.2016.09.006. Chen, K.-Y., Marschall, E.A., Sovic, M.G., Fries, A.C., Gibbs, H.L., Ludsin, S.A., 2018a. assignPOP: an R package for population assignment using genetic, nongenetic, or integrated data in a machine-learning framework. Methods Ecol. Evol. 9, 439e446. https://doi.org/10.1111/2041-210X.12897. Chen, S.F., Zhou, Y.Q., Chen, Y.R., Gu, J., 2018b. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, 884e890. https://doi.org/10.1093/ bioinformatics/bty560. Chen, X.Y., Sui, C., Liu, Y.Y., Yang, Y., Liu, P.W., Zhang, Z., Wei, J.H., 2017. Agarwood formation induced by fermentation liquid of Lasiodiplodia theobromae, the dominating fungus in wounded wood of Aquilaria sinensis. Curr. Microbiol. 74, 460e468. https://doi.org/10.1007/s00284-016-1193-7. Chikhi, R., Medvedev, P., 2014. Informed and automated k-mer size selection for genome assembly. Bioinformatics 30, 31e37. https://doi.org/10.1093/ bioinformatics/btt310. Choi, J.Y., Zaidem, M., Gutaker, R., Dorph, K., Singh, R.K., Purugganan, M.D., 2019. The complex geography of domestication of the African rice Oryza glaberrima. PLoS Genet. 15, e1007414 https://doi.org/10.1371/journal.pgen.1007414. Cotterill, F.P.D., 2005. The Upemba lechwe, Kobus anselli: an antelope new to science emphasizes the conservation importance of Katange, Democratic Republic of Congo. J. Zool. 265, 113e132. https://doi.org/10.1017/S0952836904006193. Długosz, M., Deorowicz, S., 2017. RECKONER: read error corrector based on KMC. Bioinformatics 33, 1086e1089. https://doi.org/10.1093/bioinformatics/ btw746. Ervin, J., 2003. Rapid assessment of protected area management effectiveness in four countries. Bioscience 53, 833e841. https://doi.org/10.1641/00063568(2003)053[0833:RAOPAM]2.0.CO;2. Espinoza, E.O., Cady, A.L., Kreitals, N.M., Hata, M., Cody, R.B., Blanchette, R.A., 2014. Distinguishing wild from cultivated agarwood (Aquilaria spp.) using direct analysis in real time and time of-flight mass spectrometry. Rapid Commun. Mass Spectrom. 28, 281e289. https://doi.org/10.1002/rcm.6779. Evanno, G., Regnaut, S., Goudet, J., 2005. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611e2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x. Farah, A.H., Lee, S.Y., Gao, Z., Yao, T.L., Madon, M., Mohamed, R., 2018. Genome size, molecular phylogeny, and evolutionary history of the tribe Aquilarieae (Thymelaeaceae), the natural source of agarwood. Front. Plant Sci. 9, 712. https://doi.org/10.3389/fpls.2018.00712. Girgis, H.Z., 2015. Red: an intelligent, rapid, accurate tool for detecting repeats de-novo on the genomic scale. BMC Bioinf. 16, 227. https://doi.org/10.1186/ s12859-015-0654-5. Gong, S.P., Shi, H.T., Jiang, A.W., Fong, J.J., Gaillard, D., Wang, J.C., 2017. Disappearance of endangered turtles within China’s nature reserves. Curr. Biol. 27, R170eR171. https://doi.org/10.1016/j.cub.2017.01.039. Gregory, T.R., Nicol, J.A., Tamm, H., Kullman, B., Kullman, K., Leitch, I.J., Murray, B.G., Kapraun, D.F., Greilhuber, J., Bennett, M.D., 2007. Eukaryotic genome size databases. Nucleic Acids Res. 35 (Database issue), D332eD338. https://doi.org/10.1093/nar/gkl828. Hamadeh, B., Chalak, L., d’Eeckenbrugge, G.C., Benoit, L., Joly, H.I., 2018. Evolution of almond genetic diversity and farmer practices in Lebanon: impacts of the diffusion of a graft-propagated cultivar in a traditional system based on seed-propagation. BMC Plant Biol. 18, 155. https://doi.org/10.1186/s12870018-1372-8. Hoelzel, A.R., 2015. Can DNA foil the poachers? Science 349 (6243), 34e35. https://doi.org/10.1126/science.aac6301. Huang, J.-X., Liu, X.-K., Ye, Y.-C., Liu, S.-S., Mo, L.-J., Zhuang, X.-Y., 2014. Leaf morphology and genetic diversity of Aquilaria sinensis in Dongguan and its surrounding areas. Guangdong Agric. Sci. 41 (3), 153e158. https://doi.org/10.3969/j.issn.1004-874X.2014.03.037. Jachmann, H., Berry, P.S.M., Imae, H., 1995. Tusklessness in African elephants: a future trend. Afr. J. Ecol. 33 (3), 230e235. https://doi.org/10.1111/j.1365-2028. 1995.tb00800.x. Jakobsson, M., Rosenberg, N.A., 2007. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23, 1801e1806. https://doi.org/10.1093/bioinformatics/btm233. Jim, C.Y., 2015. Cross-border itinerant poaching of agarwood in Hong Kong’s peri-urban forests. Urban For. Urban Green. 14, 420e431. https://doi.org/10. 1016/j.ufug.2015.04.007. Jombart, T., 2008. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24, 1403e1405. https://doi.org/10.1093/bioinformatics/btn129. Jombart, T., Devillard, S., Balloux, F., 2010. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 11, 94. https://doi.org/10.1186/1471-2156-11-94. Joshi, N.A., Fass, J.N., 2011. Sickle: a sliding-window, adaptive, quality-based trimming tool for FastQ files (Version 1.33) [Software]. Available at: https:// github.com/najoshi/sickle. accessed 3 March 2018. Kajitani, R., Toshimoto, K., Noguchi, H., Toyoda, A., Ogura, Y., Okuno, M., Yabana, M., Harada, M., Nagayasu, E., Maruyama, H., Kohara, Y., Fujiyama, A., Hayashi, T., Itoh, T., 2014. Efficient de novo assembly of highly heterozygous genomes from whole-genome shotgun short reads. Genome Res. 24, 1384e1395. https://doi.org/10.1101/gr.170720.113. Lee, S.Y., Mohamed, R., 2016. The origin and domestication of Aquilaria, an important agarwood-producing genus. In: Mohamed, R. (Ed.), Agarwood: Science behind the Fragrance. Springer, pp. 1e20. Lee, T.-H., Guo, H., Wang, X., Kim, C., Paterson, A.H., 2014. SNPhylo: a pipeline to construct a phylogenetic tree from huge SNP data. BMC Genom. 15, 162. https://doi.org/10.1186/1471-2164-15-162. Li, X.Y., Zeng, B.S., Qiu, Z.F., Liu, Y., Fan, C.J., 2016. Induction of adventitious roots during tissue culture of Aquilaria crassna Pierre. J. Cent. South Univ. Forestry Technol. 36 (10), 6e11. https://doi.org/10.14067/j.cnki.1673-923x.2016.10.002. Liu, H., Gale, S.W., Cheuk, M.L., Fischer, G.A., 2019a. Conservation impacts of commercial cultivation of endangered and overharvested plants. Conserv. Biol. 33, 288e299. https://doi.org/10.1111/cobi.13216. Liu, P.P., Liao, H.B., Jiang, Q.C., Tan, Z.J., Sun, H.M., Wang, R.J., Cao, H.L., Huang, X.X., Lai, S.R., Liang, D., Wang, Z.F., 2019b. Spatial genetic structure of Aquilaria sinensis in Wuguishan, Zhongshan, Guangdong province. J. Trop. Subtropical Bot. 27 (1), 65e73. https://doi.org/10.11926/jtsb.3917. Liu, Y., Chen, J., Qian, J., Lin, H., Sun, N., Huang, Z., 2018. Evolutionary analysis and structural characterization of Aquilaria sinensis sesquiterpene synthase in agarwood formation: a computational study. J. Theor. Biol. 456, 249e260. https://doi.org/10.1016/j.jtbi.2018.08.006. Manel, S., Berthier, P., Luikart, G., 2002. Detecting wildlife poaching: identifying the origin of individuals with Bayesian assignment tests and multilocus genotypes. Conserv. Biol. 16, 650e659. https://doi.org/10.1046/j.1523-1739.2002.00576.x. Manohara, T.N., 2013. Wasp-mediated seed dispersal in agarwood plant (Aquilaria malaccensis), a critically endangered and overexploited species of North East India. Curr. Sci. 105 (3), 298e299. Marin, J., Rivera, R., Varas, V., Cortes, J., Agapito, A., Chero, A., Chavez, A., Johnson, W., Orozco-terWengel, P., 2018. Genetic variation in coat colour genes MC1R and ASIP provides insights into domestication and management of South American camelids. Front. Genet. 9, 487. https://doi.org/10.3389/fgene. 2018.00487. Mei, Q.X., 2016. Fragrant MedicinesdAquilaria Sinensis. China Press of Traditional Chinese Medicine, China.

12

Z.-F. Wang et al. / Global Ecology and Conservation 22 (2020) e00958

Mohamed, R., Lee, S.Y., 2016. Keeping up appearances: agarwood grades and quality. In: Mohamed, R. (Ed.), Agarwood: Science behind the Fragrance. Springer, pp. 149e167. Muth, R.M., Bowe, J.F., 1998. Illegal harvest of renewable natural resources in North America: toward a typology of the motivations for poaching. Soc. Nat. Resour. 11, 9e24. https://doi.org/10.1080/08941929809381058. Ng, E.Y.X., Garg, K.M., Low, G.W., Chattopadhyay, B., Oh, R.R.Y., Lee, J.G.H., Rheindt, F.E., 2017. Conservation genomics identifies impact of trade in a threatened songbird. Biol. Conserv. 214, 101e108. https://doi.org/10.1016/j.biocon.2017.08.007. Pattanaik, C., Prasad, S.N., Reddy, C.S., 2008. Mangalajodi wetland: priority site for conservation. Curr. Sci. 95, 816e817. Peakall, R., Smouse, P., 2012. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research - an update. Bioinformatics 28, 2537e2539. https://doi.org/10.1093/bioinformatics/bts460. Pimm, S.L., Alibhai, S., Bergl, R., Dehgan, A., Giri, C., Jewell, Z., Joppa, L., Kays, R., Loarie, S., 2015. Emerging technologies to conserve biodiversity. Trends Ecol. Evol. 30, 685e696. https://doi.org/10.1016/j.tree.2015.08.008. Poets, A.M., Fang, Z., Clegg, M.T., Morrell, P.L., 2015. Barley landraces are characterized by geographically heterogeneous genomic origins. Genome Biol. 16, 173. https://doi.org/10.1186/s13059-015-0712-3. Pritchard, J.K., Stephens, M., Donnelly, P., 2000. Inference of population structure using multilocus genotype data. Genetics 155, 945e959. Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M.A., Bender, D., Maller, J., Sklar, P., de Bakker, P.I., Daly, M.J., Sham, P.C., 2007. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559e575. https://doi.org/10.1086/519795. Puritz, J.B., Hollenbeck, C.M., Gold, J.R., 2014. dDocent: a RADseq, variant-calling pipeline designed for population genomics of non-model organisms. PeerJ 2, e431. https://doi.org/10.7717/peerj.431. Roach, M.J., Schimdt, S.A., Borneman, A.R., 2018. Purge Haplotigs: allelic contig reassignment for third-gen diploid genome assemblies. BMC Bioinf. 19, 460. https://doi.org/10.1186/s12859-018-2485-7. Sekercioglu, C.H., Anderson, S., Akcay, E., Bilgin, R., Can, O.E., Semiz, G., Tavsanoglu, C., Yokes, M.B., Soyumert, A., Ipekdal, K., 2011. Turkey’s globally important biodiversity in crisis. Biol. Conserv. 144, 2752e2769. https://doi.org/10.1016/j.biocon.2011.06.025. Sim~ ao, F.A., Waterhouse, R.M., Ioannidis, P., Kriventseva, E.V., Zdobnov, E.M., 2015. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210e3212. https://doi.org/10.1093/bioinformatics/btv351. Singh, P., Nag, A., Parmar, R., Ghosh, S., Bhau, B.S., Sharma, R.K., 2015. Genetic diversity and population structure of endangered Aquilaria malaccensis revealed potential for future conservation. J. Genet. 94, 697e704. https://doi.org/10.1007/s12041-015-0580-3. Siragusa, E., Weese, D., Reinert, K., 2013. Fast and accurate read mapping with approximate seeds and multiple backtracking. Nucleic Acids Res. 41, e78. https://doi.org/10.1093/nar/gkt005. Stanke, M., Keller, O., Gunduz, I., Hayes, A., Waack, S., Morgenstern, B., 2006. AUGUSTUS: ab initio prediction of alternative transcripts. Nucleic Acids Res. 34, W435eW439. https://doi.org/10.1093/nar/gkl200. Steinmetz, R., Chutipong, W., Seuaturien, N., Chirngsaard, E., Khaengkhetkarn, M., 2010. Population recovery patterns of Southeast Asian ungulates after poaching. Biol. Conserv. 143, 42e51. https://doi.org/10.1016/j.biocon.2009.08.023. Subasinghe, S.M.C.U.P., Hettiarachchi, D.S., 2016. Gyrinops walla: the recently discovered agarwood-producing species in Sri Lanka. In: Mohamed, R. (Ed.), Agarwood: Science behind the Fragrance. Springer, pp. 89e102. Teh, B.T., Lim, K., Yong, C.H., Ng, C.C.Y., Rao, S.R., Rajasegaran, V., Lim, W.K., Ong, C.K., Chan, K., Cheng, V.K.Y., Soh, P.S., Swarup, S., Rozen, S.G., Nagarajan, N., Tan, P., 2017. The draft genome of tropical fruit durian (Durio zibethinus). Nat. Genet. 49, 1633e1641. https://doi.org/10.1038/ng.3972. Turjaman, M., Hidayat, A., Santoso, E., 2016. Development of agarwood induction technology using endophytic fungi. In: Mohamed, R. (Ed.), Agarwood: Science behind the Fragrance. Springer, pp. 57e72.   Vaser, R., Sovi c, I., Nagarajan, N., Siki c, M., 2017. Fast and accurate de novo genome assembly from long uncorrected reads. Genome Res. 27, 737e746. https:// doi.org/10.1101/gr.214270.116. Vurture, G.W., Sedlazeck, F.J., Nattestad, M., Underwood, C.J., Fang, H., Gurtowski, J., Schatz, M.C., 2017. GenomeScope: fast reference-free genome profiling from short reads. Bioinformatics 33, 2202e2204. https://doi.org/10.1093/bioinformatics/btx153. Walker, B.J., Abeel, T., Shea, T., Priest, M., Abouelliel, A., Sakthikumar, S., Cuomo, C.A., Zeng, Q., Wortman, J., Young, S.K., Earl, A.M., 2014. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PloS One 9, e112963. https://doi.org/10.1371/journal.pone. 0112963. Wang, A., Wang, Z., Li, Z., Li, L.M., 2018a. BAUM: improving genome assembly by adaptive unique mapping and local overlap-layout-consensus approach. Bioinformatics 34, 2019e2028. https://doi.org/10.1093/bioinformatics/bty020. Wang, Z.-F., Cao, H.-L., Cai, C.-X., Guo, Y., Wang, Z.-M., 2018b. Microsatellites development and cross-amplification for Aquilaria sinensis, an endangered agarwood-producing tree. J. Genet. 97 (Suppl. 1), e139ee145. https://doi.org/10.1007/s12041-018-1010-0. Wang, Z.F., Lian, J.Y., Ye, W.H., Cao, H.L., Zhang, Q.M., Wang, Z.M., 2016. Pollen and seed flow under different predominant winds in wind-pollinated and wind-dispersed species Engelhardia roxburghiana. Tree Genet. Genomes 12, 19. https://doi.org/10.1007/s11295-016-0973-3. Wasser, S.K., Brown, L., Mailand, C., Mondol, S., Clark, W., Laurie, C., Weir, B.S., 2015. Genetic assignment of large seizures of elephant ivory reveals Africa’s major poaching hotspots. Science 349 (6243), 84e87. https://doi.org/10.1126/science.aaa2457. Xu, Y., Zhang, Z., Wang, M., Wei, J., Chen, H., Gao, Z., Sui, C., Luo, H., Zhang, X., Yang, Y., Meng, H., Li, W., 2013. Identification of genes related to agarwood formation: transcriptome analysis of healthy and wounded tissues of Aquilaria sinensis. BMC Genom. 14, 227. https://doi.org/10.1186/1471-2164-14-227. Yang, D., Wang, J., Li, W., Dong, W., Mei, W., Dai, H., 2016. New guaiane and acorane sesquiterpenes in high quality agarwood “Qi-Nan” from Aquilaria sinensis. Phytochem. Lett. 17, 94e99. https://doi.org/10.1016/j.phytol.2016.07.017. Ye, C., Hill, C.M., Wu, S., Ruan, J., Ma, Z., 2016. DBG2OLC: efficient assembly of large genomes using long erroneous reads of the third generation sequencing technologies. Sci. Rep. 6, 31900 https://doi.org/10.1038/srep31900. Ye, W., He, X., Wu, H., Wang, L., Zhang, W., Fan, Y., Li, H., Liu, T., Gao, X., 2018. Identification and characterization of a novel sesquiterpene synthase from Aquilaria sinensis: an important gene for agarwood formation. Int. J. Biol. Macromol. 108, 884e892. https://doi.org/10.1016/j.ijbiomac.2017.10.183. Yin, Y., Jiao, L., Dong, M., Jiang, X., Zhang, S., 2016. Wood resources, identification, and utilization of agarwood in China. In: Mohamed, R. (Ed.), Agarwood: Science behind the Fragrance. Springer, pp. 21e38. Zhang, L.J., Zhang, K.H., Liang, Y.N., 2011. Experiment on introduction and cultivation of Aquilaria agallocha Roxb.in Guangdong province. Cent. S. For. Inventory Plann. 30 (3), 65e67. https://doi.org/10.3969/j.issn.1003-6075.2011.03.016. Zheng, X., Levine, D., Shen, J., Gogarten, S.M., Laurie, C., Weir, B.S., 2012. A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics 28, 3326e3328. https://doi.org/10.1093/bioinformatics/bts606.