Transcriptome analysis reveals the molecular mechanisms of response to an emergent yellow-flower disease in green Chinese prickly ash (Zanthoxylum schinifolium)

Chinese prickly ash (Zanthoxylum) is extensively used as spice and traditional medicine in eastern Asian countries. Recently, an emergent yellow-flower disease (YFD) break out in green Chinese prickly ash (Zanthoxylum schinifolium, Qinghuajiao in Chinese) at Chongqing municipality, and then leads to a sharp reduction in the yield of Qinghuajiao, and thus results in great economic losses for farmers. To address the molecular response for the emergent YFD of Qinghuajiao, we analyzed the transcriptome of 12 samples including the leaves and inflorescences of asymptomatic and symptomatic plants from three different towns at Chongqing by high-throughput RNA-Seq technique. A total of 126,550 genes and 229,643 transcripts were obtained, and 21,054 unigenes were expressed in all 12 samples. There were 56 and 164 different expressed genes (DEGs) for the AL_vs_SL (asymptomatic leaf vs symptomatic leaf) and AF_vs_SF (asymptomatic flower vs symptomatic flower) groups, respectively. The results of KEGG analysis showed that the “phenylpropanoid biosynthesis” pathway that related to plant–pathogen interaction were found in AL_vs_SL and AF_vs_SF groups, and the “Plant–pathogen interaction” found in AF_vs_SF group, implying that this Qinghuajiao YFD might cause by plant pathogen. Interestingly, we detected 33 common unigenes for the 2 groups, and almost these unigenes were up-regulated in the symptomatic plants. Moreover, most of which were homologs to virus RNA, the components of viruses, implying that this YFD was related to virus. Our results provided a primary molecular basis for the prevention and treatment of YFD of Qinghuajiao trees.

