A Comparative Transcriptome Analysis of Purple and Yellow Fleshed Potato Tubers Reveals Long Non-coding RNAs and their Targets Functioned in Anthocyanin Biosynthesis

Background: Purple eshed potato tubers accumulate large amounts of anthocyanin content, servicing as functional foods and high-value feedstock. Long non-coding RNAs (lncRNAs) have been reported to play an important role in anthocyanin synthesis by regulating gene expression in various action modes. However, the mechanism underlying anthocyanin accumulation mediated by lncRNAs in underground organs remains unclear. Results: To excavate the differentially expressed lncRNAs (DE lncRNAs) between purple and yellow eshed potato tubers, the transcriptome sequencing was performed and a total of 1421 DE lncRNAs were identied. Gene Ontology and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses showed that the target genes of these DE lncRNAs were involved in diverse biological processes and pathways for anthocyanin biosynthesis, reecting the functional diversity of the corresponding lncRNAs. A lncRNA-mRNA interaction network was created based on their correlation to investigate the regulatory relationship among them. Notably, lncRNAs like XLOC_060098 and XLOC_017372 might contribute to anthocyanin synthesis by targeting the key enzyme genes and transcription factor genes in the pathway. Conclusions: The construction of expression proling of DE lncRNAs and lncRNA-mRNA relationship network is helpful for further unraveling the molecular mechanisms of lncRNAs in anthocyanin synthesis in potato tubers, and provides theory basis for the cultivation of functional potato varieties and the improvement of nutritional quality of other underground crops. In this study, fresh tubers of purple eshed potato (‘Xisen-8’) and yellow eshed potato (‘Jin-16’) were used as experimental materials to excavate lncRNAs related to potato anthocyanin biosynthesis through transcriptome sequencing. The cis- and trans-acting target genes of the DE lncRNAs were predicted via bioinformatics analyses. Totally, 4376 lncRNAs were identied, among which 1421 were signicantly differentially expressed with 735 up-regulated and 686 down-regulated lncRNAs in ‘Xisen-8’. Chromosomal mapping analysis showed that these DE lncRNAs were unevenly distributed on potato chromosomes: the number of lncRNAs on chromosome 1 was the highest, and the number of lncRNAs on chromosome 2 was the lowest. The target genes of DE lncRNAs could be clustered into lots of GO types and enriched into multiple metabolic pathways, thus refecting the functional diversity of their corresponding lncRNAs. The co-expression network and the expression patterns of DE lncRNAs and their target genes related to anthocyanin biosynthesis in potato tubers were analyzed to reveal their potential regulatory relationship. Results showed that most anthocyanin-associated enzyme genes and transcription factor genes were predicted to be regulated by lncRNAs such as XLOC_060098 and XLOC_017372, indicating that these lncRNAs might play crucial roles in anthocyanin metabolism. The results of this study lay the foundation for further elucidation of regulatory mechanism of lncRNAs in anthocyanin biosynthesis and provide new insights for quality improvement strategies of potato and other crops. pairs with The lncRNA-mRNA interaction RNA samples were extracted from the tuber esh of ‘Jin-16’ and ‘Xisen-8’ as described above. Quantitative real-time PCR (qRT-PCR) was performed with the TB GreenTM Premix Ex TaqTM (Tli RNase H Plus) (Takara, Dalian, China) on CFX96 PCR System (Bio-Rad, USA). The specicity of these primers which designed by Primer-Blast in NCBI website (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) was tested using dissociation curve analysis. The 10 µl reaction volume samples, containing 5 µL TB Green, 1 µL diluted template, 0.4 µL 10 µM solution of each primer and 3.2 µL ddH 2 O, were used for PCR using the following cycling program: 95 ℃ for 3 min, followed by 40 cycles of 95 ℃ for 10 s, 60 ℃ for 30 s, and 72 ℃ for 20 s. The relative expression of selected lncRNAs and genes was calibrated against the reference gene EF1α using the method of 2 - ∆∆ Ct [1]. For each sample, three biological repeats and two experimental replicates were performed to make sure the results reliable.

