Transcriptome-Wide Identification and Characterization of Potato Circular RNAs in Response to Pectobacterium carotovorum Subspecies brasiliense Infection

Little information about the roles of circular RNAs (circRNAs) during potato-Pectobacterium carotovorum subsp. brasiliense (Pcb) interaction is currently available. In this study, we conducted the systematic identification of circRNAs from time series samples of potato cultivars Valor (susceptible) and BP1 (disease tolerant) infected by Pcb. A total of 2098 circRNAs were detected and about half (931, 44.38%) were intergenic circRNAs. And differential expression analysis detected 429 significantly regulated circRNAs. circRNAs play roles by regulating parental genes and sponging miRNAs. Gene Ontology (GO) enrichment of parental genes and miRNAs targeted mRNAs revealed that these differentially expressed (DE) circRNAs were involved in defense response (GO:0006952), cell wall (GO:0005199), ADP binding (GO:0043531), phosphorylation (GO:0016310), and kinase activity (GO:0016301), suggesting the roles of circRNAs in regulating potato immune response. Furthermore, weighted gene co-expression network analysis (WGCNA) found that circRNAs were closely related with coding-genes and long intergenic noncoding RNAs (lincRNAs). And together they were cultivar-specifically regulated to strengthen immune response of potato to Pcb infection, implying the roles of circRNAs in reprogramming disease responsive transcriptome. Our results will provide new insights into the potato-Pcb interaction and may lead to novel disease control strategy in the future.


Introduction
Circular RNAs (circRNAs) have recently become the intensive study object in a wide variety of organisms [1]. Different from the traditional linear RNAs, circRNAs lack 5' and 3' ends and usually are generated by head-to-tail backsplicing events [2]. Its observation was firstly described as rare events, nevertheless, in recent years, increasing evidence suggests that circRNAs are widely existent in eukaryotic cells and commonly produced by thousands of genes [3].

Genome-Wide Identification of circRNAs in Potato
Until now, circRNAs have been identified in multiple monocot and dicot plants [13]. However, our current understanding of the characteristic and biological roles of circRNAs in the potato-pathogen interaction remains limited. In this study, 13.8 billion 90 bp paired-end raw reads from two potato cultivars were used for circRNAs identification. Low quality reads were filtered and clean reads were mapped to reference genome by Burrows-Wheeler Aligner (BWA) (v0.7.15, mem - T 19), resulting in about 10.6 billion reads being uniquely mapped. The Sequence Alignment/Map (SAM) files of alignment were then scanned by circRNA identification tool CIRI2, which led to the identification of 76,882 junction reads. And 2098 circRNAs were finally detected in these thirty samples supported by at least two unique back-spliced reads, with 1404 in cultivar Valor and 1337 in cultivar BP1 (Tables S1  and S2). Among these identified circRNAs, 761 were specifically expressed in Valor, and 694 in BP1, respectively, suggesting that circRNAs are cultivar specific ( Figure 1A,B).

CircRNAs Are Differentially Expressed in Response to Pcb Infection
These two potato cultivars, Valor (susceptible) and BP1 (disease tolerant), enabled us to identity the changes of circRNAs expression during infection and trace the underlying disease resistancerelated circRNAs. Wilcoxon rank-sum test indicated that circRNAs are significantly down-regulated by Pcb infection in susceptible cultivar Valor (CK vs. 24 hpi, p = 0.04; CK vs. 72 hpi, p = 0.02; CK means control and hpi means hours post inoculation), but significantly up-regulated by Pcb infection in disease tolerant cultivar BP1 (CK vs. 12 hpi, p = 0.02; CK vs. 24 hpi, p = 0.05; CK vs. 72 hpi, p = 0.001), suggesting that they have different regulation patterns in two cultivars and likely play potential roles in immune response (Figure 2A). Among the 2098 circRNAs, 429 circRNAs had significant

