Transcriptome profiling and digital gene expression analysis of sweet potato for the identification of putative genes involved in the defense response against Fusarium oxysporum f. sp. batatas.

Sweet potato production is constrained by Fusarium wilt, which is caused by Fusarium oxysporum f. sp. batatas (Fob). The identification of genes related to disease resistance and the underlying mechanisms will contribute to improving disease resistance via sweet potato breeding programs. In the present study, we performed de novo transcriptome assembly and digital gene expression (DGE) profiling of sweet potato challenged with Fob using Illumina HiSeq technology. In total, 89,944,188 clean reads were generated from 12 samples and assembled into 101,988 unigenes with an average length of 666 bp; of these unigenes, 62,605 (61.38%) were functionally annotated in the NCBI non-redundant protein database by BLASTX with a cutoff E-value of 10-5. Clusters of Orthologous Groups (COG), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations were examined to explore the unigenes' functions. We constructed four DGE libraries for the sweet potato cultivars JinShan57 (JS57, highly resistant) and XinZhongHua (XZH, highly susceptible), which were challenged with pathogenic Fob. Genes that were differentially expressed in the four libraries were identified by comparing the transcriptomes. Various genes that were differentially expressed during defense, including chitin elicitor receptor kinase 1 (CERK), mitogen-activated protein kinase (MAPK), WRKY, NAC, MYB, and ethylene-responsive transcription factor (ERF), as well as resistance genes, pathogenesis-related genes, and genes involved in salicylic acid (SA) and jasmonic acid (JA) signaling pathways, were identified. These data represent a sequence resource for genetic and genomic studies of sweet potato that will enhance the understanding of the mechanism of disease resistance.

Introduction

Sweet potato is the 7th most important crop globally, with an annual area harvest of 8.0 million ha and a total global production of 104.5 million tons. Sweet potato is planted primarily in developing areas in Asia and Africa [1, 2]. In the past two decades, sweet potato production has decreased in Asia and greatly increased in Africa [2]. In the USA, the harvest area has increased from 33.3 kha in 2002 to 54.7 kha in 2014, and this increase has resulted in a corresponding increase in sweet potato production from 580.5 kilotons to 1.34 million tons [2].

Fusarium oxysporum, a soil-borne pathogenic fungus, is classified into more than 120 formae speciales and races, each of which infects one specific crop, such as sweet potato, tomato, banana, watermelon or pepper. The severity of the damage caused by these infections depends on the interaction between the pathogen and the host crop [3]. The production of sweet potato is constrained by Fusarium wilt, which is caused by Fusarium oxysporum f. sp. batatas (Fob) [4]. Sweet potato growth was limited in the Southeastern United States prior to the 1950s and South China prior to the 1980s, when resistant cultivars were bred and applied [5, 6]. Fusarium wilt has spread across many provinces in Central China and threatens sweet potato production [7].

During Fob infections, sweet potato secretes tyloses, cell wall material or brown matter into xylem vessels as a physiological response to impede Fob infection, leading to blockade of the vascular bundle [8, 9] and stem split and rot [5]. Colonization by Fob can induce the defense system in sweet potato. The susceptibility or resistance of sweet potato to Fusarium wilt depends on the interaction with Fob, consistent with the gene-for-gene hypothesis proposed by Flor [10] and the mode of co-evolution proposed by Jones [11]. The pathogen-associated molecular patterns (PAMPs), such as chitin, glucan, or glycoprotein, of fungal cells are recognized by pattern recognition receptors (PRRs) in plants, resulting in PAMP-triggered immunity (PTI) to prevent pathogen colonization. Successful pathogens release effector proteins to suppress PTI [11]. In response, plants have evolved resistance (R) proteins to recognize these effectors and trigger a secondary response, i.e., effector-triggered immunity (ETI), which induces the expression of downstream defense-related genes [12]. This model was elucidated based on the interaction between tomato and its pathogenic Fusarium oxysporum f. sp. lycopersici (Fol) [13]. Many R genes with identified functions have been isolated from Arabidopsis, rice, tomato, and other host plants [14-17].

Sweet potato is hexaploid with a complex genome that has not been sequenced [18]. High-throughput RNA sequencing based on next-generation sequencing technology (Illumina) has recently been applied in analyses of gene expression profiling and for the discovery of important genes involved in disease resistance and related pathways [19, 20]. Transcriptome sequencing offers a rapid approach to obtain expressed gene sequence information for sweet potato [21] and has been used to mine gene-based microsatellite markers [22, 23], elucidate the mechanism of development and metabolism in the purple tuberous roots of sweet potato [18], reveal gene expression patterns associated with tuberous root formation and development [24], and identify related genes involved in lignin and starch biosynthesis [21].

In the present study, transcriptome sequencing technology was used to perform RNA sequencing (RNA-Seq) of sweet potato cultivars inoculated with pathogenic F. oxysporum. We constructed a sufficiently large library using equally mixed total RNA. In total, 89,944,188 high-quality reads were retained and assembled into 62,605 unigenes. The unigenes were then annotated and assigned putative functions, classifications or pathways by alignment with public databases. The transcriptome data provide a global view of the gene expression in sweet potato during Fob infection, generating a resource for discovering candidate defense-related genes. Based on the transcriptome data, we constructed cDNA libraries and performed four digital gene expression (DGE) analyses using sequencing reads from the libraries to compare the gene expression profiles of a sweet potato cultivar challenged with pathogenic Fob and identify a set of genes involved in the defense response to Fob. The comparisons will provide a better understanding of the molecular mechanisms underlying the defense response of sweet potato to Fob.

