Transcriptomic de novo analysis of pitaya (Hylocereus polyrhizus) canker disease caused by Neoscytalidium dimidiatum
BMC Genomics volume 20, Article number: 10 (2019)
Canker disease caused by Neoscytalidium dimidiatum is the most serious disease that attacks the pitaya industry. One pathogenic fungus, referred to as ND8, was isolated from the wild-type red-fleshed pitaya (Hylocereus polyrhizus) of Hainan Province. In the early stages of this disease, stems show little spots and a loss of green color. These spots then gradually spread until the stems became rotten due to infection by various strains. Canker disease caused by Neoscytalidium dimidiatum poses a significant threat to pitaya commercial plantations with the growth of stems and the yields, quality of pitaya fruits. However, a lack of transcriptomic and genomic information hinders our understanding of the molecular mechanisms underlying the pitaya defense response.
We investigated the host responses of red-fleshed pitaya (H. polyrhizus) cultivars against N. dimidiatum using Illumina RNA-Seq technology. Significant expression profiles of 23 defense-related genes were further analyzed by qRT-PCR. The total read length based on RNA-Seq was 25,010,007; mean length was 744, the N50 was 1206, and the guanine-cytosine content was 44.48%. Our investigation evaluated 33,584 unigenes, of which 6209 (18.49%) and 27,375 (81.51%) were contigs and singlets, respectively. These unigenes shared a similarity of 16.62% with Vitis vinifera, 7.48% with Theobroma cacao, 6.6% with Nelumbo nucifera and 5.35% with Jatropha curcas. The assembled unigenes were annotated into non-redundant (NR, 25161 unigenes), Kyoto Encyclopedia of Genes and Genomes (KEGG, 17895 unigenes), Clusters of Orthologous Groups (COG, 10475 unigenes), InterPro (19,045 unigenes), and Swiss-Prot public protein databases (16,458 unigenes). In addition, 24 differentially expressed genes, which were mainly associated with plant pathology pathways, were analyzed in-depth.
This study provides a basis for further in-depth research on the protein function of the annotated unigene assembly with cDNA sequences.
Pitaya, also known as dragon fruit, is an important tropical and subtropical fruit tree. At present, there are two known species of pitaya: the red-fleshed pitaya (Hylocereus polyrhizus) and the white-fleshed pitaya (Hylocereus undatus). However, the most cultivated and commercialized species is the red-fleshed pitaya with its red pulp and pericarp (H. polyrhizus) that is enriched with anthocyanins [1, 2]. Pitaya belongs to the Cactaceae and its native tropical zone is Latin America, which includes Mexico, and Central and South America . Currently, it is also extensively grown throughout tropical and subtropical regions worldwide , especially in Asian countries, such as Malaysia, Vietnam, Thailand, Philippines and China. Presently, the main areas of pitaya cultivation (especially the red dragon fruits) include the Hainan, Guangxi, Guangdong, Yunnan, Fujian and Taiwan provinces of China.
Pitaya has been proven to be very adaptable to environmental conditions and disease pressures. However, the rapid spread of stem canker, a disease caused by Neoscytalidium dimidiatum among red-fleshed pitaya (H. polyrhizus), has become a serious threat to some commercial plantations (Fig. 1). The disease was first reported in Taiwan province of China, and features a distinctive canker on stems and fruits, and quickly spread to most commercial cultivation areas . Shortly thereafter, it was discovered and reported in Guangdong Province of China, Malaysia, Israel, and Florida with similar characteristics [6,7,8,9,10]. Various studies have reported the incidence of canker disease caused by Neoscytalidium dimidiatum and have provided basic information relating to the identity and isolation of pathogenic strains.
Generally, two main recognition systems are associated with the response to an attack by the microorganism. One system recognizes the pathogen-associated molecular patterns (PAMPs) via pattern-recognition receptors (PRRs), which are localized on the membrane and activate an innate immune response, or more specifically PAMP-triggered immunity (PTI) . The PAMPs are structural molecules, such as bacterial peptidoglycan and flagellin, oomycete glucans, lipopolysaccharide harpin, xylanase and chitin [12, 13]. In contrast, PRRs are cell surface-localized receptors, most of which are receptor-like kinases (RLKs) or receptor-like proteins (RLPs). The RLKs contain extracellular, transmembrane, and cytoplasmic kinase domains, whereas the RLPs lack the kinase domain . A burst of reactive oxygen species (ROS) is a typical reaction in early cellular PTI events; activation of the mitogen-activated protein kinase (MAPK) cascade and expression of immune-related genes are also associated with such events . The second system is mediated by intracellular immune receptors, which directly or indirectly recognize virulence effectors of pathogens that are secreted within host cells, thereby inducing effector-triggered immunity (ETI). In the second system, pathogen effectors are recognized by intracellular receptors, which are referred to as resistance (R) proteins, and mainly include leucine-rich repeats (LRR) or the LRR-like family, toll/interleukin-1-receptor (TIR), and serine/threonine kinases (S/TK) to activate the ETI response [14, 16]. Compared with PTI, the activation of ETI by receptors localized to the nucleus appears to be more directly associated with transcriptional regulation of the expression of defense genes . However, some ambiguous theories exist regarding PTI and ETI, as it also appears that both can be robust or weak, depending on the interaction, defense signaling pathways, receptors and environmental conditions .
The RNA-Sequencing (RNA-Seq) technique is an excellent method that permits quick and extensive comprehension of the transcriptional level of genes in a variety of plant species using next-generation sequencing technologies . Here, to investigate the genes associated with the defense response of pitaya to canker disease caused by Neoscytalidium dimidiatum in PTI and ETI systems, we used RNA-Seq technology to monitor the transcriptomic profiles of red-fleshed pitaya (H. polyrhizus) in response to fungal invasion. A comparison was made between normal and diseased stem tissuses during transcriptomic analysis, and the significant differentially expressed genes (DEGs) were analyzed in detail. These DEGs were annotated in public protein databases, suggesting that they may play important roles in PTI and ETI systems in response to the invasion of Neoscytalidium dimidiatum. The aim of the present study was to identify the key genes that encode steps in the plant pathological pathways of pitaya. This article provides a solid basis for a better understanding of the transcriptomic response of the red-fleshed pitaya (H. polyrhizus) to canker disease.
Sequencing, de novo assembly and annotation
Healthy and diseased stem tissues of red-fleshed pitaya (H. polyrhizus) (Fig. 2) were collected from the orchard of Ledong County, Hainan Province. Three normal samples (N1, N2, and N3), and three diseased samples (D1, D2, and D3), were used for transcriptional sequencing. However, N1 and D2 examples were abandoned because the results were not consistent. The unigenes assembled from the four remaining samples were then merged. In total, we evaluated 33,584 unigenes, of which 6209 (18.49%) and 27,375 (81.51%) were contigs and singlets, respectively. The total read length was 25,010,007, mean length was 744 bp, the N50 was 1206, and the guanine-cytosine (GC) content was 44.48%, based on RNA-Seq. The range of the unigene lengths in the four samples was 497–702 bp. Details relating to the assembled information from each sample and the merged results are shown in Table 1. In total, 1608 DEGs, with a fold change ≥2.00 or a fold change ≤0.5, and a false discovery rate (FDR) ≤ 0.001, were monitored. On average, the number of unigenes from D samples was 45.7% higher than other samples, reflecting a significant increase in the number of sequences resulting from fungal infection.
Owing to the lack of a reference genome for H. polyrhizus, the assembled unigenes were blast searched into the Non-redundant (NR), Kyoto Encyclopedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups (COG), InterPro, and SwissProt public protein databases using search tools. The species with the greatest number of H. polyrhizus unigenes were Vitis vinifera (16.62%), Theobroma cacao (7.48%), Nelumbo nucifera (6.6%), and Jatropha curcas (5.35%), as shown in Fig. 3. Creation of a Venn diagram (Fig. 4), showed that there were 7305 unigenes annotated in all databases, 25,161 unigenes annotated in the NR database, 17,895 unigenes annotated in the KEGG database, 10,475 unigenes annotated in the COG database, 19,045 unigenes annotated in the InterPro database, and 16,458 Unigenes annotated in the SwissProt database. These findings provide a basis for further in-depth research on protein function of the annotated unigenes, which were assembled with cDNA sequences.
Gene ontology (GO) functional annotation of DEGs
Based on GO functional annotation, approximately 1608 DEGs were classified into three essential categories, namely biological process, cellular component, and molecular function. The unigenes we identified covered a wide range of functional GO categories. Of the 1608 DEGs classified, most were associated with catalytic activity (173 DEGs, 10.75%), binding (134 DEGs, 8.33%), and transporter (30 DEGs, 1.87%) aspects of molecular function. The biological processes were associated with cellular processes (160 DEGs, 9.95%), metabolic processes (157 DEGs, 9.76%), and single-organism processes (124 DEGs, 7.71%). The cellular component processes were mainly associated with the membrane (134 DEGs, 8.33%), cell (107 DEGs, 6.65%), membrane part (106 DEGs, 6.59%), and cell part (104 DEGs, 6.47%) aspects (Fig. 5). Analysis of the GO terms of pitaya transcripts showed that these genes might have been associated with kinases (catalytic activity) and receptors (cellular and metabolic processes, and the membrane).
DEGs associated with plant-pathogen interactions
Twenty-four DEGs were associated with plant-pathogen interaction pathways; of these, 23 were up-regulated and one was down-regulated. These unigenes play different roles in the plant-pathogen interaction pathways. Unigene13603_All and Unigene11517_All were annotated in the cyclic nucleotide-gated channel (CNGC) family, which mediates Ca2+ influx from cellular stores in plants . The CNGCs are reportedly Ca2+-permeable channels that interact with the ubiquitous Ca2+ sensor calmodulin, and play important roles in the response to plant hormones, biotic and abiotic stresses, and plant immunity [20,21,22]. Unigene12122_All, Unigene16412_All, and CL2472.Contig2_All were annotated in the calcium-dependent protein kinase family (CDPK). The CDPKs are crucial sensors of changes in calcium concentration and play multiple roles in plant development and growth, abscisic acid (ABA)-mediated processes, jasmonic acid (JA) biosynthesis, plant tolerance to stress, and plant fungal stimuli interaction [23,24,25,26,27,28]. The CDPKs are phosphorylated by phosphatase and participate in the hypersensitive response (HR). Furthermore, another calcium-related gene, Unigene7838_All, was annotated in two types of Ca2+-sensing proteins, calmodulin and calmodulin-like (CaM/CML) (Fig. 6). The CaM/CML protein family plays complex and versatile roles in modulating cellular responses to various stimuli, particularly biotic stresses in grapevine . The ROS produced by respiratory burst oxidase homologs are associated with multiple signal transduction pathways in diverse biological processes in plants . In our research, an up-regulated Unigene (13565_All) was annotated in the respiratory burst oxidase homologs (Rboh) membrane-bound enzyme family, which is also known as the NADPH oxidases . In tomato plants, Rboh activity and H2O2 signaling are important components of melatonin-induced stress tolerance . This unigene may therefore encode the most critical enzyme of pitaya associated with ROS production under conditions of Neoscytalidium dimidiatum infection.
A plant-pathogen interaction pathway (Fig. 6) showed that the flagellin-sensing 2 (FLS2) family was significantly up-regulated, and included Unigene13867_All, Unigene9087_All, CL1260.Contig2_All, Unigene13234_All, Unigene24332_All, and Unigene15298_All. Furthermore, FLS2 is a well-known PRR kinase that recognizes a conserved 22 amino acid N-terminal sequence of the bacterial flagellin protein (flg22) . Previous research has shown that FLS2 contains three domains, which include an extracellular LRR domain, a transmembrane domain and a cytoplasmic kinase domain . The FLS2 receptor complex is directly coupled with Arabidopsis heterotrimeric G proteins to regulate immune signaling via both pre-activation and post-activation mechanisms . Our transcriptome data showed that LRR receptor-like serine/threonine-protein kinases (LRR-RLKs) were differentially expressed and included Unigene13867_All, CL1043.Contig1_All, and Unigene9087_All. The LRR-RLKs play pivotal roles in plant growth, development, cellular signal transduction, brassinosteroid (BR) and ABA signal pathways, and responses to biotic and abiotic stresses . The unigenes we detected that were associated with FLS2 may play important roles in the HR and defense-related gene induction in response to fungal infection.
Moreover, three LRR classes of plant disease resistance genes PRS  (Unigene8050_All, Unigene15133_All, and Unigene14674_All); a resistance PBS gene (Unigene24332_All) that encodes a serine/threonine-protein kinase ; a pto-interacting protein gene PTI (Unigene24688_All) that encodes a serine/threonine kinase ; and a receptor for the fungal elicitor ethylene-inducing xylanase (EIX) that binds with EIX to induce the HR in tobacco and tomato  (Unigene24332_All), were significantly up-regulated in pathogen interaction pathways. These DEGs may participate in the HR, defense-related gene induction or programmed cell death (Fig. 6). The expression profiles of DEGs in four samples asociated with plant-pathogen interactions are presented in the pheatmap show in Fig. 7, and in Table 2.
Predicted expression of transcription factors (TFs)
Transcription factors (TFs) are proteins that, when in contact with RNA polymerase, can confirm a transcription initiation complex and participate in the transcription initiation process. Generally, a TF contains a transcriptional activation/repression domain and a transcriptional binding domain . As a regulator, a TF regulates the expression of genes downstream of a pathway by binding with DNA sequences either by itself or by interacting with other TFs or proteins. TFs are associated with the significant mechanisms of the regulation of disease resistance genes and receptors for cellular localization and activation, as well as DNA binding and transcription activities . In the merged data of all samples in the present study, a total of 832 transcription factors were predicted, and the majority of the TF families were categorized as MYB (95, 11.42%); MYB-related (77, 9.25%); APETALA2/ethylene-responsive element binding protein (AP2-EREBP, 64, 7.69%); C3H (56, 6.73%); WRKY (52, 6.25%); bHLH (42, 5.05%); NAC (37, 4.45%); and GRAS (34, 4.09%), among others (Fig. 8).
The TFs make essential contributions to the processes of growth, development, flowering, and biotic and abiotic stresses in the plant. For instance, the family of MYB factors include the conserved MYB DNA-binding domain and appears to be involved in the regulation of the cell cycle in animals, plants and other higher eukaryotes . The AP2-EREBP genes contained the AP2 family, RAV (related to ABI3/MP) family, and ERF (ethylene response factor) family, based on the number of AP2 domains and sequence similarity . The AP2-EREBP gene reportedly plays roles in controlling flower and seed developmental processes, regulatory networks in response to hormones, pathogen attacks and environmental signals [45, 46]. Moreover, NAC TFs reportedly also take part in the plant-pathogen interactions of Arabidopsis, potato, maize, and wheat in response to wounding or pathogen attack as a transcriptional regulator . In our plant pathology interaction pathways, WRKYs (Unigene9452_All, Unigene2169_All, Unigene14674_All, Unigene2538_All, and Unigene16508_All) were significantly and differentially expressed. The WRKY proteins can interact with themselves, VQ-domain proteins, chromatin remodeling proteins, calmodulin-like proteins or some essential MAP kinases . The unigenes identified in the present study may act as downstream products of the MAPK pathway to activate the expression of defense-related genes, or recognize the fungal effector, which may result in the induction of defense-related genes or programmed cell death (Fig. 6). The differentially expressed TFs under investigation may play essential roles in growth, development, leaf senescence and responses to biotic and abiotic stimuli.
To validate the RNA-Seq data, 23 unigenes (Table 3) that were significantly expressed, most of which were associated with phytopathology, were selected for additional qRT-PCR analysis. These unigenes included one ethylene-responsive TF (Unigene12062_All); two LRR receptor-like serine/threonine-protein kinases (Unigene13867_All and Unigene9087_All); one salicylic acid-binding protein (CL2401.Contig2_All); two CDPK-related kinases (Unigene12122_All and Unigene16412_All); three disease resistance proteins (Unigene2329_All, Unigene8050_All, and Unigene15133_All); two RLPs (Unigene24332_All and Unigene15298_All); four WRKY TFs (Unigene9452_All, Unigene2169_All, Unigene2538_All, and Unigene16508_All); one auxin response factor (Unigene27293_All); one NADPH-related protein (Unigene24688_All); one lipoxygenase protein (CL1613.Contig2_All); and five differentially expressed proteins (CL2068.Contig2_All, CL1097.Contig1_All, CL1813.Contig1_All, Unigene2331_All, and Unigene9233_All). These genes were either significantly up- or down-regulated in diseased tissues. The results of the qRT-PCR assay showed that 20 of the 23 genes selected for qRT-PCR (three biological replicates) conformed to the RNA-Seq data, with an accuracy of 87%. These results (Fig. 9) indicated that the RNA-Seq findings were both stable and reliable.
The cultivation of pitaya is an emerging industry with a bright future and high commercial value, especially for red-fleshed pitaya (H. polyrhizus). However, canker disease caused by Neoscytalidium dimidiatum presents the most devastating attacks on pitaya plantations at present . The most efficient solution to prevent this disease is to cut off the infected tissues and bury them in soil away from the plantation. In addition, molecular breeding is the best way to acquire disease-resistant cultivars; however, this is not easily done with pitaya. Moreover, additional tasks need to be completed to verify which cultivar is the most resistant among a variety of pitaya species.
Plant immunity is a complicated process that is effected mainly in response to extracellular or intracellular receptors that recognize PAMPs or effectors in PTI and ETI systems [13, 14]. The defense response mechanisms and roles of the receptors have been previously reported in Arabidopsis, tomato, and rice, among others [34, 35, 37, 39, 40, 49, 50], but has been scarcely described in pitaya. In our transcriptome results, higher expression levels of unigenes were detected in H. polyrhizus and 24 DEGs associated with the plant-pathogen interaction pathway. Most of these genes are signal receptors, protein kinases, and TFs, which may play critical roles in the plant-pathogen interaction networks.
Calcium acts as an important conserved second messenger in plant immune and stress responses . Calcium is sensed by Ca2+-binding proteins in signaling networks, including calmodulin (CaM), calmodulin-like (CML), CDPKs, and calcineurin B-like proteins [19, 52]. The CNGCs are calmodulin-permeable cation transport channels that have the potential to integrate signals from cyclic nucleotide and Ca2+ signaling pathways . In our transcriptome analysis, seven unigenes (Table 2) related to Ca2+ levels and signaling pathways were annotated and showed significant changes in expression levels in diseased red-fleshed pitaya. These genes included two CNGCs related genes (Unigene13603_All and Unigene11517_All); one CaM/CML (Unigene7838_All); three CDPK genes (Unigene12122_All, Unigene16412_All, and CL2472.Contig2); and one Rboh family gene (Unigene13565_All). These genes act as receptors or kinases in the HR, cell wall reinforcement, and stomatal closure (Fig. 6).
Pathogenesis-related (PR) proteins comprise one of the major sources of plant-derived allergens, and are generally induced by various types of pathogens, such as viruses, bacteria and fungi . In pepper, PR4b reportedly interacts with LRR proteins and PR4c, as a plasma membrane-localized cysteine protease inhibitor, to trigger cell death and the defense signaling response [55, 56]. In the present study, we also monitored four differentially expressed PR homologous proteins, which were annotated in the SwissProt and InterPro public protein databases. Among these four proteins, two were up-regulated (Unigene9313_All and Unigene11745_All), and the other two were down-regulated (CL1813.Contig1_All and CL1271.Contig2_All). Further work, involving protein function verification and genetic transformation, is now being carried out to investigate their roles in the response to the fungal infection of H. polyrhizus.
A total of 832 TFs were monitored in the present study. The TF proteins are key factors in the regulatory networks that control development, metabolism, and the responses to biotic and abiotic stresses in numerous plant species. In our RNA-Seq data, MYB and MYB-related family genes comprised the largest group, which accounted for 20.67%. A MYB domain protein, COLORED1 (C1), which is associated with anthocyanin synthesis, was first identified in maize . The MYB TF family is a major player in drought and cold responses [58, 59], the phenylpropanoid metabolic pathway , and JA and ABA signal transduction pathways [61, 62]. The NAC TFs comprise a large protein family, with 138 putative NAC genes in Arabidopsis, as listed in the Plant Transcription Factor Database (PlantTFDB, http://planttfdb.cbi.pku.edu.cn/) . However, 37 (4.45%) putative NAC genes in H. polyrhizus were identified based on our RNA-Seq data. The NAC family of genes have a highly conserved N-terminal DNA-binding domain and a variable C-terminal domain. The NAC TFs participate in the response to drought, salinity, and/or low temperatures and regulate JA-signaled defense responses in various plants [64, 65]. In addition, 52 WRKY superfamily TF genes were monitored in the present study and five unigenes (Unigene 9452_All, Unigene 14674_All, Unigene 2538_All, Unigene 2169_All, and Unigene 16508_All) were differentially expressed upstream in red pitaya. The WRKY genes are quite common in plants, which contain a conserved N-terminal sequence of WRKYGQK, along with a Zn finger-like motif (C2H2 type or C2HC type) . The trend in qRT-PCR expression levels of four WRKY unigenes (Unigene 9452_All, Unigene 2169_All, Unigene 16508_All, and Unigene 2538_All) were consistent with the RNA-Seq data (Fig. 10 and Table 3). The Unigene 14674_All was an unnamed protein product, which was annotated to both WRKY and RRS1-R gene families. Pathway results suggested that the five differentially expressed WRKY unigenes were important candidate genes of H. polyrhizus involved in plant pathology interactions (Fig. 6).
Our findings may pave a fundamental step to a better understanding of the defense responses of pitaya. We have provided hypotheses for the functions of the genes identified, according to public protein annotation databases. However, specific details of the roles of these genes in pitaya defense responses pose a relevant and interesting challenge. Thus, further research is now required to clarify how significant DEGs regulate pitaya defense processes via genetic transformation, when infected by Neoscytalidium dimidiatum.
At present, transgenic research is the best method with which to clarify the functions of a gene. However, this is difficult to achieve in tropical and subtropical non-model plants. Three major difficulties are associated with transgenic research in pitaya. Firstly, there have been problems associated with the establishment of a tissue culture system. Secondly, we need to consider the stability of the tissue culture system and the resolution of issues if stability is unfavorable. Finally, the period required to screen a disease-resistant cultivar is relatively long. Nevertheless, in spite of the difficulties faced in the study of non-model plants, transgenic Cavendish bananas with resistance to Fusarium wilt tropical race 4 were reported recently following a 3-year field trial . That was a major breakthrough for the banana industry and could also imply a brighter future for the pitaya industry. Thus, the problems associated with transgenic pitaya should be solved eventually with further broad-scale research.
Plant materials and culture conditions
Healthy and diseased stem tissues of red-fleshed pitaya (H. polyrhizus) were collected from the orchard of Ledong County, Hainan Province. These plants had been growing for approximately five years in a plantation under natural conditions. Fungal spores were isolated from diseased tissues, and morphological and molecular analyses were conducted to determine the species of Neoscytalidium dimidiatum 8 (ND8). The article about isolation of ND8 has been published in Australasian Plant Pathology . The collected stem tissues were divided into two groups with similar characteristics in triplicate within each group. The stems were cut into pieces and immediately frozen in liquid nitrogen. The three replicates of the normal group were named N1, N2, and N3, whereas the three replicates of the diseased group were named D1, D2, and D3. Three biological replicates, each consisting of four stem cutting bases, were frozen in liquid nitrogen and stored at − 80 °C until further analysis. Sample tissues were transferred to Shenzhen BGI Tech Company (Shenzhen, China) for transcriptome de novo assembly and sequencing using an Illumina HiSeq system.
Total RNA extraction
Tissues (100 mg) were placed into a precooled mortar to be ground into a powder with liquid nitrogen. The ground stem powder was collected into a precooled 2 mL microcentrifuge tube containing 2 mL RNase-free cetyltrimethylammonium bromide (CTAB) extracting solution, which contained 0.1% diethyl pyrocarbonate (DEPC), 2.5% CTAB, 1.4 M NaCl, 20 mM EDTA, 100 mM Tris-HCl (pH 8.0), and 4% polyvinylpyrrolidone (PVP). The tubes were then placed in a water bath at 65 °C for 30 min to induce cell lysis, and centrifuged at 13,000 rpm for 5 min at 25–28 °C. The residue was discarded, and the supernatant was transferred into a clean 2 mL microcentrifuge tube. Trichloromethane (1 mL) was added to the supernatant solution, mixed lightly, and centrifuged at 13000 rpm for 15 min at 4 °C. The supernatant was then transferred into a clean 2 mL microcentrifuge tube, and this procedure was repeated once more. Subsequently, 8 M lithium chloride was added with RNase-free water to 0.6 volume of the supernatant and left overnight. The mixture was centrifuged for 30 min at 13000 rpm and 4 °C to obtain an RNA precipitate. RNase-free ethyl alcohol (75%) was used to wash the RNA precipitate twice. Finally, RNA was dissolved in RNase-free ddH2O and stored at − 80 °C until further use.
Construction of the cDNA library and de novo RNA-Seq
Poly (A) mRNA from the total RNA was enriched by poly-T oligo-attached magnetic beads (Invitrogen, Thermo Fisher Scientific, CN), according to the manufacturer’s protocol. Purified mRNA was then divided into short fragments by mixing with fragmentation buffer. First strand cDNA was synthesized using a random hexamer primer. Buffer, dNTPs, RNase H, and DNA polymerase I were then added to the fragmented mRNA templates to synthesize second-strand cDNA. The double-stranded cDNA was purified using magnetic beads. End repair, and the addition of 3′-end single nucleotide A (adenine), were then performed. Finally, sequencing adaptors were ligated to the fragments, which were enriched by PCR amplification. During the quality control step, a 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) and a StepOnePlus Real-Time PCR System (Applied Biosystems Inc., Foster City, CA, USA) were used to qualify and quantify sequences for the sample libraries. The constructed library products were finally prepared for sequencing using a HiSeq 2000 system (Illumina, San Diego, CA, USA) at the Beijing Genomics Institute (BGI).
Filtering of raw reads and the assembly of clean reads
Before any further analysis could be conducted, we needed to consider quality control of the data. In addition, low-quality sequences were filtered out of the raw data to reduce data noise. We defined “dirty” raw reads as those that contained adapter sequence, high levels of unknown bases, and low-quality reads. These needed to be removed before data analysis. Filtering steps were as follows: 1) Reads with adapters were removed; 2) reads in which unknown bases were more than 10% were removed; and 3) low-quality reads were removed (the percentage of low-quality bases is usually over 50% in a read, and we defined a low-quality base to be one with a sequencing quality of no more than 5). After filtering, the remaining reads were referred to as “clean reads” and used for downstream bioinformatics analysis. Clean data were obtained and stored in the FASTQ format after filtering.
Clean reads were assembled using the Trinity paired-end assembly method . The assembled clean reads, also referred to as transcripts, were clustered and this eliminated any redundancy in obtaining unigenes using the TIGR Gene Indices clustering tools (TGICL) . For the study of multi-samples, the previous step using the TGICL for downstream analysis of the unigenes was repeated. The processed unigenes were referred to as “All-Unigenes.” All-Unigenes were divided into two parts; one part included clusters of a few unigenes with more than 70 similarities among them (starting with “CL,” followed by the gene family ID); and the other included singletons (starting with “Unigene”). The clean data arising from high-throughput sequencing (RNA-Seq) results, and the processed files, have been submitted to the Gene Expression Omnibus (GEO) database at the GenBank website with accession number GSE119976. The process of preliminary bioinformatics analysis was performed according to the flow diagram in Fig. 10.
Functional annotation and the detection of differentially expressed unigenes
After assembly, H. polyrhizus pitaya All unigenes (contigs and singletons) were functionally aligned to the public protein databases, NR, GO, KEGG, COG, InterPro, and SwissProt. All unigenes were annotated according to their sequence similarity to previously annotated genes in these public databases. The phyper function of R software (https://www.r-project.org) was used for enrichment analysis, according to the results of GO and KEGG pathway gene annotations in functional analyses. The p-value was corrected by the FDR. We considered that unigenes with a FDR ≤ 0.01 were significantly expressed in the GO and KEGG pathways. Moreover, NOIseq  and PassionDis  methods were used to test DEGs. A fold change ≥2.00 or fold change ≤0.5, with a Probability ≥0.8 in NOIseq, and fold change ≥2.00 or fold change ≤0.5, with a FDR ≤ 0.001 in PassionDis of unigenes were used to define DEGs in the normal and diseased libraries. According to the results of the DEGs test, the pheatmap function of R software was used for hierarchical clustering analysis.
Verification of DEGs by quantitative reverse transcription PCR (qRT-PCR)
To investigate the stability and accuracy of RNA-Seq, 23 DEGs that were possibly associated with plant pathology, were selected for real-time quantitative RT-PCR (qRT-PCR). Total RNA extraction was performed as described in previous steps. Following treatment with DNase I (Thermo Fisher Scientific, USA), first strand cDNA synthesis was conducted according to the manufacturer’s protocols using the PrimeScript™ 1st strand cDNA synthesis kit (code No.6210A, Takara, Japan). qRT-PCR was performed using a 20-μL reaction system, containing 10 μL of SYBR® Premix Ex TaqII (Tli RNase Plus) (2×), 2 μL cDNA templates, 0.8 μL PCR forward primer (10 μM), 0.8 μL PCR reverse primer (10 μM), 0.4 μL ROX reference Dye II (50×) and 6 μL of nuclease-free water. The qRT-PCR assays (three technical replicates) were carried out using a 7500 Applied Biosystems qRT-PCR System (Life Tech, USA). The primers of all genes selected for qRT-PCR were designed using Primer Premier 6 software (www.premierbiosoft.com). The pitaya ubiquitin gene (UBQ) was used as an internal reference control for data normalization. The primers for UBQ were as follows: F: 5′-TGAATCATCCGACACCA-3′; and R: 5′-TCCTCTTCTTAGCACCACC-3′ . Reactions were carried out under the following conditions: 50 °C for 2 min, 95 °C for 10 min, followed by 40 cycles at 95 °C for 15 s, and 60 °C for 1 min. Output data were generated using the Applied Biosystems 7500 software version 2.0.6 (ABI, America). All reactions were performed with two independent biological replicates, and the expression levels calculated for each sample were based on three technical replicates. The 2-ΔΔCT method was used for qRT-PCR data analysis, in which PCT. (test) = CT (target, test) - CT (reference, check); ΔCT (calibrator) = CT (target, calibrator) - CT (reference, calibrator); and ΔΔCT = ΔCT (test) - ΔCT (calibrator). The expression ratio of the experimental and control groups was 2-ΔΔCT.
calmodulin and calmodulin-like
calcium-dependent protein kinase
cyclic nucleotide-gated channel
Clusters of Orthologous Groups
differentially expressed genes
ethylene response factor
false discovery rate
Kyoto Encyclopedia of Genes and Genomes
LRR receptor-like serine/threonine-protein kinases
mitogen-activated protein kinase
pathogen-associated molecular patterns
respiratory burst oxidase homologs
Reactive oxygen species
Mizrahi Y, Nerd A. Climbing and Columnar Cacti: New Arid Land Fruit Crops. Alexandria: Perspectives on New Crops & New Uses; 2014.
Montoya-Arroyo A, Schweiggert RM, Pineda-Castro M-L, Sramek M, Kohlus R, Carle R, Esquivel P. Characterization of cell wall polysaccharides of purple pitaya (Hylocereus sp.) pericarp. Food Hydrocoll. 2014;35:557–64.
Mizrahi Y, Nerd A, Nobel PS. Cacti as crops. Hortic Rev. 1997:291–319.
Adnan L, Osman A, Abdul Hamid A. Antioxidant activity of different extracts of red pitaya (Hylocereus polyrhizus) seed. Int J Food Prop. 2011;14(6):1171–81.
Chuang MF, Ni HF, Yang HR, Shu SL, Lai SY, Jiang YL. First report of stem canker disease of pitaya (Hylocereus undatus and H. polyrhizus) caused by Neoscytalidium dimidiatum in Taiwan. Plant Dis. 2012;96(6):906.
Lan GB, He ZF, Xi PG, Jiang ZD. First report of Brown spot disease caused by Neoscytalidium dimidiatum on Hylocereus undatus in Guangdong, Chinese Mainland. Plant Dis. 2012;96(11):1702.
Mohd MH, Salleh B, Zakaria L. Identification and molecular characterizations of Neoscytalidium dimidiatum causing stem canker of red-fleshed dragon fruit (Hylocereus polyrhizus) in Malaysia. J Phytopathol. 2013;161(11–12):841–9.
Ezra D, Liarzi O, Gat T, Hershcovich M, Dudai M. First report of internal black rot caused by Neoscytalidium dimidiatum on Hylocereus undatus (pitahaya) fruit in Israel. Plant Dis. 2013;97(11):1513.
Solsona GS, Lopez P, Palmateer AJ. First report of Neoscytalidium dimidiatum causing stem and fruit canker of Hylocereus undatus in Florida. Plant Dis. 2016;100(7).
Yi RH, Ling Lin Q, Mo JJ, Wu FF, Chen J. Fruit internal brown rot caused by Neoscytalidium dimidiatum on pitahaya in Guangdong province, China. Australas Plant Dis Notes. 2015;10(1):13.
Macho AP, Zipfel C. Plant PRRs and the activation of innate immune signaling. Mol Cell. 2014;54(2):263–72.
Gohre V, Spallek T, Haweker H, Mersmann S, Mentzel T, Boller T, de Torres M, Mansfield JW, Robatzek S. Plant pattern-recognition receptor FLS2 is directed for degradation by the bacterial ubiquitin ligase AvrPtoB. Curr Biol. 2008;18(23):1824–32.
Thomma BP, Nurnberger T, Joosten MH. Of PAMPs and effectors: the blurred PTI-ETI dichotomy. Plant Cell. 2011;23(1):4–15.
Yasuda S, Okada K, Saijo Y. A look at plant immunity through the window of the multitasking coreceptor BAK1. Curr Opin Plant Biol. 2017;38:10–8.
Boller T, Felix G. A renaissance of elicitors: perception of microbe-associated molecular patterns and danger signals by pattern-recognition receptors. Annu Rev Plant Biol. 2009;60:379–406.
Baggs E, Dagdas G, Krasileva KV. NLR diversity, helpers and integrated domains: making sense of the NLR IDentity. Curr Opin Plant Biol. 2017;38:59–67.
Peng Y, Wersch RV, Zhang Y. Convergent and divergent signaling in PAMP-triggered immunity and effector-triggered immunity. Mol Plant-Microbe Interact. 2017;31(04):403-409.
Qingzhu H, Chengjie C, Zhe C, Pengkun C, Yuewen M, Jingyu W, Jian Z, Guibing H, Jietang Z, Yonghua Q. Transcriptomic analysis reveals key genes related to Betalain biosynthesis in pulp coloration of Hylocereus polyrhizus. Front Plant Sci. 2015;6:1179.
Chiasson DM, Haage K, Sollweck K, Brachmann A, Dietrich P, Parniske M. A quantitative hypermorphic CNGC allele confers ectopic calcium flux and impairs cellular development. Elife. 2017;6.
Nawaz Z, Kakar KU, Saand MA, Shu QY. Cyclic nucleotide-gated ion channel gene family in rice, identification, characterization and experimental analysis of expression response to plant hormones, biotic and abiotic stresses. BMC Genomics. 2014;15:853.
Moeder W, Urquhart W, Ung H, Yoshioka K. The role of cyclic nucleotide-gated ion channels in plant immunity. Mol Plant. 2011;4(3):442–52.
Fischer C, DeFalco TA, Karia P, Snedden WA, Moeder W, Yoshioka K, Dietrich P. Calmodulin as a Ca2+-sensing subunit of Arabidopsis cyclic nucleotide-Gated Channel complexes. Plant Cell Physiol. 2017;58(7):1208–21.
Chen Y, Zhou X, Chang S, Chu Z, Wang H, Han S, Wang Y. Calcium-dependent protein kinase 21 phosphorylates 14-3-3 proteins in response to ABA signaling and salt stress in rice. Biochem Biophys Res Commun. 2017;493(4):1450–6.
Ding Y, Cao J, Ni L, Zhu Y, Zhang A, Tan M, Jiang M. ZmCPK11 is involved in abscisic acid-induced antioxidant defence and functions upstream of ZmMPK5 in abscisic acid signalling in maize. J Exp Bot. 2013;64(4):871–84.
Hettenhausen C, Yang DH, Baldwin IT, Wu J. Calcium-dependent protein kinases, CDPK4 and CDPK5, affect early steps of jasmonic acid biosynthesis in Nicotiana attenuata. Plant Signal Behav. 2013;8(1):e22784.
Pawelek A, Duszyn M, Swiezawska B, Szmidt-Jaworska A, Jaworski K. Transcriptional response of a novel HpCDPK1 kinase gene from Hippeastrum x hybr. to wounding and fungal infection. J Plant Physiol. 2017;216:108–17.
Wang JP, Xu YP, Munyampundu JP, Liu TY, Cai XZ. Calcium-dependent protein kinase (CDPK) and CDPK-related kinase (CRK) gene families in tomato: genome-wide identification and functional analyses in disease resistance. Mol Gen Genomics. 2016;291(2):661–76.
Weckwerth P, Ehlert B, Romeis T. ZmCPK1, a calcium-independent kinase member of the Zea mays CDPK gene family, functions as a negative regulator in cold stress signalling. Plant Cell Environ. 2015;38(3):544–58.
Vandelle E, Vannozzi A, Wong D, Danzi D, Digby AM, Dal Santo S, Astegno A. Identification, characterization, and expression analysis of calmodulin and calmodulin-like genes in grapevine (Vitis vinifera) reveal likely roles in stress responses. Plant Physiol Biochem. 2018;129:221–37.
Yamauchi T, Yoshioka M, Fukazawa A, Mori H, Nishizawa NK, Tsutsumi N, Yoshioka H, Nakazono M. An NADPH oxidase Rboh functions in Rice roots during Lysigenous Aerenchyma formation under oxygen-deficient conditions. Plant Cell. 2017;29(4):775–90.
Kaur G, Guruprasad K, Temple BRS, Shirvanyants DG, Dokholyan NV, Pati PK. Structural complexity and functional diversity of plant NADPH oxidases. Amino Acids. 2018;50(1):79–94.
Gong B, Yan Y, Wen D, Shi Q. Hydrogen peroxide produced by NADPH oxidase: a novel downstream signaling pathway in melatonin-induced stress tolerance in Solanum lycopersicum. Physiol Plant. 2017;160(4):396–409.
Li L, Li M, Yu L, Zhou Z, Liang X, Liu Z, Cai G, Gao L, Zhang X, Wang Y, et al. The FLS2-associated kinase BIK1 directly phosphorylates the NADPH oxidase RbohD to control plant immunity. Cell Host Microbe. 2014;15(3):329–38.
Chinchilla D, Bauer Z, Regenass M, Boller T, Felix G. The Arabidopsis receptor kinase FLS2 binds flg22 and determines the specificity of flagellin perception. Plant Cell. 2006;18(2):465–76.
Liang X, Ding P, Lian K, Wang J, Ma M, Li L, Li L, Li M, Zhang X, Chen S, et al. Arabidopsis heterotrimeric G proteins regulate immunity by directly coupling to the FLS2 receptor. Elife. 2016;5:e13568.
Sun J, Li L, Wang P, Zhang S, Wu J. Genome-wide characterization, evolution, and expression analysis of the leucine-rich repeat receptor-like protein kinase (LRR-RLK) gene family in Rosaceae genomes. BMC Genomics. 2017;18(1):763.
Bent AF, Kunkel BN, Dahlbeck D, Brown KL, Schmidt R, Giraudat J, Leung J, Staskawicz BJ. RPS2 of Arabidopsis thaliana: a leucine-rich repeat class of plant disease resistance genes. Science. 1994;265(5180):1856–60.
Swiderski MR, Innes RW. The Arabidopsis PBS1 resistance gene encodes a member of a novel protein kinase subfamily. Plant J. 2001;26(1):101–12.
Zhou J, Loh YT, Bressan RA, Martin GB. The tomato gene Pti1 encodes a serine/threonine kinase that is phosphorylated by Pto and is involved in the hypersensitive response. Cell. 1995;83(6):925–35.
Ron M, Avni A. The receptor for the fungal elicitor ethylene-inducing xylanase is a member of a resistance-like gene family in tomato. Plant Cell. 2004;16(6):1604–15.
Roy S. Function of MYB domain transcription factors in abiotic stress and epigenetic control of stress response in plant genome. Plant Signal Behav. 2016;11(1):e1117723.
Zhu Q, Droge-Laser W, Dixon RA, Lamb C. Transcriptional activation of plant defense genes. Curr Opin Genet Dev. 1996;6(5):624–30.
Stracke R, Werber M, Weisshaar B. The R2R3-MYB gene family in Arabidopsis thaliana. Curr Opin Plant Biol. 2001;4(5):447–56.
Chen L, Han J, Deng X, Tan S, Li L, Li L, Zhou J, Peng H, Yang G, He G, et al. Expansion and stress responses of AP2/EREBP superfamily in Brachypodium distachyon. Sci Rep. 2016;6:21623.
Li B, Li Q, Mao X, Li A, Wang J, Chang X, Hao C, Zhang X, Jing R. Two novel AP2/EREBP transcription factor genes TaPARG have pleiotropic functions on plant architecture and yield-related traits in common wheat. Front Plant Sci. 2016;7:1191.
Okamuro JK, Caster B, Villarroel R, Van Montagu M, Jofuku KD. The AP2 domain of APETALA2 defines a large new family of DNA binding proteins in Arabidopsis. Proc Natl Acad Sci U S A. 1997;94(13):7076–81.
Huang Y, Li T, Xu Z-S, Wang F, Xiong A-S. Six NAC transcription factors involved in response to TYLCV infection in resistant and susceptible tomato cultivars. Plant Physiol Biochem. 2017;120:61–74.
Chi Y, Yang Y, Zhou Y, Zhou J, Fan B, Yu JQ, Chen Z. Protein-protein interactions in the regulation of WRKY transcription factors. Mol Plant. 2013;6(2):287–300.
Katsuragi Y, Takai R, Furukawa T, Hirai H, Morimoto T, Katayama T, Murakami T, Che FS. CD2-1, the C-terminal region of Flagellin, modulates the induction of immune responses in Rice. Mol Plant-Microbe Interact. 2015;28(6):648–58.
Furukawa T, Inagaki H, Takai R, Hirai H, Che FS. Two distinct EF-Tu epitopes induce immune responses in rice and Arabidopsis. Mol Plant-Microbe Interact. 2014;27(2):113–24.
Boudsocq M, Sheen J. CDPKs in immune and stress signaling. Trends Plant Sci. 2013;18(1):30–40.
Zhang L, Du L, Poovaiah BW. Calcium signaling and biotic defense responses in plants. Plant Signal Behav. 2014;9(11):e973818.
Tunc-Ozdemir M, Tang C, Ishka MR, Brown E, Groves NR, Myers CT, Rato C, Poulsen LR, McDowell S, Miller G, et al. A cyclic nucleotide-gated channel (CNGC16) in pollen is critical for stress tolerance in pollen reproductive development. Plant Physiol. 2013;161(2):1010–20.
Sinha M, Singh RP, Kushwaha GS, Iqbal N, Singh A, Kaushik S, Kaur P, Sharma S, Singh TP. Current overview of allergens of plant pathogenesis related protein families. Sci World J. 2014;2014:543195.
Hwang IS, Choi DS, Kim NH, Kim DS, Hwang BK. Pathogenesis-related protein 4b interacts with leucine-rich repeat protein 1 to suppress PR4b-triggered cell death and defense response in pepper. Plant J. 2014;77(4):521–33.
Kim NH, Hwang BK. Pepper pathogenesis-related protein 4c is a plasma membrane-localized cysteine protease inhibitor that is required for plant cell death and defense signaling. Plant J. 2015;81(1):81–94.
Paz-Ares J, Ghosal D, Wienand U, Peterson PA, Saedler H. The regulatory c1 locus of Zea mays encodes a protein with homology to myb proto-oncogene products and with structural similarities to transcriptional activators. EMBO J. 1987;6(12):3553–8.
Baldoni E, Genga A, Cominelli E. Plant MYB transcription factors: their role in drought response mechanisms. Int J Mol Sci. 2015;16(7):15811–51.
Su LT, Li JW, Liu DQ, Zhai Y, Zhang HJ, Li XW, Zhang QL, Wang Y, Wang QY. A novel MYB transcription factor, GmMYBJ1, from soybean confers drought and cold tolerance in Arabidopsis thaliana. Gene. 2014;538(1):46–55.
Liu J, Osbourn A, Ma P. MYB transcription factors as regulators of Phenylpropanoid metabolism in plants. Mol Plant. 2015;8(5):689–708.
Chen X, Huang H, Qi T, Liu B, Song S. New perspective of the bHLH-MYB complex in jasmonate-regulated plant fertility in Arabidopsis. Plant Signal Behav. 2016;11(2):e1135280.
Wang Z, Tang J, Hu R, Wu P, Hou XL, Song XM, Xiong AS. Genome-wide analysis of the R2R3-MYB transcription factor genes in Chinese cabbage (Brassica rapa ssp. pekinensis) reveals their stress and hormone responsive patterns. BMC Genomics. 2015;16:17.
Cao L, Yu Y, Ding X, Zhu D, Yang F, Liu B, Sun X, Duan X, Yin K, Zhu Y. The Glycine soja NAC transcription factor GsNAC019 mediates the regulation of plant alkaline tolerance and ABA sensitivity. Plant Mol Biol. 2017;95(3):253–68.
Jiang H, Li H, Bu Q, Li C. The RHA2a-interacting proteins ANAC019 and ANAC055 may play a dual role in regulating ABA response and jasmonate response. Plant Signal Behav. 2009;4(5):464–6.
Bu Q, Jiang H, Li CB, Zhai Q, Zhang J, Wu X, Sun J, Xie Q, Li C. Role of the Arabidopsis thaliana NAC transcription factors ANAC019 and ANAC055 in regulating jasmonic acid-signaled defense responses. Cell Res. 2008;18(7):756–67.
Banerjee A, Roychoudhury A. WRKY proteins: signaling and regulation of expression during abiotic stress responses. Sci World J. 2015;2015:807560.
Dale J, James A, Paul JY, Khanna H, Smith M, Peraza-Echeverria S, Garcia-Bastidas F, Kema G, Waterhouse P, Mengersen K, et al. Transgenic Cavendish bananas with resistance to fusarium wilt tropical race 4. Nat Commun. 2017;8(1):1496.
Xu M, Peng Y, Qi Z, Yan Z, Yang L, He M-D, Li Q-X, Liu C-L, Ruan Y-Z, Wei S-S, et al. Identification of Neoscytalidium dimidiatum causing canker disease of pitaya in Hainan, China. Australas Plant Pathol. 2018;47(5):547–53.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.
Pertea G, Huang X, Liang F, Antonescu V, Sultana R, Karamycheva S, Lee Y, White J, Cheung F, Parvizi B, et al. TIGR gene indices clustering tools (TGICL): a software system for fast clustering of large EST datasets. Bioinformatics (Oxford, England). 2003;19(5):651–2.
Tarazona S, Garcia-Alcalde F, Dopazo J, Ferrer A, Conesa A. Differential expression in RNA-seq: a matter of depth. Genome Res. 2011;21(12):2213–23.
Audic S, Claverie JM. The significance of digital gene expression profiles. Genome Res. 1997;7(10):986–95.
Fengxia Y, Xiaopeng W, Guoli G, Jin T. Cloning and sequence analysis of housekeeping genes actin and UBQ from pitaya. Guizhou Agric Sci. 2013;41(9):4.
We would like to thank the Charlesworth Author Services for language editing services.
This work was supported by the Key Research and Development Project of Hainan Province (ZDYF2018080), the Crop Science Postgraduate Innovation Project of Hainan University Tropical Agriculture and Forestry College (ZWCX2018007), the Natural Science Foundation of China (NSFC; 31360364) and the Innovation and Research Project of the General Institutes of Higher Education of Hainan Province in 2018 (Hys2018–18).
Availability of data and materials
The datasets generated during the current study, in the form of raw sequenece reads from the four transcriptome samples of pitaya (Hylocereus polyrhizus) under normal and diseased conditions are available at the NCBI GEO database with the accession number GSE119976.
Ethics approval and consent to participate
The orchard of Ledong Country, Hainan Province clarified the permissions granted to collect H.polyrhizus plant materials for research in this article.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Xu, M., Liu, CL., Luo, J. et al. Transcriptomic de novo analysis of pitaya (Hylocereus polyrhizus) canker disease caused by Neoscytalidium dimidiatum. BMC Genomics 20, 10 (2019). https://doi.org/10.1186/s12864-018-5343-0