General view of the transcriptome results. To elucidate the molecular response for this YFD of Qinghuajiao, we collected the leaves and inflorescences of symptomatic and asymptomatic Qinghuajiao trees of the same variety from the experimental field of Fruit Research Institute of Chongqing Academy of Agricultural Sciences, in Diaojia Town, Wutan Town and Mixin Town at March 21, 2019 (Fig. 2a), and noted as DAF, Diaojia asymptomatic flower; DSF, symptomatic flower; DAL, Diaojia asymptomatic leaf; DSL, symptomatic leaf; and so on.
Through high-throughput sequencing RNA-Seq of the 12 samples of these three sample points, we obtained 80.45 Gb total clean data, and the raw and clean reads of all samples were close to each other (Table 1). By using Trinity software (https:// github. com/ trini tyrna seq/ trini tyrna seq, Version v2.8.5), we assembled all clean data from scratch, and obtained 126,550 unigenes and 229,643 transcripts. The smallest and largest length were both 201 bp and 14,553 bp for the assembled unigenes and transcripts, and the N50 length was 1095 bp for unigenes, 1524 for transcripts (Table S1). And then the clean reads of each sample were compared with the reference sequence obtained by Trinity assembly, we obtained the mapping results of each sample, and we found a total of 21,054 genes that were detectable in all 12 samples ( Fig. 2b and Table 1). Principal component analysis (PCA) showed that the 12 samples could be divided into four groups: DSF, DAF, WSF, WAF; DSL, WSL, WAL; MSF, MAF, MSL; and DAL, MAL, which was correlated to the sample location and tissue type, only DAL and MAL were separated from the groups that they were supposed to be in (Fig. 2c). Furthermore, the results of sample correlation matrix based on gene expression levels were consistent with the PCA results (Fig. 2d), which emphasized the reproducibility and reliability of our experimental samples.
The total DEGs of asymptomatic and symptomatic Qinghuajiao trees. According to the sample types, we divided the 12 samples into 4 groups: AF (DAF, MAF, WAF), SF (DSF, MSF, WSF), AL (DAL, MAL, WAL), and SL (DSL, MSL, WSL). In order to avoid the influence of environmental factors, the samples collected from three different towns were used as replicated samples. For instance, DAF, WAF, and MAF were used as replicates for AF samples. By using DESeq2 (http:// bioco nduct or. org/ packa ges/ stats/ bioc/ DESeq2/, Version 1.24.0) software, we obtained 56 and 164 differentially expressed genes (DEGs, four fold expression difference ). Interestingly, we construct the diagram of Venn for these 2 gene sets, and, among the 220 DEGs, we found that only 33 DEGs were detected in both groups (Fig. 3b), which might be related to the YFD of Qinghuajiao.
GO and KEGG analysis for DEGs. To better understand the functions of DEGs, we carried out Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis by using the free online tools of Majorbio Cloud Platform (http:// www. major bio. com). For AL_vs_SL group, 22 DEGs were related to molecular function, in which "catalytic activity" (12) and "binding" (9) were the mainly GO terms; 25 DEGs touched upon biological process, and "metabolic process" (9) and "cellular process" (8) occupied a high percentage; 27 DEGs referred to cellular component, and "membrane part" (8) and "cell part" (6) made a major contribution to this category (Fig. 4a). In AF_vs_SF group, GO terms "catalytic activity" (53) and "binding" (52) also mainly contributed to the molecular function category (119); "catalytic activity" (31), "binding" (25), and "biological regulation" (9) were the major GO terms that related to biological process (90); GO terms "membrane part" (26), "cell part" (24), and "organelle" (18) occupied higher percentage in the cellular component category (Fig. 4b).
In the KEGG annotation, 12 DEGs were divided into 10 categories for AL_vs_SL group, and 35 DEGs were referred to 24 categories for AF_vs_SF group (Table 2), among which, "Phenylpropanoid biosynthesis", "Cyanoamino acid metabolism", "Starch and sucrose metabolism", "Amino sugar and nucleotide sugar metabolism", "RNA transport", and "Basal transcription factors" were shared for both groups. Interestingly, the "Phenylpropanoid biosynthesis" pathway has been reported to be involved in response to plant pathogen. Moreover, for AF_vs_SF group, 2 DEGs were related to the "Plant-pathogen interaction" pathway ( Table 2), implying that this YFD diseases was possibly caused by plant pathogen.
The 2 groups shared 33 common DEGs. As shown in Fig. 3b, the AL_vs_SL and AF_vs_SF groups shared 33 common DEGs. More interestingly, most of the 33 common DEGs showed a consistent regulatory trend, that is to say, they were up-regulated in curled leaves and yellow flowers of the symptomatic trees (Fig. 5a), implying that these 33 DEGs might closely related to the YFD of Qinghuajiao. In addition, we further confirmed the expression level of these 33 DEGs in symptomatic trees by real time qRT PCR, the results were consistent with the transcriptome analysis (Fig. 6). Interestingly, the plant genes (NO.3, 6, 17, 20) were almost not expressed in the samples from asymptomatic trees, while were highly expressed in symptomatic trees, implying   18 . These reports are all from the plant pathogen view to elaborate disease of Chinese prickly ash. In this study, we detected the transcriptome alteration of Qinghuajiao leaves and inflorescences that encountered the YFD disease, and obtained a total of 220 DEGs, which were involved in the molecular response of Qinghuajiao to this YFD. To our knowledge, this is the first transcriptome report in response to Qinghuajiao YFD. The transcriptome profiling of cucumber leaves infected with powdery mildew (PM) revealed that the complex regulatory network for PM resistance involves plant hormone signal transduction, phenylpropanoid biosynthesis, plant-pathogen interaction and the MAPK signalling pathway 19 . Based on the KEGG pathway and GO enrichment of transcriptome data, Guo et al. showed a complex regulatory network for PM resistance in pumpkin that may also involve hormone signal transduction pathways 20 . Chen et al. reported that genes associated with plant hormone signal transduction, detoxification, phenylpropanoid biosynthesis, photosynthesis and chlorophyll metabolism were significantly affected by CMV infection in yellow passion fruit 21 . Transcriptome results showed that the moderate resistance of Pinus pinaster to Fusarium circinatum may be explained by the expression profiles pertaining to early recognition of the pathogen, the induction of pathogenesis-related proteins and the activation of complex phytohormone signaling pathways 22 . Some KEGG annotations of Qinghuajiao encountered with YFD were also belong to these previously reported to be response to plant-pathogen interaction by transcriptomic analysis. For instance, the "Plant-pathogen interaction" and "plant hormone signal transduction" pathways for AF_vs_SF group, and the "phenylpropanoid biosynthesis" for both groups. These results are important first step for insights on the molecular mechanism for the interaction between Qinghuajiao and YFD, and suggest that the phenylpropanoid biosynthesis and plant hormone signal transduction pathways were involved in the responses of Qinghuajiao to YFD.  27 . The blast results of the up-regulated 30 common genes showed that most of these genes has highly homolog with virus RNA, implying that the YFD of Qinghuajiao might be related to virus. Further study will be conducted at the isolation of the virus that resulted in YFD of Qinghuajiao, the molecular character and the transmission routes of these viruses, so as to get the methods to suppress the transmission of these viruses, and thus to inhibit the further spread of YFD. In addition, although all most virus-infected plant diseases cause economic losses, Beaver-Kanuya and Harper had reported that four viruses infecting the pollen of Prunus species might have implications for biosecurity 28 . While the symptomatic Qinghuajiao trees have intumescent stamens, which could be used to improve male sterility in other crops.
Taken together, through high-throughput RNA-Seq analysis of 12 samples from three points, we revealed the molecular response of YFD of Qinghuajiao. We found a total 220 DEGs for the symptomatic trees. The Gene Ontology (GO) analysis showed that the major related GO terms of the DEGs were "catalytic activity" and "binding" to molecular function category, "metabolic process" (9) and "cellular process" to biological process category, and "membrane part" and "cell part" to cellular component category, respectively. The Kyoto Encyclopedia of Genes and Genomes (KEGG) results showed that the DEGs were divided into 10 categories for AL_vs_SL group www.nature.com/scientificreports/ and were referred to 24 categories for AF_vs_SF group. Particularly, the "phenylpropanoid biosynthesis" pathway that related to plant-pathogen interaction were found in both groups, and the "Plant-pathogen interaction" founded in AF_vs_SF group, implying that this Qinghuajiao YFD might cause by plant pathogen. Interestingly, we found 33 shared DEGs for the AL_vs_SL and AF_vs_SF groups, and most of them were all up-regulated in the samples from diseased trees. Among which, 24 genes were homologs to virus RNA, implying that this YFD of Qinghuajiao was related to virus. Although the mechanisms underlying the response of Qinghuajiao to YFD still require further research, the knowledge obtained from this study will serve as a useful genetic resource to facilitate further investigations of YFD in Qinghuajiao, and provide many possible directions for prevention and treatment of YFD of Qinghuajiao. Additionally, our study also provides a quick and preliminary approach to identify the causes of serious diseases occurring in plants with complex genomes under natural conditions.

Materials and methods
Materials. The green Chinese prickly ash (Zanthoxylum schinifolium, Qinghuajiao) was grown in Chongqing at natural conditions. The in inflorescences and leaves of asymptomatic and symptomatic Qinghuajiao trees of same variety (Qinghuajiao) were collected from the experimental field of Fruit  High-through put RNA-Seq. The inflorescence and leaf samples were collected at March 21, 2019, and was stored at − 80 °C at Biotechnology Research Center, Southwest University. The samples were then sent to Shanghai Majorbio Bio-pharm Technology Co., Ltd. to conduct High-through put RNA-Seq as previously described and with some modification 29 . Briefly, The mRNA was enriched using the oligo(dT) magnetic beads. The sample libraries were qualified and quantified by Agilent 2100 Bioanaylzer and Nanodrop2000. The library products were sequenced via Illumina Novaseq 6000. The raw data were filtered with the FASTQ_Quality_Filter tool from the FASTX-toolkit (http:// hanno nlab. cshl. edu/ fastx_ toolk it/, Version 0.0.14). And then using Trinity software (https:// github. com/ trini tyrna seq/ trini tyrna seq, Version v2.8.5) to assemble all sample clean data from scratch 13  RNA extraction and quantitative RT-PCR. Total RNA was extracted using the plant total RNA extraction kit (Tiangen, China). First-strand cDNA was synthesized from 1 µg total RNA using a reverse transcription kit with genomic DNA remover (Takara, Japan). Gene-specific primers are designed and listed in Table S2.

Data availability
The raw sequence data supporting the results of this article are available on the free online platform of Majorbio Cloud Platform (http:// www. major bio. com). Further information and requests for data and material should be directed to and will be fulfilled by Ming Luo (luo0424@126.com). The raw sequence data supporting the results of this article are also submitted to the Sequence Read Archive (SRA), and can be found by SRR14127436, SRR14127437, SRR14127438, SRR14127439, SRR14127440, SRR14127441, SRR14127442, SRR14127443,