CircRNAs Are Differentially Expressed in Response to Pcb Infection
These two potato cultivars, Valor (susceptible) and BP1 (disease tolerant), enabled us to identity the changes of circRNAs expression during infection and trace the underlying disease resistance-related circRNAs. Wilcoxon rank-sum test indicated that circRNAs are significantly down-regulated by Pcb infection in susceptible cultivar Valor (CK vs. 24 hpi, p = 0.04; CK vs. 72 hpi, p = 0.02; CK means control and hpi means hours post inoculation), but significantly up-regulated by Pcb infection in disease tolerant cultivar BP1 (CK vs. 12 hpi, p = 0.02; CK vs. 24 hpi, p = 0.05; CK vs. 72 hpi, p = 0.001), suggesting that they have different regulation patterns in two cultivars and likely play potential roles in immune response (Figure 2A). Among the 2098 circRNAs, 429 circRNAs had significant differential expression levels in response to Pcb infection, including 275 in BP1 and 287 in Valor, and 428 between BP1 and Valor ( Figure 2B,  Figure 2B, Table S4).

CircRNAs Act as Putative Regulators of Parental Genes
It was reported that circRNAs play important roles in controlling transcription by cis-regulating their parental genes [13]. In order to reveal the functions of DE circRNAs between susceptible and disease tolerant potato cultivars, the parental genes of significantly regulated circRNAs were predicted and annotated. GO enrichment analysis showed that these GO terms related to: (1) defense response (GO:0006952); (2) substance metabolism, including pectin catabolic process (GO:0045490), structural constituent of cell wall (GO:0005199), xanthophyll metabolic process (GO:0016122); and (3) energy metabolism, including chloroplast photosystem II (GO:0030095), ADP binding (GO:0043531), were specifically enriched. The enrichment result suggested the involvement of circRNAs in controlling potato responding to pathogen infection ( Figure 3A, Table S5). These results are in accordance with those of previous studies, which suggested that plant resistance is a complex regulating network that involves in substance and energy metabolism, and defense process [21,22].

CircRNAs Act as Putative Regulators of Parental Genes
It was reported that circRNAs play important roles in controlling transcription by cis-regulating their parental genes [13]. In order to reveal the functions of DE circRNAs between susceptible and disease tolerant potato cultivars, the parental genes of significantly regulated circRNAs were predicted and annotated. GO enrichment analysis showed that these GO terms related to: (1) defense response (GO:0006952); (2) substance metabolism, including pectin catabolic process (GO:0045490), structural constituent of cell wall (GO:0005199), xanthophyll metabolic process (GO:0016122); and (3) energy metabolism, including chloroplast photosystem II (GO:0030095), ADP binding (GO:0043531), were specifically enriched. The enrichment result suggested the involvement of circRNAs in controlling potato responding to pathogen infection ( Figure 3A, Table S5). These results are in accordance with those of previous studies, which suggested that plant resistance is a complex regulating network that involves in substance and energy metabolism, and defense process [21,22]. Cell wall metabolism is an important aspect in the basal disease resistance process [23]. In our study, we found that lots of circRNAs were produced by protein-coding genes that are involved in cell wall metabolism. For example, PGSC0003DMG400021603, annotated as extensin, can generate 165 circRNAs. And PGSC0003DMG400009783, annotated as vegetative cell wall protein, can give rise to 46 circRNA (Table S3). The high abundance of circRNAs derived from the cell wall related genes implies its biological importance. Consequently, it is plausible to assume that circRNAs play important roles during the interaction between potato and pathogen Pcb. Moreover, pectin catabolic process (GO:0045490) and xanthophyll metabolic process (GO:0016122), which have been reported playing roles in disease resistance [24,25], were also significantly enriched during GO analysis. In summary, these results suggest that circRNAs enriched in these GO terms associated with metabolism and defense may be involved in the regulation of disease resistance.

CircRNAs Act as Putative miRNA Sponges
Besides regulating parental genes, circRNAs have been proposed to play important roles in miRNA function and transcription regulation by acting as competing endogenous RNAs [13]. miRNAs had been reported to participate in lots of plant stress resistance processes through regulating target genes' expression [26]. Obtaining insight into the miRNAs targeted by circRNAs will help us to further understand the functional importance of circRNAs. In this study, miRNAs Cell wall metabolism is an important aspect in the basal disease resistance process [23]. In our study, we found that lots of circRNAs were produced by protein-coding genes that are involved in cell wall metabolism. For example, PGSC0003DMG400021603, annotated as extensin, can generate 165 circRNAs. And PGSC0003DMG400009783, annotated as vegetative cell wall protein, can give rise to 46 circRNA (Table S3). The high abundance of circRNAs derived from the cell wall related genes implies its biological importance. Consequently, it is plausible to assume that circRNAs play important roles during the interaction between potato and pathogen Pcb. Moreover, pectin catabolic process (GO:0045490) and xanthophyll metabolic process (GO:0016122), which have been reported playing roles in disease resistance [24,25], were also significantly enriched during GO analysis. In summary, these results suggest that circRNAs enriched in these GO terms associated with metabolism and defense may be involved in the regulation of disease resistance.

