Cyperus esculentus is an annual herb of Cyperus genus in the family of Cyperaceae, also called “Youshacao” and “Hujianguo” in Chinese. It is native to Africa and the Mediterranean coastal countries and enjoys the reputation of "fruit of life" (http://www.371zy.com/nyjs/zzjs/16288.html). With the characteristics of high quality, high yield and high comprehensive utilization value, Cyperus esculentus is an economic crop integrating grain, oil, animal husbandry and feeding (Lu et al., 2019, Modern Agriculture, (6): 11-13). In 2016, the Ministry of Agriculture issued the National Planting Structure Adjustment Plan (2016-2020). As a new oil crop, Cyperus esculentus was recommended for demonstration and promotion in suitable areas, which is expected to become an important alternative raw material for soybean (http://www.cnfood.cn/hangyexinwen133054.html). In 2017, Zhongyousha 1, the first oil-rich and high-yield C. esculentus variety in China, was selected and bred by Oil Crops Research Institute of CAAS, which is the C. esculentus variety with the highest oil content in the Yangtze River Basin (Zhao et al., 2019, China Seed Industry, (6): 96-97). At present, domestic and foreign studies on Cyperus esculentus mainly focus on cultivation techniques (Zhang et al., 2019), nutritional composition analysis (Manek et al., 2012; Akonor et al., 2019), processing technology (Aljuhaimi et al., 2018; Liu et al., 2020).

RNA-Seq is the use of high-throughput sequencing technology to sequence all RNA reverse transcripts from cDNA libraries to reveal all transcripts expressed in specific tissues or cells and their expression levels. This technology has the advantages of short sequencing time, low cost, high accuracy and high throughput, and has become a common method for studying transcriptome (Cui et al., 2019). RNA-Seq can not only identify sequence variation of genome known species, discover new transcripts and detect fusion transcripts, but also obtain Unigene from sequencing and assembly of genome unknown species for functional annotation, classification, SSR/SNP identification and other bioinformatics analysis. RNA-Seq technology has been widely applied in the study of secondary metabolites of medicinal plants in recent years (Hui et al., 2019; Liu et al., 2019), molecular marker development of herbage (Zheng and Xie, 2020) and stress resistance of forest and grass plants (Gao et al., 2019). Zhang et al. (2018) successfully identified sucrose metabolism pathways in the root tuber of Cyperus esculentus by transcriptome analysis of oil-rich or sugar-rich storage tissues of other plants. Several conserved and significantly expressed genes may be involved in the accumulation of sugar in the root tuber. The genes specifically expressed in the root may be involved in the special molecular mechanism and transcriptional regulation of sucrose metabolism in root tubers. Jing et al. (2015) successfully cloned two segments of PEPCase gene ppc of C. esculentus by homology comparison and analyzed their evolutionary relationship with ppc gene of related species. Gao et al. (2020) found the CeDGAT1 gene of diacylglyceryl transferase 1 by analyzing the transcriptomic information of the tuber of Cyperus esculentus, and identified its function through bioinformatics analysis and qRT-PCR detection, providing a basis for the analysis of the oil biosynthesis mechanism of Cyperus esculentus tuber. In addition, the transcriptome and genome information of C. esculentus is relatively scarce, and the molecular research is still in the initial stage.

In this study, RNA-Seq technology was used to conduct transcriptome analysis on mixed samples of leaves and roots of the Cyperus esculentus. Functional annotation and classification were performed on the obtained Unigenes. And SNP loci in the transcriptome data were analyzed, laying a foundation for subsequent analysis of genetic diversity, gene location and mining and utilization of functional genes of Cyperus esculentus.

1 Results and Analysis

1.1 Transcriptome sequencing and data assembly of C. esculentus

Illumina high-throughput sequencing platform was used for transcriptome sequencing of mixed samples from leaves and roots of Cyperus esculentus. Strict quality control was carried out on the original data. After filtration, 33 988 131 pure sequences were obtained, with a total length of 10.2 Gb, in which GC content was 49.65%. The base mass values Q20 and Q30 were 98.08% and 94.27%, respectively (Table 1), indicating that the sequencing quality of C. esculentus was high and transcriptome data could be used for subsequent analysis.

After obtaining high-quality sequencing data, a total of 22 622 transcripts were assembled using Trinity software, with a total length of 25 852 194 bp, with an average length of 1 142 bp and a length of 1 906 bp for N50. The total length of Unigene was 21 274 744 bp, the average length was 1 085 bp, and the length of N50 was 1 915 bp (Table 2). The N50 lengths of both transcripts and Unigene were much larger than their average lengths, indicating that the assembly effect was good.

1.2 Functional annotation of C. esculentus Unigene

BLAST software was used to compare all Unigene with 7 common databases, and the number of Unigene with annotation information was 12 896, accounting for 65.77% of the total. The number of Unigene annotation varies greatly among different databases, among which Nr database has the largest number of Unigene notes (12 288), accounting for 62.67% of the total. The number of annotations in Swissprot and Pfam databases was 10 531 and 10 005, respectively, with the proportion of annotations exceeding 50%. COG database only has 5 237 Unigene annotated information, accounting for 26.71% of the total, which is the least annotated among all databases. In all the databases compared, the number of Unigene with annotation information≥1 000 nt was much higher than the number of Unigene with annotation information between 300 and 1 000 nt (Table 3).

1.3 GO function classification of C. esculentus Unigene

GO annotation was performed on all Unigene of Cyperus esculentus (Figure 1), and the Unigene successfully annotated could be divided into three categories, that is, cell component, molecular function and biological process, and 52 subcategories according to their functions. There are 16 891 Unigene involved in cell components, including 15 subclasses, including cell, cell part, membrane and organelle, etc. The number of Unigene was 3 500, 3 489, 2 809 and 2 472, respectively, which might be related to cell proliferation during the growth and development of the seedlings. A total of 8 000 Unigene related to molecular functions were divided into 14 subclasses, among which catalytic activity (3 657) and binding (3 298) were the most involved, accounting for 45.71% and 41.23% of the total, respectively. Biological processes are divided into 23 subgroups, including 11 711 Unigene, in which 3 761 Unigene involved in metabolic process and 3 512 Unigene involved in cellular process. The above results indicate that the leaves and roots of C. esculentus grow vigorously in seedling stage, and most Unigene participates in metabolic activities in the body.

1.4 Nr functional annotation of C. esculentus Unigene

The Unigene obtained by the assembly of C. esculentus had the largest number of annotations in the Nr database, and the species distribution of the Unigene homologous sequences annotated by Nr was analyzed (Figure 2). The number of annotated Unigene homologous sequences in Ananas comosus was 3 482 (28.36%). The next largest number of species were Elaeis guineensis (1 530) and Phoenix dactylifera (1 134), accounting for 12.46% and 9.23% respectively. 701 annotated Unigene had homologous sequences in bananas, accounting for 5.71% of the total. In addition, homologous sequences of annotated Unigene were also found in other species, such as 299 (2.43%) wild rice, 296 (2.41%) Panicum hallii, 282 (2.30%) Ensete ventricosum, 280 (2.28%) Asparagus officinalis, 273 (2.22%) Zea mays and 268 (2.18%) Setaria italica. Some Unigene was not annotated in the database due to the lack of genomic and transcriptome information of C. esculentus. The annotated Unigene of Cyperus esculentus has low similarity and poor reference with known genome information, indicating that Cyperus esculentus contains a large number of specific genes.

1.5 KOG functional classification of C. esculentus Unigene

By comparison with KOG database, a total of 7 271 C. esculentus Unigene samples were annotated (Figure 3). According to the annotation function of Unigene, it can be divided into 25 groups, and the number of Unigene in each group is different. Among them, General function prediction only contained 1 349 annotated Unigene, accounting for 16.63% of the total number, which is the largest group. Followed by posttranslational modification, protein turnover and chaperones groups and signal transduction mechanisms with the annotated Unigene of 923 (11.38%) and 747 (9.21%). The above results indicate that the gene expression of C. esculentus seedlings is relatively active, and diverse proteins are formed through different processing and modification after translation, so as to complete their normal life activities.

1.6 SNP characteristics analysis of C. esculentus

A total of 154 456 SNPs were retrieved from the transcriptome sequencing data of C. esculentus (Figure 4). According to the location of mutation locus on Unigene and gene location information on Unigene, the region of mutation locus in the genome and the influence of mutation were obtained. There were 84 843 SNPs, accounting for 56.28% of the total number of SNPs. The number of non-synonymous SNP mutations was 33 983, accounting for 22.54% of the total number. There were 12 821 SNPs (8.5%) in the 5 'UTR region and 18 733 SNPs (12.43%) in the 3' UTR region. In addition, the number of premature terminations or no terminators caused by SNP loci was 304 (0.2 %) and 77 (0.05%), respectively. These results indicated that there were a large number of SNP loci in the transcriptome of C. esculentus, which laid a foundation for the subsequent mining of SNP markers and functional genes linked to specific traits of C. esculentus.

SNP mutations could be divided into six categories, among which C:G>T:A (59 080) and T:A>C:G (51 893) had the highest frequency, 38.25% and 33.60%, respectively. The frequency of the other four types T:A>G:C (9 771), C:G>G:C (12092), T:A>A:T (10 850) and C:G>A:T (10 770) was relatively low, accounting for about 7%. Among all SNPs, the frequency of base transitions (110 973) was 71.85%, 2.55 times higher than that of base transversions (43 483) (Table 4).

2 Discussion

With the continuous recognition of economic value, social value and ecological value of C. esculentus, the planting area of C. esculentus expands rapidly, and the development of C. esculentus industry has a broad prospect. At present, C. esculentus germplasm resources in China are scarce, and varieties are seriously mixed and degraded, which seriously affect the yield and quality and become one of the constraints to the development of C. esculentus industry in China (Wang et al., 2019). Zhang (2019, China Rural Science and Technology, (4): 67-69) proposed that the development goal of C. esculentus in China is to excavate the genes related to drought resistance, oil content and fatty acids of C. esculentus, and use the combination of modern biotechnology such as gene editing and traditional breeding techniques to improve the drought resistance, yield and quality of C. esculentus. RNA-seq technology can effectively excavate and identify the gene expression and function, and complete the annotation and classification of the gene for the no reference gene. In this study, 19 607 Unigenes were obtained by transcriptome sequencing and data assembly. 12 896 Unigene were annotated, accounting for 65.77% of the total. These annotated Unigene can provide reference for studying lipid metabolism, stress response and genetic diversity of C. esculentus. In addition, there were 6 711 Unigene unannotated sequences, which were speculated to be non-coding RNA sequences, short sequences that did not contain a protein functional domain, or new genes specific to C. esculentus.

Analysis of genetic diversity of plant germplasm resources mainly relies on the identification of agronomic traits (Shan et al., 2020, https://kns.cnki.net/kcms/detail/11.1808.S.20200714.1441.004.html), and based on RAPD (Random amplified polymorphic DNA) (Sharma et al., 2018), AFLP (Amplified gragment length polymorphism) (Ipek et al., 2015), SRAP (Sequence-related amplified polymorphism) (Tang et al., 2015), SSR (Simple sequence repeats) (Lu et al., 2018), SNP (Single nucleotide polymorphism) (Li et al., 2018). Previous studies used SRAP markers to construct molecular fingerprints and analyze genetic diversity of 14 C. esculentus resources from different geographical sources. The results showed that the genetic differences of Chinese C. esculentu germplasm resources were large, which was conducive to the introduction and utilization of C. esculentus in different geographical ecological areas and the breeding of new varieties (Zhao and Wei, 2011; Zhao et al., 2013). SNP markers widely exist in plant genomes, with the advantages of rich polymorphism, large number and wide distribution (Yan et al., 2019). SNP molecular markers have been used in genetic diversity analysis, genetic basis analysis, gene mapping of important traits, marker-assisted selection breeding and variety identification. With the development of modern molecular biology and sequencing technology, transcriptome sequencing has become an important data source for SNP sequence analysis and marker development. A total of 40 919 Unigene were obtained by transcriptome sequencing of “Qianjinqiaomai 1”, and 16 152 SNP loci were detected by gene structure analysis, and the distribution of Unigene was not uniform (Zhang et al., 2020, Heilongjiang Animal Science and Veterinary Medicine, (9): 122-126). Serba et al. (2016) obtained 25 894 and 16 979 SNP markers by transcriptomic sequencing of upland genotype VS16 and lowland genotype K5 in switchgrass, respectively. There were 8 915 SNP differences between VS16 and K5, reflecting the polymorphism of ecotype variation. Ma et al. (2016) identified 28 610 SNPs by transcriptome sequencing of upland cotton (Gossypium hirsutum L.) Xuzhou142 and its hairless natural mutant Xuzhou142fl, of which 50% of the mutant Xuzhou142fl SNP overlapped with Gossypium barbadense. Therefore, the mutant Xuzhou142fl may be derived from the interspecific cross between upland cotton Xuzhou142 and Island cotton. In this study, 154 456 SNPS were identified from the transcriptome data of S. oleifera, with a frequency of 1/138 bp. The type of SNP base variation was C:G>T:A, followed by T:A>C:G. Since SNPS are mainly located in coding regions of the genome, they may be related to gene function.

Transcriptome sequencing results of this study revealed that there were abundant gene expressions in C. esculentus, and the annotated information and metabolic pathways of genes were obtained through bioinformatics analysis, laying a foundation for subsequent studies on functional genes involved in oil synthesis and abiotic stress response. In addition, abundant SNP markers were developed in this study to provide information for the subsequent construction of DNA fingerprints of C. esculentus based on SNP technology, which can be used for genetic diversity analysis of different C. esculentus germplasm materials, authenticity and purity identification of varieties, etc., to accelerate the process of breeding and variety protection of C. esculentus. It not only provides a bioinformatics basis for the development and utilization of SNP markers in C. esculentus, but also contributes to the construction of genetic linkage map, development and utilization of functional genes and molecular assisted breeding of C. esculentus.

3 Materials and Methods

3.1 Test materials

Using Cyperus esculentus 19D8 collected by Bayannur Academy of Agricultural and Animal Sciences as the test material, the mixed samples of leaves and roots of single seedlings were taken, frozen in liquid nitrogen and stored at -80℃ for later use.

3.2 RNA extraction, library construction and transcriptome sequencing

Total RNA was extracted by Trizol method with the help of BIOMICS, and the cDNA library was constructed after Nanodrop 2000 detection of RNA purity (OD_260/280>2.0) and Agilent 2100 accurate detection of RNA integrity (RIN>8.5). Sequencing was performed on Illumina NovaSeq 6000 platform after qualified library inspection.

3.3 Transcriptome sequencing data assembly

Strict quality control was carried out on a large number of raw data generated by sequencing, and high-quality sequences were obtained and base quality assessment was carried out, and sequencing data output was counted. Trinity software v2.8.5 was used to assemble the high-quality data obtained by sequencing, filter out the transcripts with low expression, and obtain reliable transcripts. Unigene and transcript length were counted respectively.

3.4 Functional annotation and classification of Unigene

BLAST software (Altschul et al., 1997) was used to compare Unigene sequences with COG, GO, KEGG, KOG, Pfam, Swissprot and Nr databases, and the annotation information of Unigene in each database was counted.

3.5 SNP prediction

The Reads and Unigene sequences of the samples were compared with the RNA-Seq alignment software STAR, and potential SNP sites were identified using GATK4 software.

Authors’ Contributions

CY and SFB were the experimental designer and executor of this study, and completed the first draft of the paper. ZC and WY completed data collation and analysis; GM and TZH participated in experimental design and analysis of experimental results; CY, SFB and YWZ are the architects and principals of the project, supervising experimental design, data analysis, paper writing and revision. All authors read and approved the final manuscript.

Acknowledgments

This study was supported by the Special Project for the Transformation of Scientific and Technological Achievements of Inner Mongolia Autonomous Region (BZC191033), the Youth Innovation Fund of Inner Mongolia Academy of Agricultural and Animal Husbandry Sciences (2020QNJJNO15), and the Variety Resources Protection Project of the Ministry of Agriculture and Rural Affairs (111821301354052249).

Reference

Akonor P.T., Tortoe C., Oduro-Yebocah C., Saka E.A., and Ewool J., 2019, Physicochemical, microstructural, and rheological characterization of tigernut (Cyperus esculentus) Starch, Int. J. Food Sci., 2019: 3830651
https://doi.org/10.1155/2019/3830651
PMid:31275955 PMCid:PMC6582905

Aljuhaimi F., Ghafoor K., Ozcan M.M., Miseckaite O.E., Babiker E., and Hussain S., 2018, The effect of solvent type and roasting processes on physico-chemical properties of tigernut (Cyperus esculentus) tuber oil, J. Oleo Sci., 67(7): 823-828
https://doi.org/10.5650/jos.ess17281
PMid:29877229

Altschul S.F., Madden T.L., Schaffer A.A., Zhang J., Zhang Z., Miller W., and Lipman D.J., 1997, Gapped BLAST and PSI-BLAST: a new generation of protein database search programs, Nucleic Acids Res., 25(17): 3389-3402
https://doi.org/10.1093/nar/25.17.3389
PMid:9254694 PMCid:PMC146917

Cui K., Wu W.W., and Diao Q.Y., 2019, Application and research progress on transcriptomics, Shengwu Jishu Tongbao (Biotechnology Bulletin), 35(7): 1-9

Gao H.J., Lü X.P., Wang R.J., Ren W., Cheng J.N., Wang Y.P., Shao K.Z., and Zhang J.L., 2019, Application of RNA-seq technology in research on herb, shrub and tree stress resistance, Caoye Xuebao (Acta Prataculturae Sinica), 28(12) :184-196

Gao Y., Chen Y., Sun Y., Wang X.Q., Gao H.L., Xue J.A., and Li R.Z., 2020, Identification and expression analysis of CeDGAT1 gene in Cyperus esculentus with a special oil crop, Shanxi Nongye Kexue (Journal of Shanxi Agricultural Sciences), 48(6): 831-835, 841

Hui F., Liu X.Y., Li Z.Y., Liu F.S., and Yang S.H., 2019, Application of transcriptome sequencing in study of medicinal plants, Zhongcaoyao (Chinese Traditional and Herbal Drugs), 50(24): 6149-6155

Ipek M., Seker M., Ipek A., and Gul M.K., 2015, Identification of molecular markers associated with fruit traits in olive and assessment of olive core collection with AFLP markers and fruit traits, Genet. Mol. Res., 14(1): 2762-2774
https://doi.org/10.4238/2015.March.31.6
PMid:25867425

Jing S.Q., Yu H., Tu Y.X., and Ebeydulla R., 2015, Cloning and sequence analysis of partial ppc gene for PEPCase from Cyperus esculentus L., Shengwu Jishu (Biotechnology), 25(1): 8-12, 16

Li Z.Y., Yu H.L., Fang Z.Y., Yang L.H., Liu Y.M., Zhuang M., Lü H.H., and Zhang Y.Y., 2018, Development of SNP in cabbage and construction of DNA fingerprinting of main varieties, Zhongguo Nongye Kexue (Scientia Agricultura Sinica), 51(14): 2771-2787

Liu H.B., Shangguan Y.N., Pan Y.C., Zhao Z.Q., Li L., and Xu D.L., 2019, Applications of RNA-Seq technology on medicinal plants, Zhongcaoyao (Chinese Traditional and Herbal Drugs), 50(21): 5346-5354

Liu H.M., Yan Y.Y., Liu X.X., Ma Y.X., and Wang X.D., 2020, Effects of various oil extraction methods on the gelatinization and retrogradation properties of starches isolated from tigernut (Cyperus esculentus) tuber meals, Int. J. Biol. Macromol., 156: 144-152
https://doi.org/10.1016/j.ijbiomac.2020.03.252
PMid:32251750

Lu X., Adedze Y.M.N., Chofong G.N., Gandeka M., Deng Z.J., Teng L.H., Zhang X.L., Sun G., Si L.T., and Li W.H., 2018, Identication of high-efficiency SSR markers for assessing watermelon genetic purity, J. Genet., 97(5): 1295-1306
https://doi.org/10.1007/s12041-018-1027-4
PMid:30555078

Manek R.V., Builders P.F., Kolling W.M., Emeje M., and Kunle O.O., 2012, Physicochemical and binder properties of starch obtained from Cyperus esculentus, AAPS PharmSciTech, 13(2): 379-388
https://doi.org/10.1208/s12249-012-9761-z
PMid:22350737 PMCid:PMC3364391

Ma Q.F., Wu M., Pei W.F., Wang X.Y., Zhai H.H., Wang W.K., Li X.L., Zhang J.F., Yu J.W., and Yu S.X., 2016, RNA-Seq-Mediated transcriptome snalysis of a fiberless mutant cotton and its possible origin based on SNP markers, PLoS One, 11(3): e0151994
https://doi.org/10.1371/journal.pone.0151994
PMid:26990639 PMCid:PMC4798417

Serba D.D., Uppalapati S.R., Krom N., Mukherjee S., Tang Y., Mysore K.S., and Saha M.C., 2016, Transcriptome analysis in switchgrass discloses ecotype difference in photosynthetic efficiency, BMC Genomics, 17(1): 1040
https://doi.org/10.1186/s12864-016-3377-8
PMid:27986076 PMCid:PMC5162099

Sharma S., Kumar P., Gambhir G., Kumar R., and Srivastava D.K., 2018, Assessment of genetic diversity in lettuce (Lactuca sativa L.) germplasm using RAPD markers, 3 Biotech, 8(1): 9
https://doi.org/10.1007/s13205-017-1039-4
PMid:29259884 PMCid:PMC5723581

Tang Y.J., Cao W.J., and Wu K.L., 2015, Genetic diversity analysis and molecular identification card construction of Chinese Cymbidium germplasms based on SRAP markers, Zhongguo Nongye Kexue (Scientia Agricultura Sinica), 48(9): 1795-1806

Wang R.Y., Wang X.S., and Xiang H., 2019, A multi-purpose novel oil crop-Cyperus s, Zhongguo Youzhi (China Oils and Fats), 44(1): 1-4

Yan T.J., Pu Y.Y., Zhang W.L., Tian Q., Wang D., Li Q., Dai S., and Ding H.F., 2019, Identification technology of genuineness and purity for common  varieties on SNP markers, Shandong Nongye Kexue (Shandong Agriculture Science), 51(12): 111-119

Zhang B., Liu J., Fan Z.X., Guo F.D., Hou L., Ma D.Y., and Wang Y., 2019, Research progress of introduction experiment and pollution free and high-efficient cultivation techniques of Cyperus esculentus, Shandong Nongye Kexue (Shandong Agricultural Sciences), 51(3): 147-150

Zhang Z.L., Liu D.T., and Ji H.Y., 2018, Sucrose metabolism in developing oil-rich tubers of Cyperus esculentus: comparative transcriptome analysis, BMC Plant Biol., 18(1): 151
https://doi.org/10.1186/s12870-018-1363-9
PMid:30041609 PMCid:PMC6056992

Zhao Y.G., Guo R.X., and Luo L.X., 2013, Construction of DNA fingerprinting and analysis of genetic diversity with SRAP markers for tigernut (Cyperus esculentus L.), Zhiwu Ziyuan Yichuan Xuebao (Journal of Plant Genetic Resources), 14(2): 222-225

Zhao Y.G., and Wei W.L., 2011, Genetic diversity analysis of tigernut (Cyperus esculentus L.) using SRAP markers, Zhongguo Youliao Zuowu Xuebao (Chinese Journal of Oil Crop Sciences), 33(4): 351-355

Zheng Y.Y., and Xie W.G., 2020, Advances in the development of molecular markers of forage grasses based on RNA-seq, Zhongguo Caodi Xuebao (Chinese Journal of Grassland), 42(1): 154-162