Materials and methods

Plant and pathogen materials, treatment and RNA extraction

Two sweet potato cultivars were used as the plant materials: Jinshan57 (JS57, highly resistant to Fusarium wilt) and XinZhongHua (XZH, highly susceptible to Fusarium wilt). JS57 was used as a highly resistant control (disease index of 14.7), whereas XZH was used as a highly susceptible control (disease index of 100) to identify resistance to Fusarium wilt [4].

The virulent pathogen fungus Fusarium oxysporum f. sp. batatas07 (F07) was used. F07 was isolated from sweet potato with Fusarium wilt and stored at -20°C.

Seedlings of sweet potato were inoculated with F07 or water (as a control). Four samples were prepared as follows: (1) JS57 inoculated with F07 (JS57-F07); (2) JS57 in water (JS57-CK); (3) XZH inoculated with F07 (XZH-F07); and (4) XZH in water (XZH-CK).

The conidia solutions were prepared as follows: F07 was grown on potato sugar agar (PSA) medium in a Petri dish for 8 d at 28°C until the mycelia grew over the surface of the medium. The PSA medium with mycelia was cut into small pieces and transferred into 250-mL conical flasks containing 100 mL of liquid Czapek medium. Each flask containing eight pieces of PSA medium with mycelia was incubated at 28°C for 8 d on a rotary shaker. The liquid medium containing the conidia and mycelia was strained through cheesecloth, and the conidia were collected into clean empty flasks. The conidia concentration of the solution in the flask was determined using a hemocytometer. Prior to incubation, the conidia suspensions were diluted in sterilized distilled water to a concentration of 1 × 107 conidia/mL.

Healthy sweet potato seedlings that were field grown for at least 28 d and had no obvious disease symptoms were freshly cut to a length of 20 cm and cultivated in glass bottles containing conidial solutions of 1 × 107 conidia/mL. Seedlings were cultivated in water as a control. The seedlings were cultivated in growth chambers at 28°C for 24 h, and the basic stems of seedlings with a length of 2 cm were collected as samples. Three independent biological replicates of each treatment were collected and immediately frozen in liquid nitrogen. Total RNA was isolated from each sample using the RNAprepPure Plant Kit (Tianen Biotech, Beijing, China), according to the manufacturer's instructions. RNA degradation was detected by electrophoresis on 1% agarose gels, and the RNA purity was measured by a NanoPhotometer spectrophotometer (IMPLEN, CA, USA). We used the Qubit RNA Assay Kit with a Qubit2.0 Fluorometer (Life Technologies, CA, USA) to determine the RNA concentration. The RNA integrity was confirmed using an RNA Nano 6000 Assay Kit on the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library preparation and transcriptome sequencing

The sequencing library was prepared using 3 μg of RNA from each sample and the NEBNextUltra™ RNA Library Prep Kit for Illumina (NEB, USA) according to the manual. Index codes were added to sequences from each sample. In brief, poly-T oligo-attached magnetic beads were applied to purify mRNA from total RNA. Then, fragmentation was conducted in the presence of divalent cations and at higher temperature in NEBNext, First Strand Synthesis Reaction Buffer (5X). We utilized random hexamer primers and M-MuLV Reverse Transcriptase to perform first-strand cDNA synthesis and DNA Polymerase I and RNase H for the second-strand cDNA synthesis. The overhangs were filled to blunt ends by exonuclease/polymerase activity. Then, the 3’ ends of the DNA fragments were adenylated, and NEBNext, Adaptors with hairpin loop structures were ligated to the DNA fragments. cDNA fragments with a size of 150–200 bp were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). The cDNA was then incubated with the USER Enzyme (NEB, USA) at 37°C for 15 min and 95°C for 5min. PCR was performed using Phusion High-Fidelity DNA polymerase with Universal PCR primers and Index (X) Primer. The PCR products were purified with an AMPure XP system. The library quality was evaluated with the Agilent Bioanalyzer 2100 system. The index-coded samples were clustered using the cBot Cluster Generation System with the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manual. The library preparations were then sequenced on an Illumina HiSeq 2500 platform. The resultant datasets are available in the NCBI Gene Expression Omnibus repository under GEO Series accession number GSE89290 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE89290).

Quality control and de novo assembly

The data produced using the Illumina HiSeq 2500 platform were transferred into raw reads (raw data) using base calling. The raw reads were initially processed via in-house Perl scripts [25], and the reads containing adapters or N (indefinite base) higher than 10%as well as low-quality reads were then removed. The remaining high-quality clean reads were used to determine the error rate, Q20, Q30, GC content and sequence duplication level. The clean data were also used to perform the following bioinformatics analyses. According to Haas et al. [25], the clean reads were concatenated into two files. The left files (read1 files) from all libraries/samples were pooled into a large left.fq file, and the right files (read2 files) were pooled into a large right.fq file. The transcriptome assembly was obtained based on the left.fq and right.fq using Trinity [26] for de novo assembly, with a minimum kmer coverage of two and default setting for all other parameters.

