Phenotype comparison and anthocyanin content determination in potato tubers of cultivar ‘Jin-16’ and ‘Xisen-8’
Tubers appearance of potato cultivar ‘Jin-16’ and ‘Xisen-8’ showed that both of the two cultivars were of long oval shape with shallow bud eyes and smooth skin. The tuber skin and flesh of cultivar ‘Xisen-8’ were observed dark purple (Fig. 1b) compared with those of ‘Jin-16’, which were yellow in color (Fig. 1a).
The anthocyanin contents of tuber flesh in ‘Jin-16’ and ‘Xisen-8’ were detected, respectively. As shown in Fig. 1c, the relative anthocyanin content of tuber flesh in ‘Xisen-8’ was significantly higher than that in ‘Jin-16’ (~ 103.5 fold), indicating that there may exist a series of key regulatory factors in ‘Xisen-8’ which promoted anthocyanin accumulation.
Sequencing and assembly of the lncRNA-Seq datasets
To identify the lncRNAs involved in the process of anthocyanin biosynthesis, the tubers of ‘Jin-16’ and ‘Xisen-8’ were used to construct lncRNA libraries and sequenced on Illumina novoseq6000 platform. After eliminating the low quality reads, Illumina adapters and reads with unidentifiable base information, the clean reads were obtained from each sample and accounted for more than 95% of the raw reads (Fig. 2a). Among these clean reads, the exonic reads made up more than 60%, followed by intergenic reads (less than 35%) and intronic reads (less than 5%) (Fig. 2b).
The filtered reads were subsequently aligned against a potato reference genome (S. tuberosum assembly 4.03) by Hisat2 software[32]. The sequence alignment results of reads mapped to the reference genome were presented in Additional file 1. Finally, 6371 non-coding transcripts corresponding 4376 lncRNAs were identified as the candidate lncRNA datasets based on the analyses of CPC2, PFAM and CNCI softwares (Additional file 2). According to the location relationship between lncRNAs and known mRNAs, lncRNAs could be classified into 4 categories: long intergenic lncRNAs (lincRNAs), antisense lncRNAs, sense overlapping lncRNAs and sense intronic lncRNAs[33]. In this study, lincRNAs contained the largest number with 3461 lncRNAs (79.1%), followed by antisense and sense overlapping lncRNAs with 705 (16.1%) and 210 (4.8%) lncRNAs, respectively. And no sense intronic lncRNAs were identified (Fig. 2c).
Construction of expression profile and Chromosomal localization of differentially expressed lncRNAs in potato tubers
After difference significance analysis, a total of 1421 DE lncRNAs were acquired between ‘Jin-16’ and ‘Xisen-8’ tuber samples. The expression levels of these 1421 DE lncRNAs were shown in Additional file 3. In this study, hierarchical clustering analysis was conducted using the FPKM value of DE lncRNAs to cluster the genes with similar expression patterns, which may play analogical roles or participate in common metabolic pathways (Fig. 3a).
Of those DE lncRNAs, 735 lncRNAs were identifed in both ‘Jin-16’ and ‘Xisen-8’ tuber samples; while 325 and 361 DE lncRNAs were found to be specifically expressed in ‘Jin-16’ (like XLOC_029213 and XLOC_053592) and ‘Xisen-8’ (like XLOC_026093 and XLOC_085889), respectively (Fig. 3b). Compared with ‘Jin-16’ tubers, there were 735 up-regulated and 686 down-regulated DE lncRNAs in ‘Xisen-8’ tubers (Fig. 3c). The identification of these DE lncRNAs is conducive to further investigation on anthocyanin synthesis mechanism in potato tubers.
The distribution of these 1421 DE lncRNAs on potato chromosomes and their expression levels were intuitively displayed in Fig. 3d. These DE lncRNAs were unevenly distributed in all potato chromosomes, which were most distributed on CHR1 (191) and least distributed on CHR2 (84). Other chromosomes, such as CHR0, 4, 6, 7, 8, 9 and 11 contained more than 100 DE lncRNAs; while CHR3, 5, 10, 12 contained less than 100 DE lncRNAs. The exact position of the 1421 DE lncRNAs on chromosomes were presented in Additional file 3. The chromosomal localization of lncRNAs contributes to obtaining complete lncRNA sequences and provides basis for better understanding of their functions.
Functional category of DE lncRNA targets by GO enrichment analysis
The co-located and co-expressed targets of all identified lncRNAs in potato tubers were presented in Additional file 4 and Additional file 5. The GO enrichment analysis of genes co-located and co-expressed with DE lncRNAs were performed to predict the potential functions of DE lncRNAs between ‘Jin-16’ and ‘Xisen-8’. In this study, the top 40 significantly enriched GO terms were shown in Fig. 4 and Additional file 6. Of the GO annotations associated with mRNAs co-located with DE lncRNAs, 22 GO terms were assigned to biological process (BP), 8 and 10 GO terms belonged to molecular function (MF) and cellular component (CC), respectively. These genes were mainly enriched in “cellular metabolic process” (GO:0044237), “phosphorus metabolic process” (GO:0006793), “hydrolase activity, acting on ester bonds” (GO:0016788), “transferase activity, transferring acyl groups” (GO:0016746), “membrane-enclosed lumen” (GO:0031974) and “organelle lumen” (GO:0043233).
For the genes co-expressed with DE lncRNAs, the top 40 significantly enriched GO terms were grouped into BP (12 GO terms) and MF (28 GO terms). And no GO term was classified in CC. In the category of BP, “biological_process” (GO:0008150), “metabolic process” (GO:0008152), and “organic substance metabolic process” (GO:0071704) including 7712 (3.26%), 5905 (2.50%), 4606 (1.95%) and 36 (2.54%) genes respectively, were the predominant GO terms. In MF category, “molecular_function” (GO:0003674), “catalytic activity” (GO:0003824) and “transferase activity” (GO:0016740) were the most representative terms and contained 9381 (3.97%), 5444 (2.30%) and 2115 (0.89%) genes, respectively. The results of GO analysis implied that genes enriched in these GO terms were probably involved in anthocyanin biosynthesis. Therefore, DE lncRNAs were the potential participants in these biological processes, mediating the regulation of anthocyanin metabolism in potato tubers through cis- and trans-regulation of the expression of their target genes.
KEGG pathway enrichment of DE lncRNA targets
To delve into the metabolic pathways in which the DE lncRNAs were involved, the genes co-located and co-expressed with these DE lncRNAs were subjected to the KEGG database. Totally, the co-located genes and co-expressed genes were identified to participate in 121 and 119 metabolic pathways, respectively (Additional file 7). According to enrichment significance, the top 20 enriched biological pathways with p value < 0.05 were presented in Fig. 5. The genes co-located and co-expressed with DE lncRNAs were both mainly involved in “DNA replication”, “Homologous recombination”, “Tyrosine metabolism”, “Fatty acid degradation”, “Carbon metabolism” and “Nucleotide excision repair”.
Except for these common pathways, the co-located targets of DE lncRNAs were primarily enriched in the pathway of “Sesquiterpenoid and triterpenoid biosynthesis”, “Zeatin biosynthesis”, “Fatty acid biosynthesis”, “Biosynthesis of secondary metabolites”, “Starch and sucrose metabolism”, etc. While, the co-expressed targets of DE lncRNAs were related to “Photosynthesis-antenna proteins”, “Plant hormone signal transduction”, “Carotenoid biosynthesis”, “RNA transport” and so on. The results of KEGG enrichment analysis further illustrated that these lncRNAs and their target genes might be involved in these specific biological pathways to mediate anthocyanin biosynthesis in potato tubers.
According to the expression patterns of the target genes involved in “Flavonoid biosynthesis” (sot00941) and “Anthocyanin biosynthesis” (sot00942), a schematic for anthocyanin synthesis pathway was proposed in potato tubers (Fig. 6). In the pelargonidin synthesis pathway, the expression levels of genes encoding CHS, F3H and DFR were significantly up-regulated, while the expression levels of CHI and ANS had no significantly difference between ‘Jin-16’ and ‘Xisen-8’. The down-regulated expression of FLS might reduce kaempferol synthesis and enhance pelargonidin production. Although the biosynthesis of caffeoyl-CoA was inhibited due to the down-regulated genes encoding HQT and C3’H, the significantly up-regulated genes such as CHS, F3’5’H, F3H and DFR might promote the preferred synthesis of delphinidin and cyanidin from caffeoyl-CoA. The DE lncRNAs, which were possibly involved in anthocyanin synthesis by regulating these enzyme genes, were listed in Table 1. Each enzyme gene was regulated by multiple lncRNAs and the same lncRNA could mediate the expression of many different genes.
Table 1
The expression of enzyme genes and their corresponding DE lncRNAs involved in anthocyanin synthesis in ‘Jin-16’ and ‘Xisen-8’ tubers.
|
Gene name
|
Annotation
|
Gene ID
|
Gene expression
|
Corresponding LncRNA ID
|
LncRNA expression
|
Xisen-8
|
Jin-16
|
Xisen-8
|
Jin-16
|
CHS
|
Chalcone synthase 1B
|
PGSC0003DMG400029620
|
0.808
|
0.000
|
XLOC_046144
|
0.943
|
0.163
|
XLOC_038153
|
3.103
|
1.442
|
XLOC_047085
|
1.337
|
0.258
|
F3H
|
Flavanone 3 beta-hydroxylase
|
PGSC0003DMG400003563
|
1.037
|
0.047
|
XLOC_038153
|
3.103
|
1.442
|
XLOC_070549
|
0.000
|
2.344
|
DFR
|
Dihydroflavonol 4-reductase
|
PGSC0003DMG400003605
|
4.644
|
0.012
|
XLOC_100183
|
1.443
|
0.057
|
XLOC_107693
|
10.329
|
64.517
|
XLOC_038153
|
3.103
|
1.442
|
XLOC_047085
|
1.337
|
0.258
|
XLOC_075817
|
32.616
|
18.880
|
F3’5’H
|
Flavonoid 3’,5’-hydroxylase
|
PGSC0003DMG400000425
|
1.785
|
0.387
|
XLOC_038153
|
3.1027
|
1.4422
|
XLOC_087286
|
2.998
|
0.409
|
XLOC_095453
|
4.643
|
0.181
|
XLOC_015689
|
4.722
|
8.089
|
XLOC_063942
|
23.342
|
10.796
|
FLS
|
flavonol synthase
|
PGSC0003DMG400014093
|
0.051
|
0.759
|
XLOC_030593
|
0.000
|
9.902
|
XLOC_085889
|
0.000
|
10.540
|
XLOC_010366
|
0.000
|
5.913
|
XLOC_025110
|
0.002
|
4.598
|
XLOC_066673
|
0.000
|
4.863
|
XLOC_004180
|
48.485
|
7.441
|
XLOC_075621
|
46.023
|
5.710
|
XLOC_058142
|
7.446
|
0.016
|
XLOC_090403
|
0.065
|
4.005
|
XLOC_040817
|
0.036
|
4.267
|
XLOC_004081
|
2.628
|
9.961
|
XLOC_012867
|
0.148
|
6.039
|
XLOC_094680
|
0.315
|
4.110
|
XLOC_107693
|
10.329
|
64.517
|
XLOC_067296
|
2.898
|
9.005
|
XLOC_043939
|
1.460
|
5.937
|
XLOC_089954
|
2.283
|
12.322
|
XLOC_049163
|
4.807
|
7.406
|
XLOC_019471
|
3.488
|
6.027
|
C3’H
|
P-coumaroyl quinate/shikimate 3’-hydroxylase
|
PGSC0003DMG400007178
|
1.546
|
3.476
|
XLOC_047468
|
1.966
|
6.560
|
XLOC_060098
|
16.316
|
2.663
|
XLOC_004180
|
48.4847
|
7.441
|
XLOC_005340
|
42.397
|
8.543
|
XLOC_015288
|
62.074
|
8.770
|
XLOC_107693
|
10.329
|
64.517
|
XLOC_096998
|
8.304
|
1.676
|
XLOC_105729
|
8.198
|
1.047
|
HQT
|
shikimate O-hydroxycinnamoyl-transferase
|
PGSC0003DMG400011189
|
8.448
|
46.876
|
XLOC_029213
|
4.977
|
0.000
|
XLOC_056281
|
2.895
|
6.228
|
XLOC_107693
|
10.329
|
64.517
|
XLOC_093630
|
0.959
|
14.693
|
XLOC_006486
|
4.765
|
1.466
|
Analysis of lncRNA-mRNA co-expression network and identification of anthocyanin synthesis related lncRNAs and genes
The interaction network of lncRNAs and their co-expressed mRNAs was constructed and visualized to investigate the potential regulatory relationship between them and explore the relevant lncRNAs and genes involved in anthocyanin synthesis (Fig. 7). The lncRNAs and mRNAs in the network were represented by squares and circles, respectively. The results illustrated that 310 DE lncRNAs and 352 target mRNAs were contained in the co-expression network, establishing 771 anthocyanin-responsive lncRNA-mRNA pairs (Additional file 8). The intricate interaction between lncRNAs and their co-expressed genes was manifested in the following aspects: on the one hand, the same lncRNA could regulate the expression of different mRNAs, such as XLOC_060098, XLOC_004180, XLOC_030593 and so on. For example, XLOC_060098 was predicted to be correlated with NAC (PGSC0003DMG400011891), ERF7 (Ethylene-responsive transcription factor 7, PGSC0003DMG401013892), UGT (UDP-glucosyltransferase, PGSC0003DMG400011492) and so on, which were previously reported as important mediators of anthocyanin metabolism. On the other hand, the expression of most mRNAs could be influenced by multiple lncRNAs. For instance, CCoAOMT5 (Caffeoyl-CoA O-methyltransferase 5, PGSC0003DMG400006448) was found to be correlated with a series of lncRNAs, like XLOC_032852, XLOC_032892, XLOC_039335, XLOC_047085, XLOC_048038, XLOC_096885 and XLOC_100659. Similarly, Peroxidase (PGSC0003DMG402025083) was also regulated by numerous lncRNAs, such as XLOC_086959, XLOC_005760, XLOC_027548, XLOC_016209 and XLOC_017372. In addition to the interaction between lncRNAs and their trans target mRNAs, lncRNAs could also interact with each other directly or indirectly. As shown in Fig. 6, XLOC_060098 was associated with XLOC_053197, XLOC_096423, XLOC_027380, XLOC_084558 and XLOC_011664. These results provide some basis for further screening candidate lncRNAs and genes to study their functions and regulatory mechanisms in anthocyanin synthesis.
Among these 352 target genes regulated by DE lncRNAs, 309 and 43 genes were up-regulated and down-regulated in ‘Xisen-8’ tubers respectively. The up-regulated genes included the mRNAs encoding CCoAOMT5, AnAT (anthocyanin acyltransferase), NAC2, UGT (UDP-glucose:glucosyltransferase), F3’H and so on, which were recognized as key positive elements in anthocyanin synthesis. While the genes coding for C3H (P-coumarate 3-hydroxylase), ARF (Auxin response factor), SPL (Squamosa promoter binding proteinlike) and HQT (Hydroxycinnamoyl-CoA: quinate hydroxycinamoyl transferase) were down-regulated. These genes have been reported to negatively regulate anthocyanin synthesis. Most genes were observed to interact with XLOC_060098, indicating that XLOC_060098 might be an important regulator participating in various biological processes during anthocyanin synthesis.
The qRT-PCR validation of anthocyanin synthesis-related DE lncRNAs and their targets in potato tubers
In order to validate the expression patterns of potential anthocyanin-associated DE lncRNAs and mRNAs, 6 DE lncRNAs and their corresponding targets, including lncRNAs and mRNAs, were selected from the co-expressed network for qRT-PCR analysis utilizing the same RNA samples as used for RNAseq. The specific primer sequences of the lncRNAs and genes were shown in Additional file 9. The results showed that the expression trends of the DE lncRNAs and their co-expressed targets examined by qRT-PCR were basically consistent with the RNAseq data (Fig. 8a). Although there existed quantitative differences in expression levels between the RNA-seq data and qRT-PCR assay, a high linear correlation (R2=0.8226) was obtained under two different detection methods, suggesting that the RNA sequencing data were credible (Fig. 8b).
As shown in Fig. 8a, the expression levels of XLOC_060098, XLOC_017372, XLOC_038153 and XLOC_046144 showed the same patterns with that of their corresponding targets, indicating that these lncRNAs positively regulated the expression of their targets, respectively. In addition, lncRNAs could also negatively regulate their target genes. For instance, the expression of XLOC_100954 and its target genes ARF8 and SPL witnessed a downward trend in ‘Xisen-8’ compared to ‘Jin-16’, which was in contrast to the expression pattern of Ubiquitin ligase. Similarly, HQT, one of the co-expressed target genes of XLOC_029213, presented a different expression pattern from XLOC_029213 and its other target genes MYC4 (Myelocytomatosis4) and AT (Acyltransferase), whose expression amounts were up-regulated in ‘Xisen-8’. The forward and inverse correlations between lncRNAs and their corresponding target genes revealed that the lncRNAs might participate in anthocyanin biosynthesis by regulating the expression of anthocyanin-related genes positively or negatively. Besides, lncRNAs, such as XLOC_017372 and XLOC_046144, could also regulate other lncRNAs via co-expression to mediate anthocyanin biosynthesis.