1. Introduction

Notopterygium incisum C. C. Ting ex H. T. Chang (Apiaceae) is an endangered perennial herbaceous medicinal plant in China (She and Watson, 2005). This species is distributed mainly in the alpine mountains of southwestern China at altitudes of 2800-5100 m (Zhou et al., 2003). In addition, it has been reported that N. incisum is an economically and ecologically important plant because of its adaptability to various environmental conditions, including cold, drought, and salt stress (Zhou et al., 2003). However, the wild natural resources of N. incisum have been declining rapidly in recent years due to habitat fragmentation and excessive exploitation, as well as a low rate of natural regeneration. Thus, this herb species has been listed in the Regulations on Conservation and Management of Wild Chinese Medicinal Material Resources (Wang and Xie, 2004). Now, N. incisum was also listed as endangered herbal species in the International Union for Conservation of Nature (IUCN) Red List (Wu and Raven, 2005). Therefore, there is an urgent need to implement conservation and management practices for this endangered species.

Most previous studies of N. incisum focused mainly on its phytochemistry (Liu et al., 2012), pharmacological action (Zhang and Shen, 2008), geographical distribution (Ma et al., 2010), systematics (Pu et al., 2000; Yang et al., 2017), phylogeography (Shahzad et al., 2017), and domestic cultivation (Dong, 2010; Fang et al., 2004). However, few studies have obtained information regarding the population genetics and conservation biology of N. incisum based on DNA molecular markers (Sun et al., 2012). Furthermore, studies that have characterized the genetic variability of N. incisum have used a limited number of molecular genetic markers (Shahzad et al., 2017; Yang et al., 2017; Zhou et al., 2010).

In recent years, RNA sequencing (RNA-seq) technology has been used as an effective molecular method for studying the evolution of species, determining differentially expressed genes, and investigating the population dynamics of organisms over time (Gu et al., 2013; Jiang et al., 2015; Rai et al., 2016). Thus, large amounts of transcriptomic and genomic information have been characterized for model and non-model species based on RNA-seq technology (Singh et al., 2016; Terol et al., 2016; Zhang et al., 2015), such as Arabidopsis, Plantago ovata, and Tinospora cordifolia, where these studies demonstrated the complexity of the growth and development processes in higher plants. In addition, transcriptome sequencing can target genomic regions by using related expressed sequence tag (EST) sequences, which is very helpful for identifying new EST-simple sequence repeat (SSR) and single nucleotide polymorphism (SNP) markers (Ma et al., 2016; Wang et al., 2017). These markers can be employed to study molecular evolution, population histories, and speciation processes in endangered species (Chen et al., 2015; Zhou et al., 2016). Moreover, RNA-seq technology can help to identify novel functional genes (Zhang et al., 2016). These techniques are particularly attractive for comparative genetic mapping and studying adaptive evolution in rare and exotic species (Li et al., 2017).

In nature, the maintenance of genetic variation is critical for the long-term conservation of threatened species (Luciano et al., 2016; Turchetto et al., 2016). However, the lack of transcriptome and genomic resources for N. incisum has greatly hindered studies of its genetic variability. In this study, we sequenced transcriptomes for the flowers, leaves, and stems of N. incisum and assembled them by Illumina paired-end sequencing. The main objectives of the present study were: (1) to characterize the transcriptome of N. incisum; (2) to develop novel EST-SSR and SNP molecular markers; and (3) to identify and annotate new functional genes, especially cold-related genes.

2. Materials and methods 2.1. Plant materials

The leaf tissues of three N. incisum natural populations were collected in Qinghai, Shaanxi, and Sichuan provinces in western China (detailed information in Table S1, Fig. 1). In addition, we have re-analyzed the transcriptome datasets of N. incisum (Jia et al., 2017) with different parameters. Leaves, stems, and flowers were sampled from one N. incisum plant in the western region of China. Total RNA was isolated separately from each tissue with an RNeasy Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Equal amounts of RNA purified from each of the three tissues were mixed together to construct cDNA libraries. The detailed transcriptome sequencing strategies were as followed by Jia et al., (2017).

2.2. De novo transcriptome assembly and functional annotation of unigenes

The raw sequencing reads were first filtered by removing invalid reads, including reads with adaptor contamination, with ambiguous "N" bases at a ratio greater than 5% and reads with more than 50% of bases with a quality lower than 20 in one sequence. Transcriptome assembly of N. incisum was performed de novo using clean reads with the Trinity program (Grabherr et al., 2011). The generated transcripts were clustered based on sequence similarities, and the longest transcript in each cluster represented the final assembled unigenes. All unigenes were then aligned with the National Community for Biotechnology Information (NCBI) nonredundant protein (Nr) database, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database using NCBI BLAST 2.2.28 + with an E value cut-off of le-5. We also used an E value cutoff of le-3 for euKaryotic Ortholog Groups (KOG) database. These analyses followed Jia et al. (2017). If the alignment results from different databases were in conflict with each other, the results from nr database were preferentially selected, followed by Swissprot, KEGG and KOG databases.