CircRNAs Act as Putative miRNA Sponges
Besides regulating parental genes, circRNAs have been proposed to play important roles in miRNA function and transcription regulation by acting as competing endogenous RNAs [13]. miRNAs had been reported to participate in lots of plant stress resistance processes through regulating target genes' expression [26]. Obtaining insight into the miRNAs targeted by circRNAs will help us to further understand the functional importance of circRNAs. In this study, miRNAs sponged by DE circRNAs were used to predict the target mRNAs. As a result, 138 circRNAs (among the 470 DE circRNAs) were found to contain miRNA binding sites, suggesting that these circRNAs may function as miRNA sponges in potato. Target prediction showed that 120 miRNAs were the targets of these 138 circRNAs (Table S6). Furthermore, 608 mRNA targets of 96 miRNAs were predicted (Table S7). We then conducted GO analysis on these targeted mRNAs. Similar to GO analysis results in the parental genes of circRNAs, defense response (GO:0006952) was specifically enriched ( Figure 3B, Table S8). Besides, phosphorylation (GO:0016310), ADP binding(GO:0043531), and kinase activity (GO:0016301) were also specifically enriched. Transcription regulations play important roles in plant response to environmental stresses, including biotic and abiotic stresses [27]. And it is well known that a large amount of predicted miRNA targets in potato are transcription factors (TFs) [28]. Similarly, in this study, many of the predicted miRNA targets are various TFs, including bZIP, WRKY, and ethylene-responsive TF (Table S7). Gu et al. [29] reported that stu-miRNA5303, the specific miRNA family of Solanaceous plants, possesses target genes that are functionally related to responsive to abiotic stress, metabolic enzymes, and uncharacterized proteins. Interestingly, in this study, stu-miRNA5303, was sponged by 151 circRNAs ( Figure S1). Thus, we speculate that these 151 circRNAs may regulate the target genes of stu-miRNA5303 in potato immune response by sponging stu-miRNA5303 family members.

mRNA-CircRNA-LincRNA Co-Expression Networks
WGCNA has been demonstrated to identify important genes associated with complex phenotypes and biological processes in plants [30]. It can identify clusters of highly co-expressed genes based on gene expression similarity [30]. To reveal the potential functions of circRNAs, a co-expression network analysis to 8928 DE genes was performed (Table S9). Totally, 26 modules which captured 6858 DE genes (5466 coding genes, 299 lincRNAs, and 1087 circRNAs) and corresponded to the clusters of 38 to 1111 highly co-expressed individuals, were identified (Table S10). These modules can be divided into three types based on their expression pattern (Figure 4). sponged by DE circRNAs were used to predict the target mRNAs. As a result, 138 circRNAs (among the 470 DE circRNAs) were found to contain miRNA binding sites, suggesting that these circRNAs may function as miRNA sponges in potato. Target prediction showed that 120 miRNAs were the targets of these 138 circRNAs (Table S6). Furthermore, 608 mRNA targets of 96 miRNAs were predicted (Table S7). We then conducted GO analysis on these targeted mRNAs. Similar to GO analysis results in the parental genes of circRNAs, defense response (GO:0006952) was specifically enriched ( Figure 3B, Table S8). Besides, phosphorylation (GO:0016310), ADP binding(GO:0043531), and kinase activity (GO:0016301) were also specifically enriched. Transcription regulations play important roles in plant response to environmental stresses, including biotic and abiotic stresses [27].
And it is well known that a large amount of predicted miRNA targets in potato are transcription factors (TFs) [28]. Similarly, in this study, many of the predicted miRNA targets are various TFs, including bZIP, WRKY, and ethylene-responsive TF (Table S7). Gu et al. [29] reported that stu-miRNA5303, the specific miRNA family of Solanaceous plants, possesses target genes that are functionally related to responsive to abiotic stress, metabolic enzymes, and uncharacterized proteins. Interestingly, in this study, stu-miRNA5303, was sponged by 151 circRNAs ( Figure S1). Thus, we speculate that these 151 circRNAs may regulate the target genes of stu-miRNA5303 in potato immune response by sponging stu-miRNA5303 family members.

