Genome-wide identification of potato lincRNAsTo identify and understand the function of lincRNAs during the progression of ZC disease, we performed RNA-seq of healthy (CT; control) and infested potato plants with psyllid vector carrying CLso [CLso(+)] and psyllids without CLso [CLso(−)] temporally at 7, 14, and 21 days after infection (DAI). In total, ~298 million paired-end reads were generated from the temporal RNA-seq datasets (nine conditions having three biological replicates in each condition), covered with ~11.4 million paired-end reads in each biological replicate. The raw sequence data were subjected to quality filtering and mapped to the reference potato genome (Solanum tuberosum (S. tuberosum) v4.04). Of these, ~99% of the reads were high quality, and 54–61% mapped uniquely to the reference genome27 (Supplementary Table 1). Transcripts were assembled using the Read mapping and transcript assembly workflow28, and lincRNAs were identified using Evolinc29. After identifying putative lincRNAs, those expressed below 0.1 transcripts per million were considered artifacts and removed from further analysis.We identified a total of ~4397 lincRNAs from the nine samples of potatoes representing healthy (CT), CLso(+), and CLso(−) conditions at 7, 14, and 21 DAI (Supplementary Table 1). Previous studies showed that lincRNAs are mostly species-specific9,30,31. For instance, only ~12% of human lincRNAs show sequence conservation with non-human species32. Consistent with this, when we compared potato lincRNAs with those of Arabidopsis and rice (downloaded from the CANTATAdb v2.0 database), we observed poor inter-species conservation (<1%) of the lncRNAs. In contrast, ~35% of potato lincRNAs were conserved with the lncRNAs of other Solanaceous plants (Fig. 1a; Supplementary Data 1a). Moreover, 4379 lincRNAs identified in this study were substantially higher than the previously reported lincRNAs under temporal infection with P. carotovorum subsp. brasiliense in the tolerant and susceptible potato cultivars17. On the contrary, about 3175 lncRNAs were detected from the apical buds of potato tubers21. Therefore, this study expanded the potato lincRNAs repertoire to include novel ones perturbed during biotic stress response.Fig. 1: Identification and characteristics of potato lincRNAs.a A Venn diagram compares potato lncRNAs available in the CANTATAdb v2.0 database with those detected in the current study. b The distribution of lincRNAs in the 12 potato chromosomes is shown in a pie chart. The fraction (%) of the highest and least number of lincRNAs is depicted. c Lengthwise distribution of the lincRNAs and comparison with protein-coding genes (PCGs) is shown via kernel density plot.Genomic attributes of potato lincRNAsThe identified lincRNAs were distributed across all chromosomes (Chr) of the potato, with the most found on Chr 1 (12.66%) and the least on Chr 5 (5.93%) (Fig. 1b; Supplementary Data 1b). The transcript length of the lincRNAs ranged from 201 to 5095 bp with an average size of 491 bp, and ~95% of them were shorter than 1 kb. The length distribution of lincRNAs and PCGs revealed that the length of lincRNAs (average length 491 bp) was significantly shorter than that of PCGs (average length 1415 bp) (p value < 0.05, Welch’s two-sided t-test) (Fig. 1c; Supplementary Data 1c). Likewise, we analyzed the structural organization of the lincRNAs. In general, lincRNAs were predominantly mono- or di-exonic, whereas >50% of PCGs contained three or more exons (Fig. 2a; Supplementary Data 2). Consistent with other plant species, including Arabidopsis maize and soybean (Glycine max)3,7,8,11, the proportion of potato lincRNAs containing a single exon (~58%) was higher than that of lincRNAs containing multiple exons.Fig. 2: Distribution of exons and orientation of potato lincRNAs.a The distribution of exons in lincRNAs and comparison with PCGs is shown in bar plot. b Types of lincRNAs based on orientation concerning their proximal PCGs are graphically depicted.Previous studies suggested that lincRNAs can influence the expression of their proximal PCGs10,17. We analyzed PCGs located in the proximal (≤50 kb) regions in the up- or downstream regions of the lincRNAs. The average distance of the lincRNAs from their nearest PCGs was ~25 kb away. About 89% of the lincRNAs were located within 50 kb regions of their nearest PCGs, and therefore, we assigned these genes as proximal PCGs hereafter. Further, we analyzed the orientation of the lincRNAs relative to their proximal PCGs. Almost half (49–51%) of the lincRNAs were found to be oriented either in the same or opposite direction respective to proximal PCGs. Among the ~51% of lincRNAs oriented in the opposite direction, convergent type (~36%) represented a major fraction as compared to divergently (~15%) oriented lincRNAs (Fig. 2b). A previous study also showed a higher representation of convergently orientated lincRNAs with respect to their associated proximal PCGs33.Functional relevance of potato lincRNAs as plausible miRNA decoysSome lincRNAs could act as miRNA decoys by mimicking the target mRNA sequences, and consequently, the target mRNAs can be protected from degradation34,35,36. For instance, in Arabidopsis, the lincRNA Induced by Phosphate Starvation1 (IPS1) serves as a decoy to ath-miR399. It interferes with the binding of ath-miR399 to its primary target PHO2 transcript and rescue from post-transcriptional silencing in governing phosphate homeostasis in Arabidopsis37. However, lincRNAs that may act as miRNA decoys in potatoes have not been explored. The 4397 potato lincRNAs identified in this study were compared for sequence homology against the known potato mature miRNA sequences in the CANTATAdb v2.0 database. Six lincRNAs (lincRNA.ID.00006654, lincRNA.ID.00002569, lincRNA.ID.00018937, lincRNA.ID.00004487, lincRNA.ID.00010417, and lincRNA.ID.00012834) were predicted to harbor complementary sequences to known miRNAs and could serve as potential miRNA decoys (target mimics). Among them, lincRNA.ID.00006654 and five PCGs harbored binding sites of stu-miR482a-5p. Likewise, lincRNA.ID.00012834 along with three PCGs and lincRNA.ID.00002569, lincRNA.ID.00018937, lincRNA.ID.00004487, and lincRNA.ID.00010417 each with two PCGs harbored binding sites of same miRNAs (Supplementary Fig. 1a; Supplementary Data 3,  8). The PCGs encoded ATP binding protein, expansin, DNA repair protein, multidrug resistance, carboxylic ester hydrolase, a ubiquitin ligase, methyladenine DNA glycosylase, pentatricopeptide repeat protein, and potassium transporter.Of the six, two lncRNAs (lincRNA.ID.00002569 and lincRNA.ID.00004487) were downregulated in CLso(+) infected tissues at 21 DAI based on the transcriptome analysis (Supplementary Fig. 1b; Supplementary Data 3). However, only one lincRNA.ID.00002569 showed significant downregulation by RT-qPCR (Supplementary Fig. 1c). This lincRNA had complementarity to a potato stu-miR8005a miRNA and similarity to the target PCGs: PGSC0003DMT400008676 (encoding Rad50 ATPase) and PGSC0003DMT400011415 (encoding Multidrug resistance pump protein). The RT-qPCR results also showed that the downregulation of lincRNA.ID.00002569 in CLso(+) infection paralleled the downregulation of PGSC0003DMT400011415 but not PGSC0003DMT400008676 (Supplementary Fig. 1; Supplementary Data 3), suggesting lincRNA.ID.00002569 could act as a potential miRNA decoy of PGSC0003DMT400011415.Differentially expressed potato lincRNAs in response to CLso and psyllidLincRNAs are important in regulating diverse biological processes via complex regulatory networks with PCGs. To examine the role of lincRNAs under pathogen infection, we identified sets of differentially expressed lincRNAs under conditions of psyllid challenges without CLso [CLso(−)/CT], psyllid carrying CLso [CLso(+)/CT] and pathogen-specific [CLso(+)/CT]/[CLso(−)/CT] temporally at 7, 14, and 21 DAI. A total of 775 unique lincRNAs were found to be differentially expressed in any of the 9 different comparisons (Supplementary Data 4). A substantially higher number of lincRNAs were found to be differentially expressed at 21 DAI preferentially under [CLso(+)/CT] followed by [CLso(+)/CT]/[CLso(−)/CT] condition (Fig. 3a; Supplementary Data 4). Furthermore, we performed principal component analysis (PCA) of the nine sets of differentially expressed lincRNAs. The first component (PC1) variation was up to 88.4%. We ruled out the extent of variation due to outliers and lack of reproducibility among the biological replicates since we used the differentially expressed values in the log2 scale. Moreover, variation due to time of sampling is unlikely as we observed a distinct and diverse differential expression profile of lincRNAs specifically at 21 DAI (Fig. 3b). The results suggest a dynamic and plausible role of lincRNAs in response to CLso and psyllid infection at the terminal stage of ZC disease in potato.Fig. 3: Differential expression profiles of lincRNAs modulated during Candidatus Liberibacter spp. (CLso) infection.a Differential expression profile of lincRNAs under CLso(−)/CT, CLso(+)/CT, and CLso(+)/CLso(−) conditions at 7 days after infection (DAI), 14 DAI, and 21 DAI is shown via heatmap. The scale represents differential expression in log2 fold change. b The similarity/difference among the nine sets of differentially expressed lincRNAs described in (a) is shown via a PCA plot. Each data point represents the differential expression status of individual lincRNA.Further, we specifically determined up- and downregulated lincRNAs among the three sets of comparisons [CLso(−)/CT, CLso(+)/CT, and CLso(+)/CLso(−)] within each stage of infection. For instance, at 7 DAI, 8 and 11 lincRNAs were specifically up- and downregulated under CLso(−)/CT condition, suggesting a psyllid-specific response. Likewise, four and six lincRNAs exhibited up- and downregulation under CLso(+)/CT condition at 7 DAI, suggesting either a psyllid or CLso-specific response. However, only one lincRNA showed upregulation under the CLso(+)/CLso(−) condition at 7 DAI (Fig. 4a; Supplementary Data 4), suggesting a CLso-specific response. Similarly, quite a few differentially expressed lincRNAs exhibiting specific expression in the three sets of comparisons were detected at 14 DAI, too (Fig. 4b; Supplementary Data 4). At 21 DAI, the number of specifically up- and downregulated lincRNAs was substantially high under CLso(+)/CT [up: 133 and down: 161] followed by CLso(+)/CLso(−) condition [up: 20 and down: 34] (Fig. 4c; Supplementary Data 4). These results suggest a dynamic regulation of lincRNAs at the subsequent stages of ZC disease in response to CLso infection preferentially at 21 DAI. Furthermore, we examined the PCGs located in the proximal regions of the specifically up- and downregulated lincRNAs under CLso(+)/CT condition at 21 DAI as they constitute a major fraction (294) of the differentially expressed lncRNAs. Intriguingly, PCGs encoding TFs (HBP and HAT), signaling transduction components (auxin/ethylene signaling genes and serine/threonine protein kinase), and those involved in stress responses, including, cytochrome P450 monooxygenase CYP736B, disease resistance protein BS2, glutathione S-transferase (GST), lipoxygenase, resistance protein PSH-RGH6 and SAUR family protein were detected at the proximal regions of the differentially expressed lncRNAs. The dynamic expression of the lincRNAs may influence the regulation of these PCGs and could impact disease tolerance.Fig. 4: Comparison of differentially expressed lincRNAs within each landmark stage.Several common and unique lncRNAs exhibiting differential expression under CLso(−)/CT, CLso(+)/CT, and CLso(+)/CLso(−) conditions at 7 DAI (a), 14 DAI (b), and 21 DAI (c) are depicted via Venn diagrams. Labels on top (black), bottom left (red), and bottom right (blue) indicate the number of differentially expressed, upregulated, and downregulated lncRNAs, respectively.Transcriptional regulatory networks of potato lincRNAs and protein-coding genesThe interactions of lincRNAs and PCGs are complex and poorly understood. However, a few studies revealed that lincRNAs can also influence the expression of their proximal PCGs4,10. We examined to identify sets of co-up and co-downregulated lincRNAs–PCGs pairs located proximal to each other under CLso(−)/CT, CLso(+)/CT, and CLso(−)/CLso(+) conditions at 7, 14, and 21 DAI. We did not observe co-differentially expressed lincRNAs and proximal PCG pairs at 7 and 14 DAI. However, at least 25 co-up and 21 co-downregulated lincRNA and proximal PCG pairs were detected under CLso(+)/CT condition at 21 DAI (Supplementary Data 9). Pearson’s correlation for the set of co-upregulated (r = 0.77, p < 0.01) lincRNA–PCG pairs was higher than the set of co-downregulated pairs (r = 0.31). Among the 46 co-regulated lincRNA–PCG pairs, a majority of them were oriented in opposite directions (29), which were further categorized into convergent (16) and divergent (13) types. Conversely, the remaining 17 lincRNA–PCG pairs were oriented in the same direction. Many of the co-differentially expressed lncRNAs, along with their proximal PCGs, were involved in biotic stress responses38,39,40,41,42,43 such as wound-inducible carboxypeptidase, SAUR, GST, F-box, and chitinases (Supplementary Data 9).Further, to examine the co-expression networks among transcripts, we performed weighted gene co-expression network analysis (WGCNA)44 to determine sets of co-upregulated transcripts under CLso(−)/CT, CLso(+)/CT, and CLso(+)/CLso(−) conditions at 7, 14, and 21 DAI. At least four distinct modules (turquoise, blue, brown, and gray) representing distinct sets of differentially co-expressed transcripts were detected (Fig. 5a, b). The turquoise module was the largest among the different modules and comprised 331 transcripts (lincRNAs and PCGs) (Supplementary Data 5). The turquoise module was significantly associated with the CLso(+) treatment at 21 DAI (r = 0.93 and p = 1e − 12) (Fig. 5b). Further, we examined the functional relevance of the co-expressed transcripts belonging to the turquoise module by determining the enrichment of gene family and gene ontology (GO) terms using the corresponding PCGs as input. Functional terms related to acyl lipid metabolism superfamily, APETALA2 (AP2)/ethylene-responsive element binding proteins (EREBP), nucleoside transporter, GST, protein kinase and small auxin upregulated RNA (SAUR) processes were significantly enriched (Supplementary Data 10). Previous studies have shown that genes encoding AP2/EREBP, GST, and SAUR play crucial roles in plant defense against pathogen invasion41,42,43. In addition, GO terms involved in cell wall modification (GO:0009827), defense response to bacteria (GO:1900424), ion transport (GO:0006811), DNA binding (GO:0003677), glutathione peroxidase activity (GO:0004602), cell wall (GO:0005618), and plasma membrane (GO:0005886) were found to be significantly enriched (Supplementary Data 10). The results suggest the possible contribution of lincRNAs in governing different biological processes involved in biotic stress responses.Fig. 5: Identification and interaction of co-expressed transcripts.a A hierarchical clustering of the transcripts exhibiting a co-expression network based on topological overlap dissimilarity is shown via dendrogram. b Different sets (modules) of co-expressed transcripts in CT, CLso(−), and CLso(+) conditions at 7 DAI, 14 DAI, and 21 DAI are shown via heatmap. Labels in each module indicate Pearson’s correlation (R) and significance level (p value). The scale on the right side indicates Pearson’s correlation (R). c A co-expression network of the topmost 24 hub transcripts belonging to the turquoise module is shown. Light-green triangles indicate lncRNAs specific to CLso. Likewise, pink and blue ovals indicate lncRNAs and PCGs, respectively.Further, we examined regulatory networks among the transcripts of either lincRNAs or PCGs within the turquoise module to identify hub transcripts. A total of 24 transcripts exhibited the most interactions among themselves, and we considered them hub transcripts. Among them, 11 transcripts belonged to lincRNAs, including seven lincRNAs specific to CLso (Fig. 5c). The remaining 13 transcripts were PCGs, including genes encoding expansin, GST, and SAUR. The role of expansin and SAUR genes was implicated in defense responses against different types of biotic and abiotic stresses43,45. Overexpression of expansin genes led to increased tolerance in response to varying abiotic stresses, such as heat, drought, and salinity45. In plants, the accumulation of reactive oxygen species and programmed cell death processes has been implicated in response to pathogen invasion. Glutathione peroxidase activity encoded by GSTs inhibits the spread of cell death in infected plants46. Likewise, the function of SAUR genes in cell wall biosynthesis has been implicated in biotic stress responses43,47. The results suggest an intricate network between coding and non-coding transcripts preferentially at the later stage (21 DAI) of ZC disease progression.Overexpression of lincRNAs confers potato tolerance to CLsoNext, we sought to characterize if lincRNAs can mediate resistance to CLso infection in potatoes. Four lincRNAs were selected with criteria of (1) co-upregulated along with their proximal PCGs, (4) belonging to the hub transcripts, and (3) progressively induced (7 and 21 DAI) during the progression of ZC disease. The respective lincRNAs were cloned in an overexpression vector under the control of a double-enhanced CaMV35S promoter (DE35S-P). Moreover, the vector harbored a GFP reporter (Fig. 6a). The constructs were used for Rhizobium rhizogenes (R. rhizogenes)-mediated potato hairy root transformation48. After 28 days post-transformation, transgenic hairy roots were analyzed to confirm the overexpression of the four lincRNAs via RT-qPCR and GFP marker fluorescence (Fig. 6b, c; Supplementary Data 6). The transgenic hairy roots showed significantly higher expression of the lincRNAs than control plants (Fig. 6c; Supplementary Data 6).Fig. 6: Overexpression of lincRNAs conferred tolerance to Candidatus Liberibacter spp.a Schematic diagram showing plasmid vector employed for overexpression (OE) of the selected lincRNAs. b Photographs showing the phenotype of the transformed potato plants at 28 d post-transformation and GFP fluorescence. Scale bars indicate 1 cm (left) and 10 µm (right). c The relative expression level of the lincRNAs between the CT and OE plants estimated via RT-qPCR is shown in the bar plot. d Likewise, the relative expression level of the proximal PCGs located either in up- and downstream regions of the lincRNAs between CT and OE plants estimated via RT-qPCR is shown in bar plots. e Relative Candidatus Liberibacter spp. (CLso) titer between CT and the OE plants estimated via qPCR is shown via bar plot. Asterisks represent statistically significant (p value ≤ 0.05) differences compared to control (CT) as estimated by Student’s t-test. The error bars represent the standard error among the three biological replicates.To examine the influence of overexpression of the lincRNAs on their proximal PCGs, we selected three PCGs, each located in the upstream and downstream regions of the lincRNA. The expression level of these proximal PCGs was compared between the control and overexpression lines via RT-qPCR. We observed a predominantly negative correlation between higher expression of the lincRNAs and their proximal PCGs expression. For instance, overexpressing lincRNA.ID.00017498 exhibited downregulation of all the six proximal PCGs compared to their controls. Interestingly, RNAi-mediated knockdown of lincRNA.ID.00017498 showed converse effects on the expression dynamics of the proximal PCGs, implicating a role of the lincRNAs in regulating PCG gene expression (Supplementary Fig. 2; Supplementary Data 7). Likewise, overexpression of lincRNA.ID.00013924 reduced the expression of the most proximal PCGs except one located in the upstream region. Conversely, the influence of overexpressing lincRNA.ID.00017902 and lincRNA.ID.00013620 represented both higher and lower expression of their proximal PCGs (Fig. 6d; Supplementary Data 6). The results suggest that lincRNAs could govern the expression level of their proximal PCGs both negatively and positively. Previous studies have also shown that lincRNAs can influence the expression of their proximal PCGs4,10. The functional annotations of the proximal PCGs, with altered expression in the lincRNA-overexpressed hairy roots, included genes encoding proteinases, ubiquitin-protein ligase components, B2 proteins, and ribosomal proteins. Many of these genes could impact plant biotic stress responses49,50,51. Next, we determined whether overexpression of candidate lincRNAs impacted potato tolerance to CLso. Overexpression of three of the four lincRNA exhibited significantly lower CLso titers in transgenic hairy roots compared to control roots transformed with an empty vector. The amplitude of CLso titer reduction was about sixfold in lincRNA.ID.0007902 and lincRNA.ID.0007498 overexpressing hairy roots, while about halffold in the lincRNA.ID.0003924 hairy roots (Fig. 6e; Supplementary Data 6). The results suggest that these three lincRNAs could also impact CLso accumulation.