In recent years, several lncRNAs have been found to take part in plant growth, ower and fruit development, stress response, and other biological processes [24,25]. For example, COOLAIR (Cold induced antisense intragenic RNA) and COLDAIR (Cold assisted intronic noncoding RNA) were reported to inhibit the owering genes to achieve rapid owering transition in Arabidopsis [26,27]. In rice, a lncRNA named LDMAR (Long-day-speci c male-fertility associated RNA) was revealed as a participant in DNA methylation process, which resulted in male sterility [28]. TAS3 (trans-acting siRNA3), which was considered as a nitrogen (N)-responsive lncRNA in Arabidopsis, has been con rmed to promote lateral root development [29]. In addition to these biological processes, lncRNAs have also been reported to be involved in anthocyanin synthesis pathways. At present, the anthocyanin synthesis associated lncRNAs were identi ed in a few plant species like sea buckthorn, apple and strawberry [21,30,31]. For instance, during fruit ripening stage in sea buckthorn, LNC1 and LNC2 mediated the synthesis of anthocyanin by regulating the expression of SPL9 and MYB114 [31]. In apple, MLNC3.2 and MLNC4.6 promoted the expression of SPL2-like and SPL33 during light-induced anthocyanin synthesis [30]. To date, the indepth researches on the lncRNAs involved in anthocyanin synthesis mostly focused on the aboveground organs. The functions and regulatory mechanism of lncRNAs during anthocyanin biosynthesis in underground organs are still obscure.
In this study, transcriptome sequencing between yellow eshed ('Jin-16') and purple eshed potato tubers ('Xisen-8') were carried out to excavate the potential lncRNAs and their target genes that associated with anthocyanin biosynthesis. Furthermore, the functional enrichment analyses were conducted on the target genes of differentially expressed lncRNAs (DE lncRNAs) to reveal the potential biological processes and pathways in which lncRNAs involved. A lncRNA-mRNA coexpression network was constructed based on their correlation to unravel the regulatory mechanism of lncRNAs during anthocyanin accumulation. These ndings replenish the molecular mechanism of non-coding RNAs in the regulation of anthocyanin synthesis, which will be of great guiding signi cance for the cultivation of new potato varieties with high added value.

Results
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 esh 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 esh in 'Jin-16' and 'Xisen-8' were detected, respectively. As shown in Fig. 1c, the relative anthocyanin content of tuber esh in 'Xisen-8' was signi cantly 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.
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 unidenti able 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 ltered 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 le 1. Finally, 6371 non-coding transcripts corresponding 4376 lncRNAs were identi ed as the candidate lncRNA datasets based on the analyses of CPC2, PFAM and CNCI softwares (Additional le 2). According to the location relationship between lncRNAs and known mRNAs, lncRNAs could be classi ed 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 identi ed (Fig. 2c).