mRNA-CircRNA-LincRNA Co-Expression Networks
WGCNA has been demonstrated to identify important genes associated with complex phenotypes and biological processes in plants [30]. It can identify clusters of highly co-expressed genes based on gene expression similarity [30]. To reveal the potential functions of circRNAs, a coexpression network analysis to 8928 DE genes was performed (Table S9). Totally, 26 modules which captured 6858 DE genes (5466 coding genes, 299 lincRNAs, and 1087 circRNAs) and corresponded to the clusters of 38 to 1111 highly co-expressed individuals, were identified (Table S10). These modules can be divided into three types based on their expression pattern (Figure 4).  Genes clustered in type I modules were characterized as suddenly down-regulated in either cultivar. In modules 1, 2, and 3, genes were down-regulated in Valor, but not in BP1 (Figure 4). Three modules enriched genes for photosynthesis and response to biotic stimulus (Table S11). For example, lots of receptor-like protein kinases (RLKs), which play important roles in activating plant innate immunity [31], were enriched in module 2. In module 2, RLKs were suddenly down-regulated in susceptible cultivar Valor, but slightly up-regulated in disease tolerant cultivar BP1 after 6 hpi, suggesting that rapid immune response in Valor were suppressed, whereas BP1 has a more rapidly response to the infection. In modules 4-6, genes that are functionally related to negative regulation of metabolism, were down-regulated in BP1, but not in Valor, which is consistent with disease tolerant phenotype of cultivar BP1 (Figure 4).
Genes clustered in 14 type II modules were characterized as only up-regulated in one cultivar, BP1 or Valor. Rapidly transcriptome reprograming is crucial for the immune response of host to pathogen infection [32]. Obviously, in BP1, genes clustered in modules (e.g., [14][15][16][17] were up-regulated immediately at 6 hpi, whereas, in Valor, genes clustered in modules (e.g., 8, 9, 10, and 11) were up-regulated after 24 hpi. The distinct regulating patterns indicated the different defense response of two cultivars. Specifically, in BP1, module 17 formed a large cluster containing 1111 individuals, and among them, 1045 genes are functionally related to defense response, such as numerous RLKs (receptor like kinases) (e.g., wall-associated kinases, WAK1; brassinosteroid insensitive 1-associated receptor kinase, BAK1. Table S12), L-phenylalanine metabolic process (GO:0006558), and cell wall thickening (GO:0052482), which were reported to have vital roles in plant growth and disease resistance ( Figure S2, Table S9) [31,33,34]. Module 16 showed a similar expression pattern as 17, and clustered genes were functionally related to response to ethylene (GO:0009723). The production of ethylene (ET) was the resulting hallmarks of pattern-triggered immunity (PTI) [35], and ET signaling is required for oxidative burst contributing to plant immunity in PTI [36]. Cuticle provides a physical barrier against pathogens' infection and therefore plays a basal role in disease resistance [37]. Module 15 contains genes functionally associated with cuticle development (GO:0042335) and cell wall organization (GO:0071555). In summary, these results indicated that the disease tolerance cultivar BP1 has a more quick and effective response to the pathogen infection, as these genes functionally related to cell wall, lignin synthesis, and ET were rapidly, robustly, and selectively transcriptionally reprogramed, which finally strengthen cell wall and launch PTI (pathogen-associated molecule patterns, PAMP triggered immunity) in cultivar BP1.
Genes clustered in type III modules showed disordered expression pattern. Nevertheless, interestingly, genes clustered in module 22 show higher expression level in BP1 than in Valor at each infection stage. These genes were functionally related to lignin metabolic process (GO:0009808), such as cellulose synthases and NAC (standing for no apical meristem (NAM), arabidopsis transcription activation factor (ATAF) and cup-shaped cotyledon (CUC)) domain-containing protein that are involved in regulation of secondary wall biosynthesis and cell death (Table S13) [38]. The cultivar-specific expression patterns suggested their possibly defense-related roles in enhancing resistance to Pcb.
Based on the co-expression relationship of circRNAs, lincRNAs, and coding genes, we inferred the potential functions of circRNAs in potato responsive to Pcb infections. Based on the co-expression analysis, we constructed a circRNA-centered (circRNA chr02:43604462|43622094) subnetwork, in which circRNA directly connected with acquired gene sets. In the subnetwork, circRNA chr02:43604462|43622094 was connected to many RLKs, TFs, and PRR (Pattern recognition receptor) coding-genes ( Figure S3). These genes were GO enriched in L-phenylalanine metabolic process and cinnamic acid metabolic process that are closely related to cell wall composition. Besides, numerous lincRNAs were also connected with chr02:43604462|43622094, suggesting that lincRNAs may also participate in regulating the rapid transcriptional reprogramming involved in plant defense response. The networks analysis revealed the relationships among circRNAs, lincRNAs, and coding genes, and suggested the possible roles of circRNAs in transcriptome reprogramming during disease response.

Validation of the Identified circRNAs
To validate the circRNAs identified from the circRNA-Seq analysis, reverse transcription RT-PCR and Sanger sequencing experiments were performed to 14 randomly selected circRNAs. A pair of divergent primers (Table S14) were designed for each circRNA, and both cDNA and gDNA (negative control) were used as template for PCR amplification ( Figure 5). Finally, among the 14 circRNAs, 13 were successfully confirmed (92.86%), suggesting the reliability of our circRNAs identification result (Table S14), as comparing to kiwifruit (68/80, 85.00%), tomato (8/11, 72.73%), and rice (10/18, 55.56%), much higher rate of circRNAs were successfully validated in our study [10,12,13]. genes, and suggested the possible roles of circRNAs in transcriptome reprogramming during disease response.

Identification and Expression Analysis of Circular RNAs
Trimmed reads were mapped to potato genome by BWA (v0.7.15, mem -T 19) [45]. The SAM file of alignment was then inspected using CIRI2 (v2.0.5) to identity circRNAs with the default parameters [46]. To compare the expression of circRNAs across time and taxa, we calculated the counts of backspliced reads for each circRNA and normalized the expression levels by comparing counts with mapped reads in corresponding sample (defined as reads per million mapped reads, RPM). DE circRNAs were defined as |log 2 (foldchanges)| ≥ 1.

Prediction of miRNA Target
The FASTA formatted circRNA sequences were used as queries and input into psRNAtarget (http://plantgrn.noble.org/psRNATarget/) to predict their sponging miRNAs [47]. Then miRNAs sponged by circRNAs were used as queries to predict their target mRNAs using TargetFinder through searching potato coding genes assembled by cufflinks [48].

Gene Co-Expression Network Analysis
The co-expression network analysis was performed to group DE genes (including coding genes, lincRNAs, and circRNAs) into modules using R package WGCNA (v1.51) [49]. In brief, we applied an optimal power function (β) as 20 to balancing the scale-free property of the co-expression network and sparsity of connection between genes. And the cutreeDynamic function was performed with parameters as deepSplit = 2, pamRespectDendro = F, minClusterSize = 30, and finally using the mergeCloseModules with cut Height set to 0.25 to merge these highly correlated clusters. Modules were visualized using Cytoscape (v3.4.0), setting the layout with edge-weighted spring embedded layout [50].

Reverse Transcription (RT)-PCR
To validate the identified circRNAs in potato, total RNA was extracted using Plant RNA Kit (Omega, London, UK) and cDNA was synthesized by reverse transcription using RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, Waltham, MA, USA). The divergent PCR primers were designed using the "out-facing" strategy to exclude linear mRNA from amplification (synthesized by Sangon Biotech Co., Ltd., Shanghai, China; Table S14). PCR products were validated by Sanger sequencing.

Conclusions and Perspectives
CircRNAs constitute a family of transcripts with unique structures and functions that are still largely unknown, especially in plants. In this study, through the expression patterns analysis and gene function annotation, we speculated that circRNA may participate in the regulation of potato immune response by regulating parental genes as well as acting as putative miRNA sponges. As GO enrichment showed that circRNAs were predicted to be involved in defense response (GO:0006952), cell wall (GO:0005199), ADP binding (GO:0043531), phosphorylation (GO:0016310), and kinase activity (GO:0016301). Furthermore, WGCNA found circRNAs being closely related with mRNAs and lincRNAs. And together they were cultivar-specifically regulated to strengthen immune response of potato to Pcb infection, implying the roles of circRNAs in reprogramming disease responsive transcriptome. The results reported in the present study lay a foundation for further studies of the roles that circRNAs play in potato immune response. The precise functions of circRNA in potato immune response need to be experimentally tested in the following researches.