2.3. Predictions of SSRs and SNPs

SSRs were mined among the unigenes using MicroSatellite software (MISA; http://pgrc.ipk-gatersleben.de/misa/misa.html). The repeat thresholds for mono-, di-, tri-, tetra-, penta-, and hexanucleotide motifs were minimums of 10, six, five, five, five, and five repeats, respectively. SSRs with flanking regions were selected to design the PCR primers using Primer 3 (http://primer3.sourceforge.net/releases.php). SNPs were identified using the SNP calling strategy implemented in the GATK2 program (McKenna et al., 2010) with the default parameters. The binary alignment map (BAM) results were processed to sort and remove duplicated reads using Picard v1.41 and SAMtools v0.1.18. SNPs with quality scores < 30 or < 5 bp from the end of the sequence read were filtered before subsequent analyses.

2.4. Survey of EST-SSR polymorphisms

Genomic DNA was extracted from these samples using the improved CTAB method (Doyle, 1987). PCR amplification was conducted in a thermocycler in a total volume of 10 μL, which contained 5 μL 2 × mixture (0.1 U Taq DNA polymerase, 500 μM dNTP, 20 mM TriseHCl, 100 mM KCl, and 3 mM MgCl₂), 0.3 M of each primer, 10-50 ng DNA template, and 3.4 μL double-distilled H₂O. The thermal program comprised the following conditions: initial denaturation for 4 min at 94 ℃, followed by 35 cycles of 45 s at 94 ℃, annealing for 45 s at 47-60 ℃ (Table 1), and extension for 1 min at 72 ℃, before a final extension for 7 min at 72 ℃. The PCR products were separated electrophoretically using a 10% nondenaturing polyacrylamide gel and then visualized with the silver staining protocol. Band sizes were determined using Quantity One (Bio-Rad Laboratories, Hercules, CA, USA) with PBR322 DNA/MspI as the molecular size standard. For the SSR data sets, we calculated the number of alleles per locus (N_a), observed heterozygosity (H_o), and expected heterozygosity (H_e) with GenAIEx 6.501 (Peakall and Smouse, 2012). In addition, CERVUS 3.0.7 (Kalinowski et al., 2007) was used to estimate the polymorphism information content (PIC) for each SSR primer.

3. Results 3.1. Illumina sequencing and de novo assembly

After data filtering, a total of 56, 759, 848 clean reads were obtained with a Q20 base value (base quality more than 20) of 96.54%. Using the Trinity de novo assembly program, 259, 314 transcripts were produced with a mean length of 791 bp and N50 length of 1247 bp. The transcripts were then clustered into 120, 716 unigenes with total lengths ranging from 200 to 10, 632 bp, and an N50 length of 850 bp. Among these unigenes, the lengths of 83, 336 (69.03%) genes ranged from 200 to 500 bp, those of 19, 624 (16.26%) ranged from 501 to 1000 bp, and 17, 756 (14.71%) were more than 1000 bp in length (Fig. 2).

3.2. Functional annotation

Sequence similarity searches were performed against public databases to obtain functional annotations of the assembled unigenes. Among 120, 716 unigenes, 65, 254 (54.05%), 24, 060 (19.93%), and 43, 795 (36.27%) consensus sequences shared homology with sequences in the Nr, Nt, and Swiss-Prot databases, respectively. In total, 70, 245 (58.19%) unigenes were successfully annotated in at least one database. The coding region sequences of 65, 965 (54.64%) unigenes were extracted according to the BLAST and ESTSCAN results, and the length distribution of CDS as showed in Fig. S1.

Based on the Nr annotations, 47, 533 unigenes were assigned GO terms using Blast2GO. Interestingly, 47 of the most highly enriched terms were related to biotic and abiotic stresses (Fig. 3). We found that 21 unigenes were significantly related to cold tolerance, five of which were highly orthologous to those in the model plant Arabidopsis. Thus, the unigenes comprising comp31895, comp42622, comp62766, comp356394, and comp77053 were homologous to cold-regulated gene (COR47) (Puhakainen et al., 2004), GRF1- interacting factor 2 (GIF2) (Lee and Kim, 2014), P-glycoprotein 20 (PGP20) (Cho and Cho, 2013), DNA helicase (RECQ4A) (Schröpfer et al., 2014), and UDP-glycosyltransferase superfamily protein (UGT80B1) (Stucky et al., 2015) in Arabidopsis, respectively. These genes are possibly related to cold acclimation, flavonoid biosynthesis, and sterol metabolic processing in alpine plants (Cho and Cho, 2013; Lee and Kim, 2014; Puhakainen et al., 2004; Schröpfer et al., 2014; Stucky et al., 2015).

The putative proteins annotated using KOG were functionally classified into categories such as cellular structure, molecular processing, and signal transduction, where 21, 343 unigenes were assigned to 26 molecular families. Among these categories, the cluster of general function prediction (3431; 16.08%) was the largest group, followed by posttranslational modification, protein turnover, chaperones (2972; 13.92%), and signal transduction mechanisms (1870; 8.76%). By contrast, only a small number of unigenes were assigned to cell motility and unnamed protein functions (eight and two unigenes, respectively) (Fig. 4).

Pathway-based analysis can help to understand the interactions among genes and metabolic biological functions. Based on the KEGG database, 20, 654 unigenes were assigned to five major categories associated with 264 KEGG pathways. In particular, 9839 (47.64%) unigenes were involved with metabolism pathways, which was the category with the highest number of unigenes, while 5475 (26.51%) unigenes were assigned to genetic information processing, 3355 (16.33%) to organismal systems, 2229 (10.79%) to cellular process, and 1899 (9.19%) to environmental information processing (Fig. S2).

3.3. Development and characterization of SSRs and SNPs

Based on the assembled unigenes, 11, 498 sequences containing 13, 149 potential microsatellite loci were identified, where 1431 of these sequences contained more than one EST-SSR and 733 ESTSSRs were present in compound form, such as dinucleotidetrinucleotide compound motifs, and dinucleotide-tetranucleotide compound motifs and so on. Among the microsatellites, dinucleotide motifs were the most abundant motifs (6749; 51.33%), followed by mononucleotide (3271; 24.88%), trinucleotide (2941; 22.37%), tetranucleotide (144; 1.1%), pentanucleotide (22; 0.17%), and hexanucleotide (22; 0.17%) motifs. Among the SSRs identified, the most abundant repeat was A/T (3161; 96.64%). AT/AT (3835; 59.19%) was the most common dinucleotide SSR. Among the trinucleotide repeats, AAG/CTT (834; 28.36%) was the most abundant, followed by ATC/ATG (462) and AAT/GTT (420) repeats (Fig. 5).

Furthermore, 85, 681 SNPs putatively connected with genes were detected in the assembled unigenes. The SNPs occurred in noncoding and coding sequences at frequencies of 64.33% (55, 115) and 35.67% (30, 566), respectively. Among the SNPs in coding regions, synonymous SNPs occurred at a frequency of 35.54% (30, 452).

3.4. Identification of polymorphic markers

Ninety-six primer pairs were randomly selected from 13, 149 microsatellites to evaluate their polymorphisms in three natural populations of N. incisum. Nineteen of 96 primer pairs yielded stable and polymorphic DNA fragments in N. incisum (Table 1). The number of alleles (N_a) ranged from three to nine, the observed heterogeneity (H_o) varied from 0 to 0.472, and the expected heterogeneity (H_e) varied from 0.323 to 0.689. The polymorphic information content (PIC) values ranged from 0.535 to 0.829 (Table 2).

4. Discussion

In this study, we characterized N. incisum transcriptomes by using RNA-seq to identify a large number of SSR and SNP markers, which will facilitate molecular evolution and conservation biology studies for this endangered herb species. In addition, a number of unigenes were assigned to a large range of GO categories, KOG classifications, and KEGG pathways (Figs. 3 and 4, and S1), thereby suggesting that the assembled unigenes represent a wide diversity of transcripts. Interestingly, we found that 21 unigenes corresponded to cold tolerance and five (comp31895, comp42622, comp62766, comp356394, and comp77053) of them had homologs (COR47, GIF2, PGP20, RECQ4A, and UGT80B1, respectively) in Arabidopsis. For example, COR47 belongs to the dehydrin protein family and it contributes to freezing stress tolerance in plants (Puhakainen et al., 2004). GIF2 is a GIF family member with an essential role in the control of cell proliferation in the leaves (Lee and Kim, 2014). PGP20 is related to stress resistance in plants (Cho and Cho, 2013). RECQ4A is generally considered to have critical roles in the regulation of homologous recombination and DNA repair (Schröpfer et al., 2014). In addition, UGT80B1 may be involved with the synthesis of steryl glucosides in Arabidopsis thaliana (Stucky et al., 2015). N. incisum is distributed mainly in high alpine mountains and the low-temperature regions of west China, so these homologous genes (comp31895, comp42622, comp62766, comp356394, and comp77053) may play important roles in cold and dry adaptation by the alpine plant N. incisum.

Furthermore, 96 EST-SSR primers were randomly selected to assess the genetic diversity of N. incisum. Many of primers that amplified discrete PCR bands from the species did not show polymorphism. After screening, nineteen (17.79%) primer pairs successfully yielded amplicons and the expected PCR products. These newly developed EST-SSR markers will be useful for future analyses of genetic diversity and molecular evolution in N. incisum and related species.

This work was co-supported by the National Natural Science Foundation of China (31470400), Shaanxi Provincial Key Laboratory Project of Department of Education (grant no. 17JS135), the programme for the Key Research and Development Plan in Shaanxi province (grant no. 2018ZDXM-SF-014), the Shaanxi Provincial Education Department Serves Local Special Projects (grant no. 2018JC032) and Public health specialty in the Department of Traditional Chinese Medicine (grants no. 2011-76, 201207002, 2012- 13, 2013-135, 201407002, 2014-76, 2015-78, 2016-44, 2017-66).

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2019.01.001.