Construction of expression pro le and Chromosomal localization of differentially expressed lncRNAs in potato tubers
After difference signi cance 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 le 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).
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 le 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 identi ed lncRNAs in potato tubers were presented in Additional le 4 and Additional le 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 signi cantly enriched GO terms were shown in Fig. 4 and Additional le 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).
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 identi ed to participate in 121 and 119 metabolic pathways, respectively (Additional le 7). According to enrichment signi cance, 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 speci c 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 signi cantly upregulated, while the expression levels of CHI and ANS had no signi cantly 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 signi cantly upregulated 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. Analysis of lncRNA-mRNA co-expression network and identi cation 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 le 8). The intricate interaction between lncRNAs and their coexpressed 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 in uenced 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 speci c primer sequences of the lncRNAs and genes were shown in Additional le 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 (R 2 =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.

Discussion
The construction of expression pro le of 1421 DE lncRNAs provides new insights into the regulation of anthocyanin biosynthesis in potato tubers Anthocyanin is a kind of water-soluble pigment with strong antioxidant capacity in plants, contributing to a variety of vibrant colors in different organs and promoting the ornamental and nutritional value of plants [30]. Purple eshed potatoes are favored by consumers because they not only contain the general nutrients in white or yellow eshed potatoes, but also enrich high level of anthocyanin content [5]. With the development of processed food and health products from purple eshed potatoes, the increase of anthocyanin content in tubers has become the aim of potato breeders and growers.
Investigating the regulatory mechanism of anthocyanin biosynthesis is of great guiding signi cance and practical application value for improving potato tuber quality and facilitating the breeding of potato germplasm resources.
As a class of non-coding RNAs, lncRNAs which were initially regarded as "transcriptional noise" [23], have been reported to play crucial roles in various biological processes such as plant development [34], fruit ripening [35,36], stress responses [37, 38] and anthocyanin biosynthesis [31]. Although thousands of lncRNAs related to anthocyanin synthesis were identi ed in fruits of several plant species [21,30,31], the lncRNAs regulating anthocyanin biosynthesis in potato tubers have not yet been reported.
In this study, two potato cultivars 'Jin-16' (the yellow eshed potato) and 'Xisen-8' (the purple eshed potato with high anthocyanin content) were used as the experimental materials to excavate the lncRNAs involved in anthocyanin biosynthesis (Fig. 1). Finally, 4376 lncRNAs were identi ed in potato tubers including 3461 lincRNAs, 705 anti-sense lncRNAs and 210 sense overlapping lncRNAs (Fig. 2c). LincRNAs accounted for the largest proportion of the total lncRNAs in potato tubers, which were similar to those identi ed in other plants [24,39], suggesting that the overwhelming majority of lncRNAs are located in intergenic regions. Totally, 1421 DE lncRNAs, including 735 up-regulated and 686 down-regulated lncRNAs, were identi ed between 'Xisen-8' and 'Jin-16' libraries by analysis of transcriptome sequencing ( Fig. 3a and 3c). LncRNAs in sea buckthorn and strawberry fruits have been reported to display expression speci city in tissue types and developmental stages [21,31]. In potato tubers, 361 DE lncRNAs were only expressed in 'Xisen-8'; while 325 DE lncRNAs were observed to be speci cally expressed in 'Jin-16', suggesting that the expression of lncRNAs also has evident varietal speci city (Fig. 3b). The result of chromosome localization showed that 1312 of the 1421 DE lncRNAs were distributed in CHR1 ~ CHR12, and 109 DE lncRNAs which might be located on the unanchored scaffolds were placed on CHR0 (Fig. 3d).
The global identi cation of DE lncRNAs between the yellow and purple eshed potatoes provided a novel perspective on the regulation mechanism of anthocyanin biosynthesis in underground organs.
The DE lncRNAs and their target genes participate in diverse biological processes to mediate anthocyanin biosynthesis in potato tubers LncRNAs perform their functions by binding to DNA, RNA or proteins via multiple different regulatory mechanisms [23]. LncRNAs can act as cis-acting or trans-acting factors [34,35], miRNA target mimicry [30], miRNA precursor[36]; and can be involved in histone modi cation [37] as well as DNA methylation [38]. Although the functions and regulation mechanisms of most lncRNAs have not been studied clearly, previous studies indicated that lncRNAs are likely to partake in diversi ed biological processes by regulating the expression of their target genes either in cis-or in trans-acting [30,40,41]. Therefore, the function of lncRNAs depends on the function of their corresponding target genes. In this research, the cis (co-located) and trans (co-expressed) target genes of DE lncRNAs between 'Jin-16' and 'Xisen-8' were identi ed and subjected to GO and KEGG pathway enrichment analyses to predict the functions of these DE lncRNAs.
For the cis-acting target genes, most of them were associated with "cellular metabolic process", "phosphorus metabolic process" and "phosphate-containing compound metabolic process" which were belonged to BP category, suggesting that these genes might be involved in phosphorus metabolism to promote anthocyanin synthesis (Fig. 4a). Numerous researches have demonstrated that the anabolism of anthocyanin in plants is usually affected by multifarious environmental factors such as temperature, light, water, nitrogen and phosphorus [30,42,43]. Phosphorus is one of the most important factors in uencing anthocyanin biosynthesis and Pi de ciency results in purple color in leaves and stems due to the anthocyanin accumulation in plants [44]. Therefore, these related DE lncRNAs and their corresponding target genes probably play a systematic role in regulating phosphorus balance and anthocyanin synthesis in plants. Anthocyanin is synthesized from phenylalanine through a series of enzymatic reactions in cytoplasm and then transported to vacuoles for storage [11]. However, the stability of anthocyanin is low at neutral pH environment and must be enhanced by glycosidylation and acylation [45]. The involvement of the genes that were signi cantly enriched in the terms of "hydrolase activity" and "transferase activity" may enhance the anthocyanin stabilization and contribute to the biosynthesis of different anthocyanins. In addition, some genes were clustered into the terms of "membrane-enclosed lumen" and "organelle lumen" in CC category, indicating that these genes might participate in repairing the cellular components (especially lumen regions) under the regulation of their corresponding lncRNAs to provide a suitable microenvironment for the synthesis and accumulation of anthocyanin. For the trans-acting target genes of DE lncRNAs, in addition to the genes involved in phosphorus metabolism, catalytic reactions and other basic metabolic processes as described above, numerous genes were also linked to "binding" like "small molecule binding", "nucleotide binding", "nucleoside phosphate binding" and so on (Fig. 4b). This result was in accordance with that in cineraria [46], and suggested that the trans-acting genes of DE lncRNAs were likely to mediate the synthesis of biomacromolecules, the generation of copolymers and the regulation of gene expression during anthocyanin synthesis in potato tubers. Notably, these trans-acting target genes were not classi ed into the terms of CC category among the top 40 signi cantly enriched GO terms, revealing that they may not be primarily involved in the synthesis of subcellular structures [47].
The KEGG enrichment analysis showed that the target genes of DE lncRNAs were signi cantly enriched in "DNA replication", "Homologous recombination" and "Nucleotide excision repair" (Fig. 5), suggesting that these genes and their corresponding DE lncRNAs might participate in cell proliferation and regulate the formation and development of potato tubers. Similar results have been found in the development of tomato owers and fruits [48]. Anthocyanin is a kind of most important secondary metabolites, and its synthesis is a quite complicated process which is based on various primary metabolic pathways, such as glucose metabolism, fatty acid metabolism, amino acid metabolism and so on [49]. It has been reported that starch degradation could supply abundant substrates for anthocyanin synthesis in tuberous roots of purple eshed sweet potato [50]. Therefore, the genes enriched in "Biosynthesis of secondary metabolites", "Starch and sucrose metabolism", "Fatty acid biosynthesis and degradation", "Glycolysis/Gluconeogenesis" and different amino acid metabolism like "Glycine, Serine and Threonine metabolism" are possibly to play necessary roles in anthocyanin synthesis.
In addition, the genes involved in "plant hormone signal transduction" may also take part in the regulation of anthocyanin metabolism (Fig. 5b). DELLA protein, a key factor in gibberellin signal transduction pathway, has been revealed to regulate anthocyanin synthesis [51]. These results help to predict the functions of DE lncRNAs and provide a basis for the further study of the regulation mechanism of lncRNAs in anthocyanin biosynthesis.
The DE lncRNAs contribute to anthocyanin biosynthesis by targeting the key enzyme genes and transcription factor genes in potato tubers Numerous genes encoding enzymes and transcription factors have been recognized to be involved in the anthocyanin metabolic process [52][53][54]. Enzyme genes such as CHS, CHI, F3H, DFR, ANS, ANT and UFGT play indispensable roles in synthesis, modi cation and accumulation of anthocyanins in different plants [11,55,56]. During anthocyanin biosynthesis, p-Coumaroyl-CoA was synthesized from phenylalanine via a series of enzymes and competed by CHS and HQT [57]. In 'Xisen-8' tubers, the expression levels of CHS and F3H were up-regulated; while HQT was down-regulated, enhancing the production of dihydrokaempferol which is the substrate for synthesizng pelargonidin and avonol[58] (Fig. 6). Some researchers have proposed that DFR could catalyze DHK (dihydrokaempferol), DHQ (dihydroquercetin) and DHM (Dihydromynicetin) to synthesize different anthocyanins [59]. . During this process, the expression of these anthocyanin-related enzyme genes would be regulated by multiple lncRNAs in different plant species [21,30]. TCONS_01039552, a lncRNA in sea buckthorn fruit, was observed to regulate the expression of F3H [31]. In strawberry fruit, lncRNAs TRINITY_DN48515_c0_g3_i1 and TRINITY_DN1328_c0_g1_i1 were predicted to positively and negatively correlated with CHI and CHS, respectively [21]. In potato tubers, multiple lncRNAs such as XLOC_046144, XLOC_038153, XLOC_047085, XLOC_060098, XLOC_070549 were also predicted to in uence the anthocyanin content by regulating the expression of enzyme genes displayed in Fig. 6 (Table 1)

. Besides, CCoAOMT and
AnAT, which have been reported to be involved in methylation and acylation during anthocyanin metabolism respectively[63, 64], were also predicted to be regulated by lncRNAs like XLOC060098 and XLOC038153 in potato tubers (Fig. 8).
Additionally, a number of transcription factors that mediate anthocyanin metabolism are often regulated by lncRNAs as well [30]. In this study, the transcription factor genes such as MYB, SPL, ARF, NAC and ERF were identi ed to be targeted by different lncRNAs through the analysis of co-expression network ( Fig. 7 and Fig. 8). MYB transcription factor usually interacts with bHLH and WD40 to form the ternary WBM complex which subsequently regulates the expression of related enzyme genes like F3'5'H, DFR and UFGT, thereby participating in anthocyanin synthesis [15,[65][66][67]. In sea buckthorn fruit, LNC2 was identi ed as the endogenous target mimic (eTM) of miR828a, thereby inducing MYB114 expression [31]. It has been found in rice that lncRNA TWISTED LEAF could affect OsMYB60 expression and leaf phenotype [5]. Similarly, in potato tubers, the changes of MYB expression under the guidance of the endogenous lncRNAs, such as XLOC_017372, XLOC_029213, and XLOC_038153, might in uence the anthocyanin accumulation and tuber colors (Additional le 8).
SPL, which was considered as the target of miR156, has been reported to negatively regulate anthocyanin biosynthesis under the action of lncRNAs in different plant species [11,16,68]. In banana, the lncRNA could modulate the expression of SPL2-like gene [68]. LNC1 was predicted to be a decoy for miR156a and reduce SPL9 expression, thus increasing anthocyanin content in sea buckhorn fruit [31]. In apple, MLNC3.2 and MLNC4.6 were recognized as eTMs for miRNA156a and facilitated the expression of the SPL2-like and SPL33 by preventing cleavage of these two SPL genes [30]. They also proposed that the regulatory mechanism by which lncRNAs act as eTMs for miR156 may be conserved and universal in plants, and these lncRNAs may in uence various aspects of plant development [30]. In 'Xisen-8', the down-regulated XLOC_100954 reduced the expression of SPL and promoted anthocyanin production (Fig. 7, Fig. 8). Therefore, XLOC_100954 might function as the eTM for miR156 in potato and prevent breakage of SPL by miR156. However, the speci c regulation mechanism of lncRNAs on SPL needs further validation in potato tubers.
ARF was regarded as a key regulator of auxin responsive signaling pathway and plays a crucial role in anthocyanin biosynthesis in plants [69]. MdARF13 was proved to suppress MdDFR expression to decrease anthocyanin accumulation [70]. Recently, 6 PeARFs and 5 lncRNAs were predicted to be targeted by pei-miR160a in populus [42] . But the regulatory relationship between lncRNAs and ARFs has not been clari ed. In potato tubers, ARF8 was identi ed as the coexpressed target of XLOC_100954 (Fig. 7). The low expression of XLOC_100954 affected ARF8 expression and thus increased the anthocyanin content in purple eshed potato tubers (Fig. 8). Other anthocyanin-related transcription factor genes like NAC2 and ERF4 were also predicted to be regulated by XLOC_060098 and XLOC_017372, respectively. Although several lncRNAs were predicted to regulate the expression of NAC and EFR in other plants [71], the speci c action mechanism remains unclear.
The regulatory relationship between speci c lncRNAs and mRNAs is quite complicated due to the involvement of many other factors like miRNAs and other genes. In this study, DE lncRNAs and their potential target genes were identi ed and predicted by comparing transcriptome data from purple and yellow eshed potato tubers. The target genes co-expressed with XLOC_060098 and XLOC_017372 were the most, indicating that XLOC_060098 and XLOC_017372 may be involved in the most biological processes and play crucial roles during anthocyanin metabolism. These results provide a basis for further investigation of the function of potato lncRNAs in anthocyanin synthesis. The analysis of the regulation mechanism of anthocyanin biosynthesis in potato tubers enriches the existing function mechanism of anthocyanin biosynthesis in different plant organs especially the underground organs.

Conclusions
In this study, fresh tubers of purple eshed potato ('Xisen-8') and yellow eshed potato ('Jin-16') were used as experimental materials to excavate lncRNAs related to potato anthocyanin biosynthesis through transcriptome sequencing. The cis-and trans-acting target genes of the DE lncRNAs were predicted via bioinformatics analyses. Totally, 4376 lncRNAs were identi ed, among which 1421 were signi cantly differentially expressed with 735 up-regulated and 686 down-regulated lncRNAs in 'Xisen-8'. Chromosomal mapping analysis showed that these DE lncRNAs were unevenly distributed on potato chromosomes: the number of lncRNAs on chromosome 1 was the highest, and the number of lncRNAs on chromosome 2 was the lowest. The target genes of DE lncRNAs could be clustered into lots of GO types and enriched into multiple metabolic pathways, thus refecting the functional diversity of their corresponding lncRNAs. The co-expression network and the expression patterns of DE lncRNAs and their target genes related to anthocyanin biosynthesis in potato tubers were analyzed to reveal their potential regulatory relationship. Results showed that most anthocyanin-associated enzyme genes and transcription factor genes were predicted to be regulated by lncRNAs such as XLOC_060098 and XLOC_017372, indicating that these lncRNAs might play crucial roles in anthocyanin metabolism. The results of this study lay the foundation for further elucidation of regulatory mechanism of lncRNAs in anthocyanin biosynthesis and provide new insights for quality improvement strategies of potato and other crops.

Determination of anthocyanin contents in potato tubers
Three potato tubers with similar sizes (5-6 cm in length) were selected from 'Jin-16' and 'Xisen-8', respectively. Anthocyanin was extracted according to the method used by Wang et al. [70]. The potato esh from each tuber was ground into powder and then exposed to HCl-methanol solution (1:99 by volume) at 4 ℃ for 6-8 h under darkness until the tissues were completely decolorized. After centrifuging at 12000 rpm for 10 min, the absorbance values of supernatants were determined at 530 nm using a UV-2450 spectrophotometer (Shimadzu, Kyoto, Japan). Each sample had three replicates to ensure the results reliable.
Total RNA extraction, library preparation and lncRNA sequencing Total RNA was isolated from 0.5 g potato esh samples using the Quick RNA Isolation Kit (Huayueyang, Beijing, China).
Subsequently, electrophoresis was performed with 1% agarose gel to monitor the presence of RNA degradation and DNA contamination. The purity and concentration of RNA samples were measured using Nanodrop 1000 spectrophotometer (Thermo Scienti c). After integrity testing by Agilent 2100 BioAnalyzer (Agilent Technologies), the total RNA samples were used for the construction of lncRNA libraries and validation of deep sequencing results.
Strand-speci c library was constructed using the deoxy-UTP (dUTP) strand-marking method [72]. Three biological replicates were set for each potato cultivar. Therefore, 6 strand-speci c libraries were constructed, including Jin-16_1, Jin-16_2, Jin-16_3, Xisen-8_1, Xisen-8_2 and Xisen-8_3. After quantitative and qualitative determination of all libraries, lncRNA deep sequencing were carried out on an Illumina novaseq 6000 platform by Novogene Bioinformatics Technology Co. Ltd. (Beijing, China). The obtained raw reads were processed by getting rid of the low quality reads, the reads with sequencing adapters and poly-N sequences. The clean reads were acquired and aligned to a potato reference genome (S. tuberosum assembly 4.03) using Hisat2 software [32]. The mapped reads were spliced and assembled into transcripts using Stringtie [73] and Cuffmerge software.

Identi cation of differentially expressed lncRNAs
The identi cation of lncRNAs was based on their structural characteristics and non-coding functional characteristics. The transcripts with less than two exons and transcripts less than 200 nt in length were removed. CPC2, PFAM and CNCI tools were employed to estimate the coding potential of the spliced transcripts. Based on the comprehensive analysis results of the three softwares, the transcripts without coding potential were regarded as candidate lncRNAs. The differential expression analysis of identi ed lncRNAs was performed using DESeq [18] and edgeR [19]. The lncRNAs with q value less than 0.05 were designated as differentially expressed lncRNAs. Prediction and functional enrichment of lncRNA target genes LncRNAs regulate target genes through a variety of mechanisms, and the most common ways of acting on downstream target genes are co-location and co-expression regulation. The coding genes which located in 100 kb upstream and downstream of lncRNAs were searched and classi ed into co-located target genes. The co-expressed target genes were predicted by correlation analysis between lncRNAs and mRNAs with the correlation coe cient above 0.95.
To predict the main functions of lncRNAs, the functional enrichment analyses were performed on the co-located and coexpressed target genes, respectively. The sequences of the target genes were blasted to the Gene Ontology (GO) database to obtain their GO annotation information. Then GO enrichment analysis was implemented by GOseq software [74]. The GO terms with p-adjust < 0.05 were considered signi cantly enriched.
The gene sequences were also mapped to the KEGG database to detect the main pathways in which they participate.
KOBAS (2.0) software [75] was used to perform the pathway enrichment analysis of target genes of DE lncRNAs. The pathways with p value < 0.05 were regarded as the signi cant enriched pathways.
Construction of co-expression network between DE lncRNAs and their co-expressed target genes According to the predicted co-expression relationship between DE lncRNAs and their co-expressed mRNAs related to anthocyanin biosynthesis, lncRNA-mRNA pairs with correlation coe cient greater than 0.98 were screened out. The

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The datasets generated and/or analyzed for this work were deposited in the NCBI Sequence Read Archive under the Bioproject accession PRJNA729884, available from https://dataview.ncbi.nlm.nih.gov/object/PRJNA729884? reviewer=ntlkjmravag9c9ousg57ps9k86.     Gene Ontology classi cation of the co-located target genes (a) and co-expressed target genes (b) of DE lncRNAs. GO functions were presented in X-axis; the number and the percentage of annotated target genes in each GO term were presented in left and right Y-axis, respectively.

Figure 7
Co-expression network between DE lncRNAs and their co-expressed targets related to anthocyanin synthesis. Nodes represented DE lncRNAs and mRNAs in potato, edges indicated pairwise correlation. LncRNAs and mRNAs were represented as squares and circles, respectively. The different colors of the nodes indicated the numbers of lncRNA or mRNA involved in the interaction: the darker the color, the more interactions were involved. The colors were in order from dark to light: purple, dark green, light green and yellow.