Gene functional annotation

The functions of the assembled transcripts were annotated by searching their sequences against the following databases:

NR (NCBI non-redundant protein sequences, http://www.ncbi.nlm.nih.gov) with an E-value cutoff of 1e-5;

Nt (NCBI non-redundant nucleotide sequences, http://www.ncbi.nlm.nih.gov) with an E-value cutoff of 1e-5;

PFAM (Protein family, http://pfam.sanger.ac.uk/) with an E-value cutoff of 1e-2;

KOG/COG (Clusters of Orthologous Groups of proteins, http://www.ncbi.nlm.nih.gov/COG/) with an E-value cutoff of 1e-3; and

Swiss-Prot (A manually annotated and reviewed protein sequence database, http://www.ebi.ac.uk/uniprot/) with an E-value cutoff of 1e-5.

Based on the annotations in NR and PFAM, Blast2GO (v2.5) was used to obtain the GO annotations (http://www.geneontology.org) according to the molecular function, biological process and cellular component ontologies [27]. The Automatic Annotation Server and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (http://www.genome.jp/kegg) were used to generate pathway assignments with an E-value cutoff of 1e-5.

Differential expression analysis

The analysis of differential expression between the four groups was performed with the DESeq R package (1.10.1). The resulting P-values were adjusted according to Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P-value <0.05, as determined by DESeq, were regarded as differentially expressed.

Experimental validation

Five genes were examined using real-time quantitative PCR (Q-PCR) to confirm their expression levels using FPKM, and the primers (Table 1) were designed using software Primer 5 and synthesized by Invitrogen. The primer sequences are listed in Table 1. A total of 1 μg of total RNA was used to synthesize the cDNA using the Prime Script RT Reagent Kit (Takara, DaLian, China) in a 20-μL reaction mixture. qRT-PCR was performed on a StepOnePlus™ Real-Time PCR System (Thermo Fisher Scientific), and the sweet potato GAPDH gene [28] was used as the endogenous reference control. All reference and selected genes were measured in triplicate.

Table 1

Sequences of primers used in the qRT-PCR analysis.

Sample	Gene ID	Primers	Primer sequences
JS57	C52297	F	CCCTTCATTCCCTTTCCACC
		R	CGAGTTCTTCGTCGCCGTTA
	C51550	F	TTTGTTTCTTCCACCTCTGT
		R	ATGCACGAATATCGCCTTAG
	C47391	F	TGTTGGGTCTGATAGTATTG
		R	CTCTTCCCTTTATTCTTCTG
XZH	C55420	F	CAGCCCAATAATCCAAGCAC
		R	GCATCACATCGGACCCTAAC
	C56184	F	TACCAATCAAAGGCTCACAA
		R	CCACCTAAAGGCACCAAGAC
Reference	GAPDH	F	GCAGGAACCCGGAAGAGATT
		R	GCAGCCTTGTCCTTGTCAGTG

Each reaction had a total volume of 20 μL and contained 2 μL of cDNA, 0.2 μL (10 μmol L-1) of gene-specific primers and 0.1 μL of CXR SYBR Green Master Mix. The real-time quantitative PCR program was 95°C for 10 min followed by 40 cycles of95°C for 15 s, 55°C for 20 s and 72°C for 30 s. The data were processed with Excel 2007(Microsoft, WA, USA).

Results

RNA-Seq and de novo transcriptome assembly

To comprehensively generate the transcriptome and obtain insights into the molecular mechanisms involved in the resistance of sweet potato to Fusarium wilt, the total RNA from the stems of the highly resistant cultivar JS57 and highly susceptible cultivar XZH was isolated. The RNAs were equally mixed to create an RNA pool for the cDNA library, and a database of the sweet potato stem transcriptome was then generated through Illumina HiSeqTM 2500 sequencing. The deep sequencing produced 91,366,006 sequence reads with a length of 100 bp each and an average GC content of 45.10%. After removing adaptors, primer sequences, poly-A tails and short and low-quality sequences, 89,944,188 (98.44%, 11.24 Gbp) high-quality reads with Q≥20 were retained and used for the assembly.

Due to the absence of a sequenced genome, all trimmed and cleaned reads were aligned and assembled using Trinity [29], yielding a set of 146,859 transcript sequences that were longer than 201 bp with an N50 length of 1371 bp and 101,988 unigenes with an average length of 666 bp and an N50 length of 1079 bp. The length of the assembled unigenes varied from 201 to 15,576 bp, and more than 36.6% of these unigenes were longer than 501 bp, consistent with previously reported transcriptome data from sweet potato (Tables 2 and 3, Fig 1) [24, 30, 31].

Table 2

Frequency of the length distribution of the assembled transcripts and unigenes.

Transcript length interval	200–500 bp	500–1 kb	1–2 kb	>2 kb	Total
Number of transcripts	78031	30741	25195	12892	146859
Number of unigenes	64695	18719	12437	6137	101988

Table 3

Length distribution of the assembled transcripts and unigenes.

	Min Length	Mean Length	Median Length	Max Length	N50	N90	Total Nucleotides
Transcripts	201	812	459	15576	1371	314	119187574
Unigenes	201	666	367	15576	1079	267	67942174

Fig 1

Length distribution of the assembled transcripts and unigenes in sweet potato.

The horizontal axis represents the length intervals of the transcripts and unigenes, and the vertical axis represents the number of transcripts and unigenes.

The lack of a reference sweet potato genome sequence is a barrier to identifying candidate genes involved in disease resistance. The overview of the gene expression profile of sweet potato challenged with Fob obtained in this study offers abundant sequence data for identifying candidate genes associated with defense against Fusarium wilt.

Functional annotation

To annotate the obtained unigenes, we performed a BLASTX search against the NR NCBI protein databases with a cutoff E-value of 10−5 based on sequence similarity. In total, 62,605 unigenes were detected (Table 4) that had good comparability with known gene sequences in at least one database, corresponding to approximately 61.38% of the total unigenes; 8,683 (8.51%) unigenes were annotated in all databases, including NR, NCBI non-redundant nucleotide sequences (NT), KO, Swiss-Prot, Protein Family (PFAM), Gene Ontology (GO) and Eukaryotic Ortholog Groups (KOG). As shown in Fig 2, the species that provided the best BLASTX matches (first hit) was Nicotiana tomentosiformis (19.1%), with more than 10,000 genes having high homology, followed by Solanum lycopersicum, Naegleria gruberi, Coffea canephora and Solanum tuberosum. However, more than half of the unigenes matched other species.

Table 4

Statistics of unigenesannotated in public databases.

Search item	Number of unigenes	Percentage(%)
Annotated in NR	56,362	55.26
Annotated in NT	32,361	31.73
Annotated in KO	19,578	19.19
Annotated in Swiss-Prot	38,491	37.74
Annotated in PFAM	38,592	37.83
Annotated in GO	40,918	40.12
Annotated in KOG	22,613	22.17
Annotated in all databases	8,683	8.51
Annotated in at least one database	62,605	61.38
Total unigenes	101,988	100

Fig 2

Species distribution of the BLASTX matches of the transcriptome unigenes.

This figure shows the species distribution of the unigene BLASTX matches against the NR protein database (cutoff E-value of E<10−5) and the proportions of each species.

GO and KOG classifications

GO classifications offer an ontology of defined terms representing the properties of gene products [18]. According to the results of a BLASTX search against the NR protein database, 40,918 genes (40.1%) were assigned to at least one GO term and were categorized into 46 classes. To better review the GO classifications involving these genes, each GO term was further clustered into its parent term. Ultimately, the classes were clustered into 45 functional groups and three major classifications (Fig 3), and the term “binding” (GO: 0005488) in molecular function represented the largest group (23,444 unigenes), followed by “cellular process” (GO: 0009987) in biological process (23,339 unigenes), “metabolic process” (GO: 0008152) in biological process (23,146 unigenes), and “catalytic” (GO: 0003824) in molecular function (19,536). The smallest groups included the terms “hormone secretion” (GO: 0046879) in biological process (2), “symplast” (GO: 0055044) in cellular component (5), “nucleoid” (GO: 0009295) in cellular component (7) and “synapse” (GO: 0045202) in cellular component (11) (Fig 3).

Fig 3

Histogram of GO classifications.

The percentages indicate the proportion of unigenes with the GO annotations.

A phylogenetic classification was performed using COG databases to classify the potential functions. In total, 22,613 genes were aligned to 26 functional classes (Fig 4). The two largest cluster groups were “general function prediction only” (3,374, 14.9%) and “posttranslational modification, protein turnover, chaperones” (3,299, 14.6%). Interestingly, 122 genes were aligned to the “defense mechanisms” cluster, which is helpful for understanding the mechanisms of resistance to Fusarium wilt in sweet potato.

Fig 4

Histogram of the COG classification in sweet potato.

The horizontal axis represents the names of the COGs, and the vertical axis represents the proportion of the group of the total number. In total, 22,613 sequences had COG classifications across 25 categories.

KEGG pathway mapping

The KEGG database is widely used as a reference database of pathway networks for the integration and interpretation of large-scale datasets generated using high-throughput sequencing technology and is helpful for obtaining a better understanding of the functions of genes that are up- or down-regulated in response to the Fusarium wilt pathogen [18]. A BLASTX search against the KEGG protein database was performed using a cutoff E-value of e-5. In total, 19,578 distinct gene sequences were assigned to 276 KEGG pathways. The pathways most represented by unique sequences included translation (2,560 genes), signal transduction (2,089 genes) and carbohydrate metabolism (1,963 genes) (Fig 5). The KEGG annotations reveal the specific enrichment of pathways and the mechanisms involved in the response to Fob infection in sweet potato.

Fig 5

Histogram of KEGG clusters in sweet potato.

Gene expression profile of sweet potato in response to Fob

The library was prepared using 3 μg of RNA from each sample, and sequencing was performed on an Illumina HiSeq 2500 platform. The number of unambiguous tags in each library was normalized to the expected number of fragments per kilobase of transcript per million base pairs sequenced (FPKM), and the differentially expressed genes were then compared between different libraries. The DESeq R package (1.10.1) was used to analyze the differential expression between two groups using samples with three biological replicates [32]. DESeq provides a statistical approach for determining differential expression of digital gene expression data using a model based on the negative binomial distribution. Genes with an adjusted P-value <0.05 were considered differentially expressed.

A comparison of the differentially expressed genes between the libraries revealed that many genes were up- or down-regulated in both sweet potato cultivars after inoculation with Fob. Notably, the number of up-regulated genes was higher than the number of down-regulated genes. In the JS57 cultivar, the expression of 6,454 genes was affected by F07 inoculation (JS57-F07 library); of these genes, 3,843 and 2,611 were up-regulated and down-regulated, respectively, compared with the control (JS57-CK library). In the XZH cultivar, the expression of 6,826 genes was affected by inoculation with F07 (XZH-F07 library); of these genes, 5,552 and 1,274 genes were up-regulated and down-regulated, respectively, compared with the control (XZH-CK).

Experimental validation by qRT-PCR

To confirm the reliability of the Illumina Solexa sequencing data and the differential DGE results, five unigenes were selected for quantitative RT-PCR assays: c47391-g1, c51550-g2, c52297-g1, c54420-g2 and c56184-g1. The results revealed the same expression tendency observed in the DGE analysis in both cultivars (Fig 6). qRT-PCR indicated that c52297-g1, c51550-g2, c54420-g2 and c56184-g1 were up-regulated by 4.88-, 2.30-, 3.85- and 2.12-fold after infection with Fob, respectively, whereas these were found to be up-regulation by 1.40-, 1.62-, 1.71- and 3.1-fold according to the DGE results, respectively. The unigene c47391-g1 was down-regulated 33.19-fold according to the qRT-PCR results but by 1.93-fold according to the DGE results.

Fig 6

qRT-PCR validation of the Fob-induced fold-changes detected using RNA-Seq.

Standard error bars are shown for the fold-changes determined by qRT-PCR.

Identification of related genes involved in the defense response

We identified a set of differentially expressed genes involved in the defense response by comparing the commonalities and differences in the libraries of the highly resistant and highly susceptible cultivars upon Fob infection. These genes, include chitin elicitor receptor kinase 1 (CERK), mitogen-activated protein kinase (MAPK), WRKY, NAC, MYB, and ethylene-responsive transcription factor (ERF), as well as resistance (R) genes, pathogenesis-related (PR) genes and genes involved in the salicylic acid (SA) and jasmonic acid (JA) signaling pathways, and the total number of defense-related genes and the differentially expressed genes are shown in Table 5.

Table 5

Number of defense-related and differentially expressed genes identified in the response to Fob infection.

Related genes	Gene numbers involved	Significant differential expression	Ratio of differentially expressed genes to total (%)
CERK 1	5	2	40.0
MAPK	66	6	9.1
JA signaling	14	2	14.3
SA signaling	25	4	16.0
ERF	40	14	35.0
WRKY	123	8	6.5
NAC	20	2	10.0
MYB	22	1	4.5
R genes	94	3	3.2
PR genes	57	7	12.3

CERK 1 genes

CERk1 is a PRR in host plants that plays a crucial role in chitin perception in the immune response against fungal pathogens [33]. We identified two differentially expressed CERK1 genes (Table 6): c60044-g1 (putative CERK1-like) and c62516-g1 (putative CERK1). c60044-g1 was up-regulated in both JS57-F07 and XZH-F07. Interestingly, c62516-g1 was down-regulated in the highly resistant cultivar JS57-F07, suggesting a complicated interaction between sweet potato and Fob.

Table 6

Differentially expressed genes involved in the defense response against Fob infection.

	Gene ID	Blast annotation	JS57-F07 vs JS57-CK	XZH-F07 vs XZH-CK
CERK	c60044-g1	CERK1 [A. thaliana]	3.2	1.4
	c62516-g1	CERK1-like [N. tomentosiformis]	-4.1
MAPK	c51382-g1	Wound-induced MAPK [N. glutinosa]		0.9
	c54598-g1	MAPK 20 isoform X1 [N. tomentosiformis]		0.7
	c58290-g1	MAPK [I. batatas]		1.3
	c58750-g5	MAPK 9 [N. tomentosiformis]	0.8	1.1
	c65078-g2	MAPK 6 isoform X1 [N. tomentosiformis]	-1.2
	c66053-g1	Receptor protein kinase [R. communis]	-1.8
JA signaling	c43470-g1	Jasmonic acid 2 [S. lycopersicum]	3.0
	c52478-g1	Jasmonate O-methyltransferase [B. rapa]		-0.9
SA signaling	c56489-g2	Salicylic acid-binding protein 2 [N. benthamiana]	-1.3
	c57992-g1	Salicylic acid-binding protein 2-like [N. tomentosiformis]	-3.2
	c58625-g1	UDP-glycosyltransferase 74F2-like [N. tomentosiformis]		-2.0
	c65549-g1	UDP-glycosyltransferase 74F2-like [N. tomentosiformis]	1.3	-0.7
ERF	c33996-g1	Ethylene-responsive-element-binding protein [S. lycopersicum]	2.6
	c45744-g1	Ethylene response factor 3 [S. lycopersicum]		0.8
	c51416-g2	Ethylene response factor 4 [S. lycopersicum]		1.2
	c52634-g1	Ethylene response factor 2 [N. tomentosiformis]		1.3
	c53488-g1	Ethylene response factor 13 [A. thaliana]	1.2	0.9
	c55705-g1	Ethylene response 2 [A. thaliana]	2.2
	c57791-g1	Ethylene-responsive factor [C. chinense]	3.6	0.6
	c58358-g1	Ethylene response factor 13 [C. arabica]	1.9
	c58425-g1	Ethylene receptor homolog [N. tabacum]	2.0
	c60985-g1	Ethylene-responsive transcription factor ERF027 [A. thaliana]		1.6
	c61066-g1	Ethylene-responsive-element-binding factor 5 [P. hybrida]		-1.0
	c61066-g2	Ethylene response factor 14 [A. deliciosa]	2.6
	c64410-g1	Ethylene response factor 1 [I. batatas]		0.7
	c65518-g1	Ethylene-insensitive protein 3 [N. tabacum]	0.7	0.5
WRKY	c54420-g2	WRKY61 [N. tomentosiformis]		1.7
	c56135-g1	WRKY75 [N. tomentosiformis]		1.6
	c57247-g1	WRKY22 [C. roseus]	1.8	1.6
	c60589-g1	WRKY71 [A. thaliana]		1.3
	c64728-g1	WRKY6 [A. thaliana]		1.1
	c46250-g1	WRKY56 [A. thaliana]		-2.1
	c50623-g1	WRKY48-like [S. tuberosum]	-1.4
	c33599-g2	Protein containing a WRKY-binding domain	3.4
NAC	c51550-g2	NAC29 [A. thaliana]		2.1
	c52292-g3	NAC7-like [M. acuminata]	-3.7
MYB	c5221-g1	Ipomoea batatas IbMYB1-2b	4.6
PR protein	c35430-g1	PR-1 type-like [P. dactylifera]		2.1
	c40413-g1	PR class 4 [Prunus dulcis]	3.5
	c42566-g1	PR-4 [P. mume]	3.2
	c47344-g1	PR-4A-like isoform X1 [C. arietinum]	-1.5
	c54564-g1	PR-1 [N. tomentosiformis]	-1.5
	c55964-g1	PR-10 [N. benthamiana]	5.9	2.6
	c58341-g1	PR-1 type-like [P. dactylifera]	0.8
R gene	c56584-g1	NBS-coding resistance protein [A. thaliana]	9.8
	c66373-g2	Late blight resistance protein homolog R1B-14 [N. tomentosiformis]	-1.8
	c56383-g1	Resistance gene analog [I. batatas]		2.4

Note: Significant differences (FDR≤0.001) in relative levels and |log2.Fold-change|≥1 are shown in boldface. FDR, false discovery rate. log2.Fold-change, log2-fold-changes using the ratio base 2 logarithm.INF, the adjusted read count of the gene was 0 in the control sample.

MAPK signaling genes

MAPK signaling cascades comprise an MAP kinase kinase kinase (MPKKK), an MAP kinase kinase (MPKK or MAP2K), and an MAP kinase(MPK) that are transduced from the receptor to downstream targets as a signal to sequentially activate the expression of defense-related genes [34]. We identified 66 MAPK-related genes, of which six were differentially expressed. MAPKs have attracted much attention because these genes are involved in the induction and coordination of hormone-dependent and hormone-independent plant responses against necrotrophic fungal pathogens [34].

Genes associated with the SA and JA signaling pathways

The plant hormones SA and JA are major defensive players in the regulation of signaling networks underlying resistance to pathogens [35]. SA plays an extensive signaling role in plants, particularly in the defense against biotrophic/hemibiotrophic pathogens, by inducing cell death, activating the expression of defense-related genes, and provoking systemic acquired resistance (SAR) [36, 37]. We identified 20 genes that were associated with the SA pathway (Table 5), and four of these genes were differentially expressed. Two unigenes, c21763-g1 and c39777-g1, were up-regulated in XZH-F07, and both were annotated as S-adenosyl-L-methioninesalicylic acid carboxyl methyltransferase (SAMT), which has been implicated in disease resistance responses in rice [38]. The unigene c58625-g1 was down-regulated in XZH-F07 and annotated as UDP-glycosyltransferase (UDPG), which might influence the accumulation of free SA, leading to changes in the resistance of the host plant [39].

JA-dependent signaling is effective against necrotrophic pathogens [35], and we identified two differentially expressed JA metabolism-related genes that were annotated as jasmonic acid 2 and jasmonate O-methyltransferase. These genes might play an important role in defense against the necrotrophic phase of Fob infection.

WRKY transcription factor genes

In total, 123 WRKY genes were identified. Of these 123 WRKY genes, ten (8.1%) were significantly differentially expressed in at least one comparison (Tables 4 and 5). The WRKY transcription factors belong to a large family of transcriptional regulators that are almost exclusively expressed in plants. These transcription factors have diverse biological functions and play important roles in the defense response and hormone-controlled processes by acting as transcriptional activators or repressors that regulate major changes in plants during the early phases of root colonization with fungi [40-42].

Five unigenes (c54420-g2, c56135-g1, c57247-g1, c60589-g1, and c64728-g1) were exclusively up-regulated in XZH-F07 and were identified as putative WRKY61, WRKY75, WRKY22, WRKY71 and WRKY6, respectively, all of which are related to the defense response. In Arabidopsis, the WRKY61 gene plays a role in reducing the symptom severity of turnip crinkle virus disease and is associated with plant immunity pathways [43]. WRKY 75 plays a role in decreasing the severity of bacterial soft rot disease and activates a subset of defense-related genes in Chinese cabbage, similar to Arabidopsis [44]. WRKY22, WRKY71 and WRKY6 are all involved in rice defense responses. WRKY22 also plays a role in the resistance response to blast [45], and WRKY 71 functions as a transcriptional regulator upstream of OsNPR1 and OsPR1b in defense signaling pathways [46]. WRKY6 positively regulates defense responses by increasing SA concentrations through activation of the SA biosynthesis-related gene OsICS1 [47]. The unigene c46250-g1 was annotated as putative WRKY56 and was exclusively down-regulated in XZH-F07; however, the function of the WRKY56 gene remains unclear. The unigene c50623-g1 was annotated as a putative WRKY48-like gene, and its expression level was significantly lower in JS57-CK than in XZH-CK. WRKY48 was identified as a negative regulator of PR gene expression and basal resistance against the bacterial pathogen P. syringae in Arabidopsis [48]. The unigene c33599-g2 was annotated as a putative WRKY gene because its coding protein contains a WRKY-binding domain, and the expression level of WRKY in JS57-F07 was 3.4-fold higher than that in ZXH-F07, indicating an important role in the defense response.

NAC transcription factor genes

NAC transcription factors are a large family of transcriptional regulators that play crucial roles in plant immune responses, basal defense, and SAR [49]. We identified 20 NAC transcription factor genes, and only two (10%) of these genes were significantly differentially expressed in at least one comparison (Table 5). The unigene c51550-g2 was annotated AS NAC29, and its expression level increased in XZH-F07. NAC29 plays important roles in senescence processes and responses to salt and drought stresses in wheat [50]; however, there are no studies of the function of this gene in disease resistance. The other unigene, c52292-g3, was annotated as a NAC7-like gene, and its expression level was 3.7-fold lower in JS57-CK than in XZH-CK, indicating that this gene is also a negative regulator of the defense response.

ERF genes

ERFs are responsible for generating tolerance to stress in plants. Most ERFs are transcriptional activators, although some act as transcriptional repressors [51]. We identified 40 ERFs, of which 14 were differentially expressed. These were all up-regulated, except forc61066-g1 (putative ERF 5), indicating that most ERFs are activated in response to Fob invasion. The gene c61066-g1 was annotated as ERF5 and was down-regulated in XZH-F07, indicating that this gene may be a negative regulator of disease resistance in sweet potato. This finding is consistent with results in Arabidopsis indicating that ERF5 negatively regulates chitin signaling and plant defense against the fungal pathogen Alternaria brassicicola [52].

MYB transcription factor genes

The transcription factor MYB family of proteins is large and functionally diverse in all eukaryotes and plays crucial roles in the interaction of regulatory networks that control development, metabolism and the response to biotic or abiotic stresses [53]. The unigene c5221-g1 was annotated as IbMYB1 and was expressed at a 4.6-fold higher level in Js57-CK than in XZH-CK, indicating its important potential function in the defense response. In sweet potato, MYB1 (IbMYB1) stimulates the accumulation of anthocyanin pigmentation in the tuberous roots of sweet potato and the leaves of tobacco and arabidopsis [54-56]. IbMYB1 also improves tolerance and enhances secondary metabolism in transgenic potatoes [57].

PR genes

PR proteins are induced in plants upon attack by various pathogens and stresses and perform diverse functions, including plant defense, disease resistance, and antifungal activity [58]. PR proteins play important roles in hypersensitive responses and most likely contribute to SAR [59]. We identified 57 PR genes, of which seven were differentially expressed in at least one comparison and were annotated as PR-1-, PR-4- and PR-10-type proteins (Table 5). PR-1-type proteins appear to be hallmarks of activation of the hypersensitive response and SA-dependent SAR pathways [60]. The PR-4 protein is one of the less well-studied of the 17 PR protein groups and is considered an endochitinase with chitinase activity [61] and a bifunctional enzyme with both RNase and DNase activities [62]. We identified three genes that were annotated as PR-4 type proteins, but further studies are required to better understand their functions. The gene c55964-g1 was annotated as PR-10 and was up-regulated 5.9-foldin JS57-F07, which is markedly higher than the 2.6-fold increase in expression observed in XZH-F07, indicating its important function in the defense against Fob in sweet potato and consistent with its function in peanut [63].

R genes

According to the gene-for-gene hypothesis, R proteins recognize and interact with effectors and trigger ETI [12]. We identified 94 R genes in the transcriptome, most of which contain a highly conserved domain characterized by an NBS-LRR structure (data not shown). However, only three of these genes were differentially expressed (Table 5). Interestingly, the unigene c56584-g1 was identified as a putative NBS-coding R protein that was up-regulated 9.8-fold in JS57-F07, indicating a crucial role in the defense response.

Discussion

Sweet potato is an important crop worldwide and represents a non-model plant with a large, complex genome [18]. The lack of a reference genome sequence is a barrier to the identification of candidate genes involved in disease resistance. The application of high-throughput sequencing technology can generate a large amount of transcriptome data, which is helpful for gaining an understanding of the mechanism of resistance and identifying related genes.

In the present study, we first characterized the transcriptome profile of sweet potato challenged with Fob. Illumina sequencing generated 89,944,188 clean, high-quality reads that were assembled into 101,988 unigenes with an average length of 666 bp, and 61.38% of these were functionally annotated in the NR protein database. Undoubtedly, the results of our study provide an overview of the gene expression profile of sweet potato challenged with Fob and offer abundant sequence data for the discovery of candidate genes involved in the resistance to Fusarium wilt. A large proportion of the unigenes (18,574 unigenes, 18.2%) are longer than 1 kb, facilitating the discovery of defense-related genes.

Many genes were up- and down-regulated in sweet potato when challenged with Fob. A Venn diagram illustrating the overlap in the differentially expressed genes in the comparisons of JS57-F07 with JS57-CK and XZH-F07 with XZH-CK is shown in Fig 7. There were more up-regulated than down-regulated genes, with 5,552 up-regulated and 1,274 down-regulated genes in XZH and 3,843 up-regulated and 2,611 down-regulated genes in JS57. In contrast, in banana infected with Foc TR4, the numbers of up-regulated and down-regulated genes were almost the same [64]. In the two comparisons, 2,438 overlapping differentially expressed up-regulated genes were identified, accounting for 43.91% (2438/5552) and 63.44% (2438/3843) in each comparison. However, only 290 differentially expressed down-regulated genes were detected, accounting for 11.10% (290/2611) and 22.76% (290/1274) in each comparison.

Fig 7

Venn diagram showing overlapping differentially expressed genes in the comparisons of JS57-F07 with JS57-CK and XZH-F07 with XZH-CK.

Plants thrive in certain environments and are continuously challenged by various forms of biotic stresses; thus, plants have evolved sophisticated mechanisms to sense biotic aggressors, activate complex phytohormones signaling networks, stimulate plant immune system responses, and regulate the expression of defense-related genes [34]. We identified a set of defense-related genes, including transcription factors, PR genes, R genes, and genes involved in the SA and JA signaling pathways, in the transcriptome.

PAMPs of fungal pathogens, such as chitin, chitosan and endopolygalacturonase, are perceived by PRRs in host plants and activate PTI [34]. CERK1 is an important PRR that plays a crucial role in chitin perception in Arabidopsis thaliana in the immune response against fungal pathogens [33]. We identified two CERK1 genes that were differentially expressed in the examined libraries. Intriguingly, c62516-g1 (putative CERK1-like) was down-regulated 4.1-fold in JS57-F07, indicating that effectors of Fob might suppress the expression of c62516-g1 during colonization, thereby restraining PTI in host plants. However, because JS57 is highly resistant to Fob, this plant must have R genes that recognize the effector and trigger ETI to activate defense responses.

The recognition of pathogens by PRRs in plants triggers MAPK cascades and activates SA, JA or ET signaling pathways, facilitating the mounting of defense responses against invading pathogens and inducing the expression of numerous defense-related genes [34]. SA plays a crucial role in plant-pathogen interactions by activating defense responses and is also involved in inducing SAR [65]. SA activates plant immune responses against biotrophic pathogens, whereas JA activates plant immune responses against necrotrophic pathogens [35]. The JA and SA signaling pathways might interact antagonistically [66] or positively [67]. The crosstalk between SA and JA is controlled by a novel function of cytosolic NPR1 [68] and WRKY transcription factors [69]. As a hemibiotrophic pathogen, the lifestyle of Fob includes an initial biotrophic phase (host cell remains alive) and a necrotrophic phase (host tissue is killed) [70]. The crosstalk between the JA and SA signaling pathways might be more complicated in sweet potato during Fob invasion, and further research is required to better understand the interactions during signal transduction.

Effectors enable pathogens to overcome PTI in host plants, and these factors can be recognized by the corresponding R genes of the host plants, leading to ETI [11]. R proteins play a pivotal role in inducing ETI responses. Most R genes encode NB-LRR proteins [11]. We identified 94 R genes in the transcriptome, and three of these genes were differentially expressed. Interestingly, c56584-g1 (putative NBS-coding resistance protein) was up-regulated 9.8-fold in JS57-F07, indicating its important role in defense responses against Fob. Cloning of the R gene from sweet potato and the effector gene from Fob and further studies of the interaction between these genes are required to elucidate the defense mechanism of sweet potato.

Conclusions

In the present study, we de novo constructed and characterized the transcriptomes of sweet potato challenged with Fob and identified 101,988 unigenes, generating a broad survey of genes involved in disease resistance. We performed DGE profiling analyses of sweet potato challenged with Fob and calculated the numbers of differentially expressed genes between libraries, identifying several up- and down-regulated genes. A set of candidate genes in sweet potato involved in the disease response to Fusarium wilt were identified, including WRKY and NAC transcription factors, R genes, PR genes and SA pathway-related genes. Some of these genes, such as c56584-g1 (putative R gene) and c55964-g1 (putative PR-10), were significantly up-regulated and implicated in the defense response against Fob. These data provide an understanding of the molecular mechanisms of disease resistance. This sequence resource will be valuable for genetic and genomic studies and will accelerate breeding programs for sweet potato with Fusarium wilt resistance.