http://orcid.org/0000-0003-3620-8295
Figures
Abstract
The pear is an important temperate fruit worldwide that is produced by a group of species in the genus Pyrus. Callery pear (Pyrus calleryana Decne) is characterized by high resistance to multiple diseases, good adaptability, and high ornamental value, and is therefore widely planted in pear orchards for edible fruit production or as stock. Plant pathogens are a major threat to pear yield. Black spot disease, caused by the filamentous fungus Alternaria alternata, is one of the most serious diseases in pear. Elucidation of resistant genes to black spot disease is extremely important for understanding the underlying mechanisms as well as for the development of resistant cultivars. In this study, high-throughput single-strand RNA-sequencing was used to compare the transcriptome profiles of Callery pear leaves before and after A. alternata incubation for 7 days. The analysis yielded 73.3 Gb of clean data that were mapped onto the reference genome of the Chinese pear, and differentially expressed gene(DEG)s were identified with |log2FC| ≥ 1. Functional annotation demonstrated that black spot disease promoted great changes in the overall metabolism, and enrichment analysis of gene ontology terms showed that most of them are closely linked to signalling network and photosynthesis. Specifically, the genes included mainly transcription factors and genes involved in calcium signalling and ethylene and jasmonate pathways. Eight members of the ethylene response factor transcription factor gene family Group IX, including ERF1, ERF7, and ERF105, were up-regulated to 2.03–3.37-fold compared with CK, suggesting their role in the defence response to pathogen infection. Additionally, multiple transcription factors involved in biotic stresses, such as NAC78, NAC2, MYB44, and bHLH28, were up-regulated. Furthermore, we identified 144 long non-coding (lnc)RNAs, providing new insight into the involvement of lncRNAs in the response to black spot disease. Our study provides valuable data on the molecular genetics and functional genomic mechanisms of resistance to black spot disease in Callery pear. A good understanding of the molecular response to this disease will allow the development of durable and environmentally friendly control strategies.
Citation: Kan J, Liu T, Ma N, Li H, Li X, Wang J, et al. (2017) Transcriptome analysis of Callery pear (Pyrus calleryana) reveals a comprehensive signalling network in response to Alternaria alternata. PLoS ONE 12(9): e0184988. https://doi.org/10.1371/journal.pone.0184988
Editor: Yun Zheng, Kunming University of Science and Technology, CHINA
Received: March 21, 2017; Accepted: September 5, 2017; Published: September 21, 2017
Copyright: © 2017 Kan et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All RNA-seq files are available from the NCBI database (accession number PRJNA393405).
Funding: This study was funded by the Jiangsu Agricultural Science and Technology Independent Innovation Fund of China (CX(15)1023) and the Outstanding Youth Fund of Jiangsu Province (BK20160016). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The pear is an important fruit grown in Asia, Europe, and all other temperate regions worldwide, produced by a group of species in the genus Pyrus. Pyrus belongs to the subfamily Spiraeoideae in the family Rosaceae, which contains multiple fruits and flowers like apple, strawberry, and rose. The pear is a fruit with a delicate, pleasant taste and smooth texture and thus has a high customer acceptance. Approximately 22 million tons of pear fruit are produced yearly worldwide, representing 4% of all fruit production (FAOSTAT, http://faostat.fao.org/). The Japanese pear (Pyrus pyrifolia Nakai), European pear (P. communis L.), and Chinese pear (P. bretschneideri Rehd. and P. ussuriensis Maxim.) are three major edible species commercially grown for fruit production. In China, other species such as P. hopeiensis, P. sinkiangensis, P. armeniacifolia, P. pseudopashia, P. xerophila, P. betulifolia Bunge, P. phaeocarpa Rehder, Callery pear (P. calleryana Decne), and P. pashia Buch are cultivated regionally on a small scale. These pears are also planted as stock for grafting to enhance the resistance and tolerance of scion pear to biotic and abiotic stresses.
Callery pear is characterized by high resistance to multiple diseases, high adaptability, and high ornamental value. As a species native to China, it is widely used as graft stock in south China. In the early 20th century, the Callery pear was introduced into the US as graft stock in an aim to promote resistance to fire blight disease in the European pear, which was used as scion. In the following decades, the Callery pear was also used as an ornamental plant. In total, 155 thousand Callery pear seedlings were sold in the US in 1998, worth 30 million dollars [1]. Notably, as the Callery pear is well adapted to the US climate, it is now becoming an invasive plant.
Plant pathogens are a major threat to economically important crops worldwide. Black spot disease (BSD), caused by the filamentous fungus Alternaria alternata Keissler (previously Alternaria kikuchiana Tanaka), is among the most serious diseases that decrease pear production. Pear BSD was first reported by Tanaka in 1933 [2] and was later identified in France [3]. Fungi in the genus Alternaria produce host-selective toxins, the structures of which have been determined in Japanese pear, strawberry, tangerine, apple, tomato, and rough lemon pathotypes [4,5]. The fungus mainly infects leaves, fruits, and new shoots, and usually causes pear necrosis and defoliation, shortens the productive period, and decreases the post-harvest storage period. In 1988, a mutant of Japanese pear resistant to BSD was selected from the susceptible ‘Nijisseiki’ by chronic gamma irradiation [6]. Susceptibility to BSD is controlled by the single dominant gene A, and susceptible cultivars are heterozygous. To date, the homozygous dominant (A/A) genotype has not been identified in Japanese pear.
Although the genomes and transcriptomes of the Chinese pear and the European pear are publicly available, genomic information of Callery pear is lacking. To date, only around 5200 protein sequences are searchable in the NCBI GenBank (Dec 1st, 2016). Xu et al. (2015) reported 78,695 unigenes in Callery pear identified by high-throughput RNA sequencing (RNA-seq). Comparison with the reference genome and gene models available on the pear genome website (http://peargenome.njau.edu.cn) identified 218167 SNPs, 23248 indels, and 18322 alternative splicing events. Although the resistance of Callery pear to BSD was identified decades ago, its molecular mechanism remains unclear. This study aimed to assemble mRNAs and lncRNAs with reference to a related pear genome in order to identify DEGs in response to A. alternata.
Materials and methods
Plant materials
The Callery Pear seeds were collected from the Pear Germplasm Resource Center, Institute of Horticulture, Jiangsu Academy of Agricultural Sciences. The original source of these seeds were the inbred offspring of collected wild race from Shengzhou City, Zhejiang Province, China (E120.58, N30.01) in 2008. Seeds were rinsed with water, coated with wet sand, and stored in a refrigerator at 4°C on November 4, 2015. On March 14, 2016, the seeds were germinated in 30-cm-diameter plastic pots, which contained a 1:1:1 mixture of peat moss, perlite, and vermiculite. The seedlings were grown in the greenhouse at 25°C with a 14-h light/10-h dark photoperiod under a photon flux of 300 μmol m-2 s-1. The seedlings were irrigated with ¼ Murashige and Skoog nutrient solution under routine management. Seedlings at the nine-leaf stage were used for experiments.
Inoculation of A. alternata
Seedlings were inoculated with spores of a virulent strain of A. alternata provided by Dr Youzhou Liu (Institute of Plant Protection, Jiangsu Academy of Agricultural Sciences, Nanjing, China). The isolate was cultured in potato dextrose broth without shaking at 25°C for 10 days. Mycelial mats were washed with tap water and maintained at 25°C in the dark. The spores formed were collected, suspended in distilled water, and diluted to approximately 5 × 105 spores/mL by using a haemocytometer. The spore suspension was sprayed with a glass atomizer onto two young leaves detached from each plant. The inoculated leaves were incubated in a moist chamber at 25°C for 48 h. For each seedling, 50 μL conidia were daubed on the leave surfaces, and control leaves of each plant were mock-inoculated with sterile water. After 7 days, 3 leaves from 3 seedlings were pooled. The leaves were weighed, rinsed in water to remove contaminants, quickly dried with paper towels, quickly frozen in liquid nitrogen, and stored at –70°C until RNA extraction.
RNA library construction and sequencing
From two replicate samples refers to before and after incubation, total RNA was extracted using the Total RNA Kit (Tiangen, Beijing, China) and digested with DNase I (Takara, Japan) after quality testing on the Agilent 2100 BioAnalyzer. The RNA was purified with Dynabeads Oligo dT (Life, USA), and 100 ng of purified mRNA was used to construct RNA libraries using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA) according to the manufacturer’s instructions. Qubit testing, 2% gel testing, and high-sensitivity DNA chip testing were used to guarantee library quality. A 10 ng RNA library was cluster generated in a cBot instrument using the TruSeq PE Cluster Kit (Illumina, USA) and sequenced on Illumina Hiseq™ 2500 instrument by Nanjing Huasequen Biotechnologies Co, Ltd., China. Raw data of each sample were used for further analysis.
Transcriptome assembly
Trimmomatic (v0.32) was used to trim low-quality reads in raw data to obtain clean reads [7]. Reads with ambiguous N and adapter sequences (forward adapter: GATCGGAAGAGCACACGTCTGAACTCCAG, reverse adapter: GATCGGAAGAGCGTCGTGTAGGGAAAGAG) were trimmed, and low-quality bases with a Phred score lower than 20 on each end of the read were eliminated. Additionally, reads of fewer than 50 nt were eliminated. The genomic sequences and gene structure annotation of Chinese pear (P. bretschneideri Rehd.) were used as a reference for read mapping and transcriptome assembly. Specifically, clean reads were mapped to the pear genome using Bowtie software (v2.1.0) [8], and gene expression was determined as RPKM (Reads Per Kilo bases per Million reads [9] using SAMtools (v0.1.19) [10] in reference to annotated gene models with a GFF format file.
Digital gene expression tag profiling
For each sample, the reads were mapped onto the pear gene models using Bowtie software (v2.1.0). Gene expression was compared pair-wise between all time points for each cultivar using the EBSeq software [11]. The log2 of fold change (log2FC), p-value and q-value were calculated for all genes, and only transcripts with |log2FC| ≥ 1 were considered to be differentially expressed. It is noteworthy that for very low RPKM values (<0.01), we considered these values as 0.01 in DEG analysis for convenience.
Annotation and classification of genes
Gene annotations were obtained by BLASTP searching against the NR, KOG, KEGG, and TrEMBL databases with an E-value cut-off of 10−5. The protein functional annotation of the hits with the highest sequence similarity was chosen to represent the potential function of the sequence.
Real-time quantitative (q)PCR validation
Ten genes responsive to BSD identified by RNA-seq were chosen for experimental validation by qPCR: NAC78 (Pbr040795), NAC2 (Pbr027956), MYB44 (Pbr015309 and Pbr022028), ERF1 (Pbr023899 and Pbr012024), ERF5 (Pbr023902), ERF105 (Pbr012685), bHLH28 (Pbr018411), and TCP15 (Pbr039105). The primers for these genes were designed by Beacon Designer 7 software. qPCR assays were performed with three biological and three technical replicates as previously reported. Total RNA was isolated from inoculated and control leaves using a TaKaRa MiniBEST Plant RNA Extraction Kit (9769), according to the manufacturer’s instructions. The RNA was used to synthesize first-strand cDNA using M-MLV reverse transcriptase (TaKaRa, Japan). The SYBR Premix Ex Taq II reagent (Takara, Japan) with SYBR Green I as the fluorescent dye was used for qPCR on an ABI 7300 real-time PCR system (Applied Biosystems, Foster City, CA, USA). Each reaction contained 10 μL 2× SYBR Premix Ex Taq II Reagent, 1.0 μL cDNA sample, and 500 nM gene-specific primers in a final volume of 20 μL. RNA levels were expressed relative to that of the actin gene (forward: GGAATGGTCAAGGCTGGGTT; reverse: CAAAGCATCTGTGAGGTCACG) as 2–ΔΔCt, where Ct is the cycle threshold, according to a previous workflow [12,13].
Results
Morphological changes in response to A. alternata and confirmation of disease establishment
Upon inoculation with 5 × 105 spores/mL, all pear seedlings showed typical symptoms of black spots on the leaves after 7 days of incubation. PCR on cDNA generated from RNA extracted from inoculated Caller pear leaves with primers to detect the pathogen yielded an approximately 400-bp fragment, which was consistent with the 398-bp fragment in previous reports [14]. In addition, qPCR results showed that the expression of PcNAC1 increased 1.86-fold after incubation, confirming successful disease establishment (data not shown). Thus, a time period of 7 days after incubation was selected for characterization of gene expression to explore the possible molecular mechanism of the resistance response to BSD.
Sequencing and transcriptome assembly
To obtain a global overview of the Callery pear response to BSD, we constructed 4 RNA libraries containing mRNAs and lncRNAs, two replicates of which were constructed before incubation and two replicates were constructed after incubation. RNA sequencing produced 512.83 million 150-bp paired-end reads in total. After adapter removal and filtering steps, we obtained 73.3 Gb of clean data with more than 95.59% paired-end reads having a Phred score higher than Q30. All RNA-seq files are available from the NCBI database (accession number PRJNA393405).
High-quality filtered reads were used for reference-guided transcriptome assembly. The genome of the Chinese White pear was selected as the reference because it is closely related to Callery pear. The ratio of reads mapped to the Chinese pear genome was moderate, with a mapping rate between 61.54% and 65.06% (Table 1). This mapping ratio was close to that of 64.58% in a previous report [15]. For each sample, the mapped reads were assembled to genes according to the updated reference gene models, and the expression values were calculated for each gene and each isoform. The Chinese pear genome sequencing project predicted 42812 protein-coding genes, of which, on average, gene models consisted of transcript lengths of 2776 bp, coding lengths of 1172 bp, and 4.7 exons per gene. In this study, we obtained expression of 42812 genes (S1 Table). The average length of the assembled genes was longer than pervious study that adopted de novo assembly method.
Identification of lncRNAs in Callery pear
The distinction of protein-coding from non-coding transcripts relies on the homologous alignment with well-studied proteins. Mapped reads that overlapped with coding gene models of Chinese pear were eliminated; the rest of the mapped reads that contained non-coding representative characteristics of other species were identified as lncRNA candidates. A large number of predicted and validated plant lncRNAs can be searched in public databases, such as the NONCODE [16] and PNRD [17] databases. In addition, the CPC (Coding Potential Calculator), CNCI (Coding-Non-Coding Index), and PLEK (predictor of long non-coding RNAs and messenger RNAs based on an improved k-mer scheme) pipelines are software programs designed to identify lncRNAs via computational approaches. CPC is a support vector machine-based classifier that can discriminate coding from noncoding transcripts with high accuracy [18]. CNCI was developed to effectively classify protein-coding or non-coding transcripts by profiling adjoining nucleotide triplets independent of known annotations [19]. PLEK is an effective predictor of lncRNAs and mRNAs using an improved k-mer scheme and a support vector machine algorithm [20]. Approximately 1730 novel transcripts were identified to be potential lncRNA transcripts by the intersection results of CPC, CNCI and PLEK software (S2 Table). The lncRNAs were shorter than coding RNAs, with an average length of 564 bp. The lengths of 409 lncRNAs fell into the range of 200 bp to 300 bp, 343 lncRNAs were 300 bp to 400 bp, and 261 lncRNAs were 400 bp to 500 bp. Among the latter, only 27 lncRNAs were longer than 2000 bp (Fig 1).
DEGs in response to A. alternata
To explore the molecular response to BSD, DEG analysis was performed. For each time point in two genotypes of Callery pear, three biological replicates were included. The p-value and q-value were calculated to measure the significance and reliability to reduce bias. In this study, DEGs were defined as those genes having a fold-change > 2 after eliminating bias on the basis of the three biological replicates (S3 Table). In total, 2090 DEGs were identified in response to A. alternate (Fig 2).
Annotation of DEGs
The gene ontology (GO) consortium was set up to produce a dynamic, controlled vocabulary applicable to all eukaryotes even as the knowledge of gene and protein roles in cells is accumulating and changing [21]. In this study, GO enrichment analysis was carried out to clarify hierarchical relationships pertaining to molecular functions, cellular components, and biological processes of the 1191 DEGs in response to A. alternata.
To better understand the biological functions of genes, we conducted KEGG pathway enrichment analysis. Based on sequence homology, DEGs were classified into 34 functional groups in 4 categories (Fig 3). In the report by Xu et al. (2015) on the response of Callery pear to salt stress, around 2100 genes were classified into the Signal Transduction functional group of the Environmental Information Processing category based on KEGG analysis. This study identified around 900 DEGs in this category, indicating the response of Callery pear to BSD shares some common features with that to abiotic stress. In the Metabolism category, genes involved in Carbohydrate Metabolism, Energy Metabolism, and Amino Acid Metabolism were the most abundant. In Genetic Information Processing category, over 900 genes were assigned to Translation and around 700 genes were assigned to Folding, Sorting and Degradation.
Furthermore, KOG terms were assigned to Callery pear genes based on the best BLAST hits from the KOG database. Out of the 2090 DEGs in total, 1066 DEGs (51%) could be assigned KOG terms. The three most represented categories were KOG1187 (serine/threonine protein kinase, 19 DEGs), KOG0540 (3-methylcrotonyl-CoA carboxylase, non-biotin containing subunit/acetyl-CoA carboxylase carboxyl transferase, subunit beta, 15 DEGs) and KOG1192 (UDP-glucuronosyl and UDP-glucosyl transferase, 14 DEGs).
GO and KEGG enrichment analysis of DEGs
The above-mentioned GO annotation of DEGs was mainly focused on the accumulated number of DEGs in a distinct entry, ignoring the number of all genes in this species. Next, we performed KEGG enrichment analysis to unveil major pathways involved in the response to BSD based on the ratio of DEGs among all genes in each entry and on the p-value. In total, 1191 DEGs were assigned GO items and thus, were used in cluster frequency calculation. Correspondingly, 15348 were used as the denominator in the calculation as 15348 genes in Callery pear could be assigned GO items. The GO enrichment analysis showed that genes linked to electron transport, structural molecule activity, and cell surface were three most significant pathways in the response to A. alternata (Fig 4).
Similarly, 632 DEGs were annotated as functional enzymes involved in 23 pathways in the KEGG enrichment analysis (Fig 5). Interestingly, 23 out of 25 antenna protein (ko00196)-coding genes were differentially expressed, among 64 genes involved the photosynthesis pathway (ko00195), suggesting that photosynthesis is drastically deregulated in response to A. alternata. Carbon-related (ko01210, ko00190, and ko00710) and fatty acid metabolism-related (ko00061, ko01212, and ko00062) pathways were also differentially expressed in response to A. alternata. In addition, DEGs were also found in the flavonoid biosynthesis pathway.
LncRNAs involved in the response of Callery pear to A. alternata
lncRNAs regulate target genes via cis- or trans-regulatory effects. Differentially expressed lncRNAs were selected for target prediction using the same criteria as DEGs. The genes that could be paired with lncRNAs transcribed within a 20-kbp window upstream or downstream of lncRNAs were considered as cis target genes [22]. RNAplex software was used to screen trans-acting target genes [23]. We identified 144 differentially expressed lncRNAs.
Interestingly, GO enrichment analysis results of targets of differentially expressed lncRNAs revealed their strong correlation with the DEGs. Ten functional groups, including Structural molecule activity, Pathogenesis, Oxidoreductase activity, Metabolism, Lyase activity, Intracellular, Electron transport, Cytoplasm, Chromosome, Biosynthesis, were also found in the GO enrichment analysis of Callery pear DEGs in response to A. alternata (Fig 6).
qPCR validation
To validate RNA-seq data, 9 genes were selected for analysis by qPCR. The correlation between the results was evaluated in R language using the RNA-seq fold-change values and relative expression levels quantified by qPCR. The Pearson correlation coefficient (r = 0.88) indicated that the gene expression levels obtained by qPCR and RNA-seq correlated, confirming the accuracy and reproducibility of the RNA-seq results (Table 2, Fig 7).
Error bars represent standard errors from three independent biological replicates.
Discussion
Improvement of comparative transcriptome analysis of Callery pear
The current study had several innovations in experimental design and bioinformatics analysis as compared with the previous report on the transcriptome of Callery pear in response to salt stress by Xu et al. (2015). In the report by Xu et al., the mRNAs were first enriched and RNAs without poly(A) tails were eliminated by chemical methods before sequencing. In this study, only rRNAs and small RNAs were eliminated and mRNAs were not enriched, which means that a large amount of lncRNA fragments were retained in the RNA libraries. Moreover, Xu et al. assembled reads into genes by the de novo method. In contrast, as Callery pear is closely related to Chinese White pear, the genome of the latter was used for transcriptome assembly, resulting a more precise dataset of gene models.
Signalling in Callery pear in response to A. alternata
Five major categories of pathogenicity genes in A. alternata comprise signalling genes, toxin production- and detoxification-related genes, genes encoding metabolic enzymes, genes involved in the generation of specific infection structures, and genes encoding cell wall-degrading enzymes [24]. The signalling response of plants to stress involves a comprehensive series of molecular expression changes to achieve effective regulation of the response. First, receptor proteins in the cell membrane perceive environmental stimulation, inducing second messengers such as Ca2+, ROS, and phosphoinositide. These second messengers regulate the intracellular content of Ca2+, which bind intracellular Ca2+ receptors to activate the phosphorylation cascade. Some genes subsequently induced participate in plant hormone synthesis pathways and influence the contents of abscisic acid (ABA), salicylic acid (SA), jasmonic acid (JA), and ethylene (ET), which in turn activate downstream signalling pathways. In response to salt stress, Callery pear activates Ca2+ signalling pathways mediated by calcineurin B-like proteins and their interacting kinases [15]. Annotation of DEGs revealed many genes encoding second messengers, plant hormone synthetic proteins, and transcription factors. CML27 (also termed CaM27, Pbr038859) was up-regulated 2.88-fold in response to A. alternate by RNA-seq technology; this gene is also involved in calcium signalling pathways in Callery pear plants in response to salt stress [15].
Activation of ET biosynthesis in response to A. alternata
ET is an important plant hormone that is transported as its precursor, 1-aminoeyclopropane-l-carboxylic acid (ACC), throughout the plant. ACC synthase and ACC oxidase are two key enzymes that control ET biosynthesis in higher plants. ACC is produced from S-adenosylmethionine (SAM) by ACC synthase and is converted to ET by ACC oxidase (Arimura et al. 2005). SAM-dependent methyltransferases are crucial in metabolism, cell signalling, and epigenetic programming via methylation of small molecules [25]. The SAM-coding genes Pbr010088 and Pbr016813 were up-regulated 1.26- and 1.82-fold by RNA-seq technology, respectively.
Activation of JA and SA pathways in response to A. alternata
JAs are plant signalling molecules that are systemically activated in response to various biotic and abiotic stresses [26]. Activation of the JA pathway increases the resistance of host plants to some pathogens, while activation of the SA pathway mediates resistance to other pathogens [27]. In this study, dozens of genes involved in JA metabolism were up- or down-regulated. Two genes (Pbr004541 and Pbr004568) encoding LOX2, which is involved in JA biosynthesis upon wounding or methyl jasmonate treatment in Arabidopsis [28], were increased over 1.5-fold in Callery pear in response to A. alternata. Lipid peroxidation is common to all biological systems and appears in the interaction of plants with pathogens, insects, and abiotic stresses, and at distinct stages of development. As ubiquitous proteins in plants, linoleate:oxygen oxidoreductases (LOXs) contribute to the oxygenation of linoleic acid leading to (13S)-hydroperoxy derivatives of poly-unsaturated fatty acids [29].
Pbr015509, which encodes PR protein B, was up-regulated 1.11-fold in this study. JA can induce major signalling components in the JA pathways, including F-box proteins and PR proteins. Thionins are small (5 kDa), usually basic, cysteine-rich peptides that are classified into the PR-13 family of PR proteins [30]. The expression of leaf thionins increases upon fungal infection in barley [31] and A. thaliana [32]. We identified only a few genes involved in SA signalling, indicating that SA signalling plays a limited role in the response of Callery pear to BSD.
Transcription factors potentially involved in the resistance of Callery pear to BSD
In plants, transcription factors play a central role in synchronizing cell metabolism by self-regulation or by regulating the transcription of downstream target genes. Various transcription factors involved in biotic stress and plant development, including NAC, MYB, ERF, bHLH, and TCP, were found to be potentially involved in the response of Callery pear to A. alternata. We validated their expression by qPCR.
The ET-responsive transcription factor (ERF) family, which is part of the AP2/ERF superfamily, is a large gene family of transcription factors of about 60 to 70 amino acids in Arabidopsis [33]. Among over 100 members in Arabidopsis and rice, the ERF genes in group IX have often been associated with defence gene expression in response to pathogen infection [34]. Eight members of Group IX were up-regulated 1.03–2.37-fold in this study, including homologs of AtERF1 (Pbr023899, Pbr035787, Pbr012024, Pbr035775), AtERF5 (Pbr023902), AtERF9 (Pbr016222), and AtERF105 (Pbr012685), whereas no members were found to be down-regulated. ERF1 is associated with the molecular response and resistance to fungal infection. Overexpression of ERF1 in Arabidopsis enhanced resistance to two necrotrophic fungi, while it reduced resistance to Pseudomonas syringae [35]. ERF5 might act as a positive regulator of JA-mediated defence to Botrytis cinerea. Expression of ERF105 increases strongly in response to cold treatment as well as cold stress-related hormones in Arabidopsis [36]. Besides Group IX, Pbr024467, which is homologous to AtERF7 in Group VIII and plays an important role in the ABA response in plants [37], was up-regulated 1.44-fold in our study. Up-regulation of these genes was validated by qPCR, illustrating that ERF genes play an important role in the response of Callery pear to A. alternata and suggesting their involvement in the interaction networks between with ET, JA, and ABA signalling pathways.
MYB and MYB-like protein members were suggested to be involved in calcium signalling pathways in the salt-stress response in pear as indicated by high-throughput RNA-seq; however, their involvement was not verified [15]. In the current study, several proteins homologous to MYB4, MYB44, MYB308, MYB-like protein APL, MYB-like protein J, and MYB-like protein G were identified as being expressed in response to A. alternata. MYB4 (Pbr024492), which is orthologous to At1g22640, an SA- and JA-inducible transcription factor [38], was increased 16.61-fold. In rice, MYB4 represents a crucial node in the cross-talk of stress signalling cascades through the activation of multiple components, thus improving tolerance/resistance to multiple stress conditions such as drought, salt, UV, ozone, viruses, bacteria, and fungi [39]. In Arabidopsis, MYB44 belongs to the R2R3 MYB subgroup transcription factor family, which can be induced by ABA, methyl jasmonate, ET, and abiotic stresses such as dehydration, low temperature, and salinity [40]. qPCR results showed that two genes encoding MYB44 were up-regulated 0.90- and 1.47-fold as possible pathogen-responsive genes.
qPCR results confirmed that the NAC2 gene, which is mainly activated by salt stress in Arabidopsis [41], was up-regulated over 3.19-fold in Callery pear by A. alternata. According to qPCR results, NAC78 increased 0.82-fold; NAC78 has been implicated as a positive regulator of the expression of core 26S [42,43]. Moreover, it is involved in the regulation of gene expression related to flavonoid biosynthesis [44], which might explain why some genes involved in flavonoid biosynthesis were differentially expressed according to KEGG enrichment analysis.
Moreover, TEOSINTE BRANCHED1-CYCLOIDEA-PCF transcription factors (TCP), which play a vital role in the regulation of developmental processes and hormone responses, were differentially expressed in Callery pear upon BSD. TCP15 is a repressor of auxin biosynthesis, and regulates cytokinin and auxin responses during gynoecium development in Arabidopsis [45]. Additionally, TCP15 can modulate the expression of transcription factors involved in the induction of anthocyanin biosynthesis genes, thus inhibiting anthocyanin accumulation during exposure to high light intensity [46]. TCP15 (Pbr039105) was up-regulated 2.34-fold by RNA-seq and up-regulated 1.73-fold by qPCR. This suggested that TCP15 might act as positive regulator of the plant response to biotic stress.
Genes linked to pathogenicity factors
A large number of genes in plants are linked to pathogenicity factors of fungal infection. For example, the plant FRP1 gene is required for pathogenesis. The fungus Fusarium oxysporum completely loses pathogenicity in the absence of the FRP1 gene in tomato, and fully regains virulence upon re-introduction of this gene. F-box proteins interact with S-phase kinase-associated protein 1 (Skp1) and facilitate targeting of proteins to the SCF—ubiquitination complex. SCF mediates ubiquitination in the infectious programme of F. oxysporum [24]. Five genes (Pbr005806, Pbr023291, Pbr034809, Pbr007369, and Pbr007370) encoding SCF ubiquitin ligase were up-regulated 1.22–2.19 fold after A. alternata invasion. In addition, 9 genes coding F-box proteins were identified, 8 of which were up-regulated.
Conclusion
In conclusion, we performed a comprehensive analysis of the molecular response to A. alternata, the causative agent of BSD, using ssRNA-seq. In total, 1946 DEGs were identified. Additionally, 144 long non-coding RNAs were identified that might be involved in the response to A. alternata. This assembled and annotated transcriptome provides a valuable genomic resource for further investigation of the molecular mechanisms of resistance to BSD in Callery pear. A good understanding of molecular response to pathogens will allow the development of durable and environmentally friendly strategies for disease control.
Supporting information
S1 Table. The gene expression of Callery pear.
https://doi.org/10.1371/journal.pone.0184988.s001
(XLSX)
S2 Table. The expression of lncRNAs in Callery pear.
https://doi.org/10.1371/journal.pone.0184988.s002
(XLSX)
S3 Table. The expression and annotation of differentially expressed gene in response to fungus Alternaria alternate.
https://doi.org/10.1371/journal.pone.0184988.s003
(XLSX)
Acknowledgments
This study was funded by the Jiangsu Agricultural Science and Technology Independent Innovation Fund of China (CX(15)1023) and the Outstanding Youth Fund of Jiangsu Province (BK20160016). We thank the staff of Nanjing Huasequen Biotechnologies Co, Ltd., China for help with the bioinformatics analysis.
References
- 1. Li Y, Cheng ZM, Smith WA, Ellis DR, Chen YQ, Zheng XL, et al. (2004) Invasive ornamental plants: Problems, challenges, and molecular tools to neutralize their invasiveness. Critical Reviews in Plant Sciences 23: 381–389.
- 2. Simmons EG, Roberts RG (1993) Alternaria Themes and Variations (73). Mycotaxon 48: 109–140.
- 3. Baudry A, Morzieres J, Larue P (1993) First report of Japanese pear black spot caused by Alternaria kikuchiana in France. Plant Disease 77.
- 4. Tsuge T, Harimoto Y, Akimitsu K, Ohtani K, Kodama M, Akagi Y, et al. (2013) Host-selective toxins produced by the plant pathogenic fungus Alternaria alternata. Fems Microbiology Reviews 37: 44–66. pmid:22846083
- 5. Walton JD (1996) Host-selective toxins: Agents of compatibility. Plant Cell 8: 1723–1733. pmid:8914323
- 6. Sanada T, Nishida T, Ikeda F (1988) Resistant Mutant to Black Spot Disease of Japanese Pear Nijisseiki Induced by Gamma-Rays. Journal of the Japanese Society for Horticultural Science 57: 159–166.
- 7. Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30: 2114–2120. pmid:24695404
- 8. Langmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nature Methods 9: 357–U354. pmid:22388286
- 9. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B (2008) Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods 5: 621–628. pmid:18516045
- 10. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. (2009) The Sequence Alignment/Map format and SAMtools. Bioinformatics 25: 2078–2079. pmid:19505943
- 11. Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BMG, et al. (2013) EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments (vol 29, pg 1035, 2013). Bioinformatics 29: 2073–2073.
- 12. Pfaffl MW (2001) A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Research 29.
- 13. Tan HW, Song XM, Duan WK, Wang Y, Hou XL (2015) Genome-wide analysis of the SBP-box gene family in Chinese cabbage (Brassica rapa subsp. pekinensis). Genome 58: 463–477. pmid:26599708
- 14. Li Y, Chen X, Zhang A, Yang X, Chen Y, Yao J (2016) Molecular Detection of Alternaria alternata Causing Black Spot of Dangshan Pear. Chinese Agricultural Science Bulletin 32: 150–154.
- 15. Xu Y, Li X, Lin J, Wang Z, Yang Q, Chang Y (2015) Transcriptome sequencing and analysis of major genes involved in calcium signaling pathways in pear plants (Pyrus calleryana Decne.). BMC Genomics 16: 738. pmid:26424153
- 16. Bu DC, Yu KT, Sun SL, Xie CY, Skogerbo G, Miao RY, et al. (2012) NONCODE v3.0: integrative annotation of long noncoding RNAs. Nucleic Acids Research 40: D210–D215. pmid:22135294
- 17. Yi X, Zhang ZH, Ling Y, Xu WY, Su Z (2015) PNRD: a plant non-coding RNA database. Nucleic Acids Research 43: D982–D989. pmid:25398903
- 18. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, et al. (2007) CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Research 35: W345–W349. pmid:17631615
- 19. Sun L, Luo HT, Bu DC, Zhao GG, Yu KT, Zhang CH, et al. (2013) Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Research 41.
- 20. Li A, Zhang J, Zhou Z (2014) PLEK: a tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinformatics 15: 311. pmid:25239089
- 21. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. (2000) Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 25: 25–29. pmid:10802651
- 22. Han L, Zhang KL, Shi ZD, Zhang JX, Zhu JL, Zhu SJ, et al. (2012) LncRNA profile of glioblastoma reveals the potential role of lncRNAs in contributing to glioblastoma pathogenesis. International Journal of Oncology 40: 2004–2012. pmid:22446686
- 23. Tafer H, Hofacker IL (2008) RNAplex: a fast tool for RNARNA interaction search. Bioinformatics 24: 2657–2663. pmid:18434344
- 24. Duyvesteijn RGE, van Wijk R, Boer Y, Rep M, Cornelissen BJC, Haring MA (2005) Frp1 is a Fusarium oxysporum F-box protein required for pathogenicity on tomato. Molecular Microbiology 57: 1051–1063. pmid:16091043
- 25. Grove TL, Benner JS, Radle MI, Ahlum JH, Landgraf BJ, Krebs C, et al. (2011) A radically different mechanism for S-adenosylmethionine-dependent methyltransferases. Science 332: 604–607. pmid:21415317
- 26. Turner JG, Ellis C, Devoto A (2002) The jasmonate signal pathway. Plant Cell 14: S153–S164. pmid:12045275
- 27. Koornneef A, Pieterse CMJ (2008) Cross talk in defense signaling. Plant Physiology 146: 839–844. pmid:18316638
- 28. Bell E, Creelman RA, Mullet JE (1995) A Chloroplast Lipoxygenase Is Required for Wound-Induced Jasmonic Acid Accumulation in Arabidopsis. Proceedings of the National Academy of Sciences of the United States of America 92: 8675–8679. pmid:7567995
- 29. Feussner I, Wasternack C (2002) The lipoxygenase pathway. Annual Review of Plant Biology 53: 275–297. pmid:12221977
- 30. Sels J, Mathys J, De Coninck BMA, Cammue BPA, De Bolle MFC (2008) Plant pathogenesis-related (PR) proteins: A focus on PR peptides. Plant Physiology and Biochemistry 46: 941–950. pmid:18674922
- 31. Bohlmann H, Clausen S, Behnke S, Giese H, Hiller C, Reimannphilipp U, et al. (1988) Leaf-Specific Thionins of Barley—a Novel Class of Cell-Wall Proteins Toxic to Plant-Pathogenic Fungi and Possibly Involved in the Defense-Mechanism of Plants. Embo Journal 7: 1559–1565. pmid:16453847
- 32. Epple P, Apel K, Bohlmann H (1995) An Arabidopsis-Thaliana Thionin Gene Is Inducible Via a Signal-Transduction Pathway Different from That for Pathogenesis-Related Proteins. Plant Physiology 109: 813–820. pmid:8552715
- 33. Riechmann JL, Heard J, Martin G, Reuber L, Jiang C, Keddie J, et al. (2000) Arabidopsis transcription factors: genome-wide comparative analysis among eukaryotes. Science 290: 2105–2110. pmid:11118137
- 34. Nakano T, Suzuki K, Fujimura T, Shinshi H (2006) Genome-wide analysis of the ERF gene family in Arabidopsis and rice. Plant Physiology 140: 411–432. pmid:16407444
- 35. Berrocal-Lobo M, Molina A, Solano R (2002) Constitutive expression of ETHYLENE-RESPONSE-FACTOR1 in Arabidopsis confers resistance to several necrotrophic fungi. Plant Journal 29: 23–32. pmid:12060224
- 36. Bolt S, Zuther E, Zintl S, Hincha DK, Schmulling T (2017) ERF105 is a transcription factor gene of Arabidopsis thaliana required for freezing tolerance and cold acclimation. Plant Cell Environ 40: 108–120. pmid:27723941
- 37. Song CP, Agarwal M, Ohta M, Guo Y, Halfter U, Wang P, et al. (2005) Role of an Arabidopsis AP2/EREBP-type transcriptional repressor in abscisic acid and drought stress responses. Plant Cell 17: 2384–2396. pmid:15994908
- 38. Gil-Amado JA, Gomez-Jimenez MC (2013) Transcriptome analysis of mature fruit abscission control in olive. Plant Cell Physiol 54: 244–269. pmid:23292600
- 39. Vannini C, Iriti M, Bracale M, Locatelli F, Faoro F, Croce P, et al. (2006) The ectopic expression of the rice Osmyb4 gene in Arabidopsis increases tolerance to abiotic, environmental and biotic stresses. Physiological and Molecular Plant Pathology 69: 26–42.
- 40. Jung C, Seo JS, Han SW, Koo YJ, Kim CH, Song SI, et al. (2008) Overexpression of AtMYB44 enhances stomatal closure to confer abiotic stress tolerance in transgenic Arabidopsis. Plant Physiology 146: 623–635. pmid:18162593
- 41. Balazadeh S, Siddiqui H, Allu AD, Matallana-Ramirez LP, Caldana C, Mehrnia M, et al. (2010) A gene regulatory network controlled by the NAC transcription factor ANAC092/AtNAC2/ORE1 during salt-promoted senescence. Plant Journal 62: 250–264. pmid:20113437
- 42. Gladman NP, Marshall RS, Lee KH, Vierstra RD (2016) The Proteasome Stress Regulon Is Controlled by a Pair of NAC Transcription Factors in Arabidopsis. Plant Cell 28: 1279–1296. pmid:27194708
- 43. Yabuta Y, Osada R, Morishita T, Nishizawa-Yokoi A, Tamoi M, Maruta T, et al. (2011) Involvement of Arabidopsis NAC transcription factor in the regulation of 20S and 26S proteasomes. Plant Science 181: 421–427. pmid:21889048
- 44. Morishita T, Kojima Y, Maruta T, Nishizawa-Yokoi A, Yabuta Y, Shigeoka S (2009) Arabidopsis NAC Transcription Factor, ANAC078, Regulates Flavonoid Biosynthesis under High-light. Plant and Cell Physiology 50: 2210–2222. pmid:19887540
- 45. Lucero LE, Uberti-Manassero NG, Arce AL, Colombatti F, Alemano SG, Gonzalez DH (2015) TCP15 modulates cytokinin and auxin responses during gynoecium development in Arabidopsis. Plant J 84: 267–282. pmid:26303297
- 46. Viola IL, Camoirano A, Gonzalez DH (2016) Redox-Dependent Modulation of Anthocyanin Biosynthesis by the TCP Transcription Factor TCP15 during Exposure to High Light Intensity Conditions in Arabidopsis. Plant Physiol 170: 74–85. pmid:26574599