Integrating long noncoding RNAs and mRNAs expression profiles of response to Plasmodiophora brassicae infection in Pakchoi (Brassica campestris ssp. chinensis Makino).

The biotrophic protist Plasmodiophora brassicae causes serious damage to Brassicaceae crops grown worldwide. However, the molecular mechanism of the Brassica rapa response remains has not been determined. Long noncoding RNA and mRNA expression profiles in response to Plasmodiophora brassicae infection were investigated using RNA-seq on the Chinese cabbage inbred line C22 infected with P. brassicae. Approximately 5,193 mRNAs were significantly differentially expressed, among which 1,345 were upregulated and 3,848 were downregulated. The GO enrichment analysis shows that most of these mRNAs are related to the defense response. Meanwhile, 114 significantly differentially expressed lncRNAs were identified, including 31 upregulated and 83 downregulated. Furthermore, a total of 2,344 interaction relationships were detected between 1,725 mRNAs and 103 lncRNAs with a correlation coefficient greater than 0.8. We also found 15 P. brassicaerelated mRNAs and 16 lncRNA interactions within the correlation network. The functional annotation showed that 15 mRNAs belong to defense response proteins (66.67%), protein phosphorylation (13.33%), root hair cell differentiation (13.33%) and regulation of salicylic acid biosynthetic process (6.67%). KEGG annotation showed that the vast majority of these genes are involved in the biosynthesis of secondary metabolism pathways and plant-pathogen interactions. These results provide a new perspective on lncRNA-mRNA network function and help to elucidate the molecular mechanism of P. brassicae infection.

Introduction

Clubroot, a soil-borne disease, has caused considerable damage to Brassicaceae crops [1, 2]. This disease is caused by the protist Plasmodiophora brassicae (P. brassicae), which can survive for up to 20 years in soil [3]. The two stages of P. brassicae, root-hair infection and cortical infection, play an important role in the infection process and make it difficult to control [4]. The Pakchoi (Brassica campestris ssp. chinensis Makino), also called non-heading Chinese cabbage, is one of the most important Brassica vegetable crops in China and other eastern Asian countries. Most Pakchoi cultivars are highly susceptible to the P. brassicae.

To date, considerable progress has been made in cultivating clubroot resistant (CR) crops. Genetic analysis and QTL mapping have identified some CR genes or loci in Brassica crops: CRa [5], Crr1a and Crr1b [6], CRb [7], Crr2 [8], Crr3 [9, 10], Crr4 [11], CRc and CRk [12], Rcr1 [13, 14], PbBa3.1 and PbBa3.3 [15], QS_B1.1 [16], and Pb-Br8 [17]. Three loci for clubroot resistance, Rcr4, Rcr8, Rcr9, have been revealed by Genotyping-by-sequencing, but they cannot be distinguished from the abovementioned loci [18]. Among them, CRa, Crr1a and CRb have been cloned. CRa and Crr1a contain Toll-interleukin receptor (TIR)—nucleotide-binding (NB)—leucine -rich repeats (LRRs) and CRb contains NB-LRRs, which are known to be responsible for race-specific resistance in higher plants [19, 20]. However, these genes or loci have been demonstrated to be responsible for race-dependent resistance [21], and the molecular mechanism of the Brassica rapa responsehas not been determined.

To date, a number of transcriptome sequencing projects have been employed to explore the molecular basis of the interaction between Brassica crops and P. brassicae. Twenty protein spots that were observed with changes in expression played a role in lignin synthesis, cytokinin synthesis, calcium steady-state, glycolysis, and oxygen activity in Brassica napus [22]. Then, the signaling and metabolic activity of jasmonate acid (JA) and ethylene (ET) were found to be upregulated significantly in resistant populations while genes involved in salicylic acid metabolic (SA) and signaling pathways were generally not elevated at 15 days post inoculation (dpi) [13]. Moreover, genes associated with pathogen-associated molecular patterns (PAMPs) and effector recognition, calcium ion influx, hormone signaling, pathogenesis-related (PR) genes, transcription factors, and cell wall modification showed different expression patterns between CR and clubroot-susceptible (CS) lines in Brassica rapa [23]. PR genes are involved in SA signaling which is important to clubroot resistance at the early stage after inoculation. In addition, it was proven that response changes in transcript levels under P. brassicae infection were primarily activated at the primary stage between Broccoli (Brassica oleracea var. italica) and wild Cabbage (Brassica macrocarpa Guss.) [24]. By comparing the transcriptome landscape between CS and CR Chinese cabbage lines, Jia et al. (2017) confirmed that the differentially expressed genes related to disease-resistance in CR lines enriched in calcium ion influx, glucosinolate biosynthesis, cell wall thickening, SA homeostasis, chitin metabolism and PR pathway. The upregulated genes in CS lines were mostly related to cell cycle control, cell division and energy production and conversion [2]. In addition, the Indole acetic acid (IAA) and cytokinin-related genes were found to affect the root swelling in clubroot development [2, 25].

LncRNAs are a set of RNA transcripts (>200 nt in length) which have no protein-coding ability. During the past several decades, a small number of long noncoding RNAs (lncRNAs) have been identified and shown to mediate various biological processes in plants [26], such as biotic and abiotic stress responses [27, 28]. In plant-pathogen interactions, some lncRNAs have been identified and shown to respond to (1) stripe rust pathogen stress in wheat [24]; (2) Fusarium oxysporum infection [29] and Pseudomonas syringe pv tomato DC3000 (ELF18-induced lncRNA) [30] in Arabidopsis thaliana; (3) Pectobacterium carotovorum in potato [31]; (4) Phytophthora infestans (lncRNA 16397) [32], tomato yellow leaf curl virus [33] and Phytophthora infestans in tomato (lncRNA23468) [34]; and (5) Sclerotinia sclerotiorum in Brassica napus [35]. The B. rapa and B. napus genome has a large number of lncRNAs [36, 37]. In addition, lncRNAs are demonstrated with the ability to be expressed broadly across many developmental times and in different tissue types [37]. However, only a few lncRNAs coexpressed with genes of temperature expression patterns were reported in Brassica rapa [36].

In this study, we first conducted a comprehensive analysis of intergrating long noncoding RNAs and mRNA expression profiles of response to Plasmodiophora brassicae infection in Brassica rapa L. and identified a great number of significant differentially expressed genes and some lncRNAs. The regulatory network of mRNA and lncRNA helps to elucide the Brassica rapa responses during P. brassicae infection and breeding of resistant CR cultivars.

Materials and methods

Ethics statement

This study was carried out in a phytotron. No specific permissions were required. The study did not involve any endangered or protected species.

Sample collection

Pakchoi inbred line CS22 is a cold tolerant type and susceptible to the 7th physiology race of Plasmodiophora Brassicae by using the inoculation method of Williams [38]. The pathogen was propagated on CS22 named CS22A, and the clubs in infected roots were stored at -20°C until required. All plants were sown in a growth chamber at 25/20°C (day/night) with a photoperiod of 14h containing. The CS22A plants were inoculated in a pot containing 5×106 spores per gram of dry soil. The root tissue samples were obtained by 6 weeks post inoculation. No infected root samples of CS22 were the control. For each treatment, the samples were immediately frozen in liquid nitrogen and then stored at −80°C until use. All plant materials examined in this study were obtained from Shanghai Academy of Agricultural Sciences.

RNA extraction, library construction, and sequencing

Total RNA was extracted from each root tissue sample using the mirVana miRNA Isolation Kit (Ambion) following the manufacturer’s protocol. RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The samples with RNA Integrity Number (RIN) ≥ 7 were subjected to the subsequent analysis. The libraries were constructed using TruSeq Stranded Total RNA with Ribo-Zero Gold according to the manufacturer’s instructions. The main steps of library construction and sequencing are as follows: (1) removing rRNA from total RNA, (2) breaking RNA into fragments, (3) RNA fragments are reverse-transcribed into cDNA, (4) adapter sequences are added to cDNA, and suitable fragment sizes are selected for the next step, and (5) PCR amplification. Then these libraries were sequenced in the Illumina HiSeqTM 2500 sequencing platform and 150 bp paired-end reads were generated.

Data filtering and transcriptome assembly

The RNA-seq data sets were analyzed as previously described [39]. High quality clean data were kept for downstream analysis after we use Trimmomatic v0.32 with ‘LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50’ to remove low-quality reads from the raw data, such as the reads containing adapters, the reads containing over 10% of poly (N), and low-quality reads (> 50% of the bases having Phred quality scores <10). Basic information of clean data was calculated, such as read number, base contents, Phred score (Q30) and GC content. Brassica rapa reference genome and gene model annotation files, which were downloaded from the Genome Database (http://brassicadb.org/). First, the index of the reference genome was built with Hisat-build, and then paired-end clean reads were aligned to the reference genome using Hisat with default parameters [40]. Second, the sam format result from hisat2 was translated to bam format by samtools [41], and then bam files from each library were assembled with Stringtie [42]. Stringtie was run with ‘-library-type fr-firststrand’, and other parameters were set as default. Last, each library results from Stringtie were merged to a final genome transcript feature file by cuffmerge [43].

Pipeline for LncRNA identify

To obtain putative lncRNAs, assembled novel transcripts were filtered following the steps according to the assembly results. (1) First, cuffcompare was used to compare the assembly transcript and reference transcript one by one [43], only transcripts annotated as “i”, “u”, “x”, and “o” representing a transfrag falling entirely within a reference intron, unknown intergenic transcript, exonic overlap with reference on the opposite strand, and generic exonic overlap with a reference transcript, respectively, were retained. (2) Second, the transcripts with a length of above 200 bp and with an exon number of more than 1 were kept for the next step. (3) Finally, four different methods were used to identify the coding potential of new transcripts, namely, Coding Potential Calculator (CPC) [44], Coding-Non-Coding Index (CNCI) [45], PLEK [46] and Pfam [47]. The methods were used to assess the coding potential of the remaining transcripts from step 2. Transcripts that were likely to contain a known protein-coding domain removed. Only transcripts considered to be lncRNAs via four methods will be kept for downstream analysis.

Identification of differentially expressed mRNA and lncRNA

Express [48] and bowtie2 [49] were used to calculate FPKM scores for the lncRNAs and coding genes in each library. Differentially expressed lncRNAs and mRNAs between any two libraries were identified by DESeq (release 3.2) [50]. P value < 0.05 and an absolute value of the fold change ≥ 2 were used as a threshold to evaluate the statistical significance of lncRNA and mRNA expression differences.

Quantitative real-time PCR validation

To validate the credibility of the findings of RNA analysis, mRNAs and lncRNAs were randomly selected for real-time PCR. Total RNA was collected from the root tissue samples of the two groups using TruSeq Stranded Total RNA LT—(with Ribo-Zero Plant). The SuperScript III First-Strand Synthesis System was used to reverse the transcription to cDNA. Quantitative RT-PCR was conducted in a ViiA 7 Real-time PCR System (Applied Biosystems) using PowerUp™ SYBR Green Master Mix (Applied Biosystems, Carlsbad, CA, USA). Tubulin beta-6 (TUB6) was used as an internal control to normalize the data [51]. The primers used in qRT-PCR and cDNA synthesis were designed in the laboratory and synthesized by OEBiotech (Shanghai OEBiotech. Co., Ltd, Shanghai, China) based on the sequences. Primers are listed in S1 Table. The reaction conditions were as follows: incubation at 95°C for 10 min, followed by 40 cycles of 95°C for 10 s and 60°C for 1 min. The relative expression levels were calculated using the 2−ΔΔCt method and were normalized to TUB6, as an endogenous reference transcript.

Functional enrichment of differentially expressed mRNA

The Gene Ontology (GO) database (http://www.geneontology.org) is a description database that was usually applied to elucidate the genetic regulatory network of interest by forming hierarchical categories according to the molecular function, biological process, and cellular component. The Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) is the main public database about pathways. GO annotation and pathway analysis were used to study the effects of all significant differentially expressed mRNAs. The p value is calculated using hypergeometric test. R program packages were used to elucidate the GO and KEGG for targets of significant differential enrichment. GO and KEGG terms with P values < 0.05 were recognized as significant enrichment.

Target gene prediction

The function of lncRNAs is mainly realized by cis acting on target genes. The basic principle of cis-acting target gene prediction holds that the function of lncRNA is related to the protein-coding genes adjacent to its coordinates; therefore, the mRNA adjacent to lncRNA is selected as its target gene. Target gene analysis method: Pearson correlation coefficients of lncRNA and mRNA ≥ 0.8 were required. LncRNA is determined as regulator if it is within 100 k upstream or within 100 k downstream of mRNAs.

LncRNA-mRNA co-expression network construction

According to the differentially expressed lncRNA and mRNA results, we constructed a regulatory network to identify the relationships between lncRNA genes and mRNA genes. The Pearson correlation test was used to calculate the correlation between differential lncRNA and mRNA expression data. Pearson’s correlation coefficients equal to or greater than 0.8 and a P value less than 0.05 were considered to be lncRNA-mRNA pairs. Arranging from small to large according to the p-value, we chose 600 top results to construct the regulatory network, and the lncRNA-mRNA pairs associated with disease resistance were also used for the network construction. Cytoscape software (Cytoscape Consortium, San Diego, CA, USA) was used to present the lncRNA-mRNA regulatory network relationship.

Accession number

The RNA-seq datasets used in this study can be found in the NCBI Gene Expression Omnibus under accession number: PRJNA528807.

Results

Overview of RNA sequencing

To elucidate lncRNA and mRNA expression patterns in response to Plasmodiophora brassicae infection in Brassica rapa L., 6 libraries were constructed from control and clubroot tissues (CS22A) (Fig 1) for three biological replicates and sequenced using the Illumina platform [52]. The raw data obtained from RNA-seq are available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA). A total of 296.14 million and 295.29 million raw data reads were obtained in the control and clubroot libraries, respectively. The number of reads after quality filtering of the above two libraries were 289.57 million and 289.73 million respectively. Approximately 61.45% of the reads were mapped to the reference genome (BRAD database http://brassicadb.org/) [53]. Q30 (reads with an average quality score > 30) reads were more than 93% and the GC content of all sequencing libraries were less than 52%. In this study, we identified 38,483 mRNAs and 1,492 lncRNAs in the control and clubroot libraries. Gene expression levels were calculated using the FPKM (fragments per kilobase of exon model per million mapped) method [54]. Among the identified lncRNAs, 659 lncRNAs were known based on the database of CANTATAdb 2.0 (http://cantata.amu.edu.pl/) [55]. The mRNA expression level varied from 0 to 9,527.5 among all libraries, with an average value of 20.2. The lncRNA expression level varied from 0 to 65,092.9 among the six libraries, with an average value of 218.5. Principle component analysis (PCA) of the transcriptome and lncRNA FPKM values for all samples showed more separation between control and clubroot samples, respectively. This result showed that the sequencing data could be used for further analysis.

Fig 1

Clubroot symptoms of CS22.

Plants were inoculated with the 7th physiology race of P. brassicae (CS22A), while the control was not subjected to pathogen inoculation.

mRNA and lncRNA expression profiles in Pakchoi

Compared to the control samples, 5,193 mRNAs were observed to be significantly differentially expressed (fold change ≥ 2 and P ≤0.05), including 1,345 upregulated and 3,848 downregulated in CS22A. In total, a number of 114 significantiy differentially expressed lncRNAs were identified, including 31 upregulated and 83 downregulated. The number of downregulated mRNAs and lncRNAs was higher than the number of upregulated. Clustering analysis of the top 40 most significantly differentially expressed mRNAs between control and CS22A is shown with a heatmap (Fig 2A, S2 Table), and the heatmap of the top 40 significantly differentially expressed lncRNAs is shown in Fig 2B

Fig 2

Heatmap of 40 significantly differentially expressed mRNAs (A) and lncRNAs (B).

The length distribution and categorization of identified lncRNAs were also analyzed (Fig 3). The length of lncRNAs ranged from 200 to 4,483 bp, with an average length of 658 bp. The most abundant lncRNAs were between 200–400 bp. The number of lncRNAs decreased as the length increased. The lncRNA lengths were mostly less than 2,000 bp. lncRNAs were categorized into four groups, intronic, intergenic, sense and antisense based on their location on the genome [56, 57]. The majority of lncRNAs (55.16%) were intergenic and located in intergenic regions. The rates of lncRNAs were 2.41%, 27.82% and 14.61% for intronic, sense and antisense, respectively. Because lncRNAs encode small RNAs, the sequences of the lncRNAs were mapped to small RNA precursors. Twenty-five small RNA families were mapped to fifteen lncRNAs.

Fig 3

Length distribution and categorization of identified lncRNAs.

(A) The length distribution of identified lncRNAs. X-axis: the length of LncRNAs. (B) The rate of lncRNAs based on their location on the genome.

To confirm the expression level of differentially expressed RNAs identified from the RNA sequencing data, qRT-PCR analysis was used to assay the expression level of 10 randomly selected differentially expressed RNAs and lncRNAs. The trend of expression changes of these select genes based on the qRT-PCR was similar to the sequencing data, which suggested that the RNA-seq data were reliable (Fig 4).

Fig 4

Randomly selected differentially expressed RNAs were analyzed using qRT-PCR.

The expression level was normalized using Y-axis: the relative expression of selected genes compared with control as indicated. Data are shown as the mean ± standard deviation of three independent experiments.

Functional annotation analysis of significant differentially expressed mRNAs

GO enrichment analysis was conducted on the significantly differentially expressed mRNAs (fold change ≥ 2 and P ≤ 0.05) to gain more insights into the function of these mRNAs which can be divided into three main functional groups (Fig 5). In biological processes, the top 40 GO terms of the upregulated mRNAs showed that the majority of the functions related to the defense response to bacterium (GO:0042742), the defense response to fungus (GO:0050832), the response to wounding (GO:0009611), the response to jasmonic acid (GO:0009753) and the response to toxic substances (GO:0009636) (Fig 5, S3 Table). GO categories of the downregulated genes were shown to be closely related to defense response (GO:0006952), defense response to bacterium (GO:0042742), auxin-activated signaling pathway (GO:0009734), response to wounding (GO:0009611), response to auxin (GO:0009733) and response to jasmonic acid (GO:0009753) (Fig 5, S3 Table). It can be assumed that the genes or proteins that the mRNAs code for are involved in the reaction. In the cellular component, upregulated genes were mapped to membrane (GO:0016020), thylakoid (GO:0009579 and GO:0044436) and membrane protein complex (GO:0098796), while downregulated genes were mapped to intrinsic component of membrane (GO:0031224), integral component of membrane (GO:0016021), and membrane (GO:0044425 and GO:0016020). Regarding the molecular function, the enriched GO terms targeted by upregulated genes included catalytic activity (GO:0003824), oxidoreductase activity (GO:0016491), cofactor binding (GO:0048037) and transporter activity (GO:0005215), the enriched GO terms targeted by downregulated genes included catalytic activity (GO:0003824), transferase activity (GO:0016740) and oxidoreductase activity (GO:0016491).

Fig 5

Heatmaps of significantly differentially expressed mRNAs classified by biological process.

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the significantly differentially expressed mRNAs indicated that the top 3 KEGG terms for downregulated mRNAs were associated with plant hormone signal transduction (ko04075), MAPK (mitogen-activated protein kinase) signaling pathway (ko04016) and ABC (ATP-binding cassette transporters) transporters (ko02010), while the top 3 KEGG terms for upregulated mRNAs were associated with biofilm formation (ko02026), drug metabolism (ko00982) and metabolism of xenobiotics by cytochrome (ko00980). The plant hormone signal transduction pathway included 8 plant hormones that contained such acides as jasmonic acid and salicylic acid (https://www.genome.jp/dbget-bin/www_bget?ko04075). Therefore, these genes likely play an important role in the interaction in the infected process.

Target analysis for cis-regulated lncRNAs and their function annotation in significantly differentially expressed lncRNAs

Previous studies reported that lncRNAs regulated neighboring or overlapping genes and might show linked function or co-expression with their target genes [58-60]. Significant differentially expressed mRNAs located within 100 kb windows upstream or downstream of the lncRNAs were used to calculate the Pearson correlation coefficient for further analysis. A total of 2,344 interaction relationships (1,479 positive and 865 negative correlation) were detected between 1,725 mRNA and 103 lncRNA with a correlation coefficient greater than 0.8. The GO analysis was based on biological processes for all potential target mRNAs. Functional analysis showed that the upregulated co-expressed mRNAs of the neighboring lncRNAs were enriched in 39 GO terms in biological processes, and many of the GO terms were closely related to the regulation of gene expression (Table 1). The downregulated mRNAs of the potential lncRNA targets showed that the term GO enrichment was closely related to the response to stimulus (GO:0050896), response to stress (GO:0006950), defense response (GO:0006952) and response to biotic stimulus (GO:0009607) (Table 2). Regarding the molecular function, the enriched GO terms targeted by upregulated genes included oxidoreductase activity (GO:0016491), cofactor binding (GO:0048037), and transmembrane transporter activity (GO:0022857), and the enriched GO terms targeted by downregulated genes included catalytic activity (GO:0003824), transferase activity (GO:0016740), and oxidoreductase activity (GO:0016491). In the cellular component, the upregulated gene showed that the majority of the function related to membrane protein complex (GO:0098796) and thylakoid (GO:0044436 and GO:0009579). GO categories of the downregulated genes were shown to be closely related to integral component of membrane (GO:0016021), intrinsic component of membrane (GO:0031224), and membrane (GO:0044425 and GO:0016020). The expression levels of the downregulated mRNAs for the above four GO terms (40 mRNAs) and 16 lncRNAs that regulated these RNAs are shown in Fig 6 and S4 Table, and all these mRNAs and lncRNAs were significantly differentially expressed between the control and clubroot groups. These results suggest that the principal functions of these lncRNAs may be the regulation of gene expression and play an important role in the clubroot infection process.

Table 1

Gene Ontology (GO) enrichment for the significantly upregulated co-expressed mRNAs of the neighboring lncRNAs.

GO ID	Term	Annotated	Significant	Expected	Classic Fisher
GO:0015979	photosynthesis	122	9	1.65	3.90E-05
GO:0006091	generation of precursor metabolites and energy	75	7	1.01	6.50E-05
GO:0009765	photosynthesis, light harvesting	34	5	0.46	8.50E-05
GO:0019684	photosynthesis, light reaction	51	5	0.69	0.0006
GO:0009414	response to water deprivation	9	2	0.12	0.0061
GO:0055114	oxidation-reduction process	1431	30	19.31	0.0088
GO:0019438	aromatic compound biosynthetic process	1660	33	22.4	0.0129
GO:0018130	heterocycle biosynthetic process	1667	33	22.49	0.0136
GO:0006811	ion transport	485	13	6.54	0.0141
GO:0005985	sucrose metabolic process	14	2	0.19	0.0148
GO:1901362	organic cyclic compound biosynthetic process	1677	33	22.63	0.0148
GO:0006355	regulation of transcription, DNA-template	1507	30	20.33	0.0175
GO:1903506	regulation of nucleic acid-templated transcription	1507	30	20.33	0.0175
GO:2001141	regulation of RNA biosynthetic process	1507	30	20.33	0.0175
GO:0051252	regulation of RNA metabolic process	1512	30	20.4	0.0183
GO:0010556	regulation of macromolecule biosynthetic	1513	30	20.42	0.0184
GO:2000112	regulation of cellular macromolecule biosynthetic	1513	30	20.42	0.0184
GO:0009889	regulation of biosynthetic process	1515	30	20.44	0.0187
GO:0031326	regulation of cellular biosynthetic process	1515	30	20.44	0.0187
GO:0019219	regulation of nucleobase-containing compound metabolic process	1517	30	20.47	0.019
GO:0034654	nucleobase-containing compound biosynthetic	1649	32	22.25	0.0198
GO:0051171	regulation of nitrogen compound metabolic	1538	30	20.75	0.0226
GO:0080090	regulation of primary metabolic process	1540	30	20.78	0.023
GO:0009415	response to water	18	2	0.24	0.024
GO:0031323	regulation of cellular metabolic process	1549	30	20.9	0.0247
GO:0006487	protein N-linked glycosylation	2	1	0.03	0.0268
GO:0009072	aromatic amino acid family metabolic process	2	1	0.03	0.0268
GO:0006351	transcription, DNA-templated	1630	31	21.99	0.0281
GO:0032774	RNA biosynthetic process	1630	31	21.99	0.0281
GO:0097659	nucleic acid-templated transcription	1630	31	21.99	0.0281
GO:0010468	regulation of gene expression	1599	30	21.58	0.0361
GO:0006772	thiamine metabolic process	3	1	0.04	0.0399
GO:0009228	thiamine biosynthetic process	3	1	0.04	0.0399
GO:0042723	thiamine-containing compound metabolic process	3	1	0.04	0.0399
GO:0042724	thiamine-containing compound biosynthetic	3	1	0.04	0.0399
GO:0044070	regulation of anion transport	3	1	0.04	0.0399
GO:0060255	regulation of macromolecule metabolic process	1624	30	21.91	0.0431
GO:0019222	regulation of metabolic process	1635	30	22.06	0.0466
GO:0006812	cation transport	350	9	4.72	0.0475

Annotated: number of genes that are annotated with the GO-term.

Significant: number of genes belonging to the term that are annotated with the GO-term.

Expected: an estimate of the number of genes a node of size annotated would have if the significant genes were to be randomly selected from the gene universe.

Classic Fisher: p-values computed by Fisher’s exact test

Table 2

Gene Ontology enrichment for the significantly downregulated co-expressed mRNAs of the neighboring lncRNAs.

GO ID	Term	Annotated	Significant	Expected	classicFisher
GO:0006979	response to oxidative stress	136	20	5.22	2.20E-07
GO:0006950	response to stress	552	42	21.18	1.60E-05
GO:0032989	cellular component morphogenesis	10	4	0.38	0.00037
GO:0055114	oxidation-reduction process	1431	78	54.9	0.00073
GO:0050896	response to stimulus	1190	67	45.65	0.00078
GO:0048869	cellular developmental process	14	4	0.54	0.00157
GO:0001558	regulation of cell growth	8	3	0.31	0.00272
GO:0009826	unidimensional cell growth	8	3	0.31	0.00272
GO:0042814	monopolar cell growth	8	3	0.31	0.00272
GO:0051510	regulation of unidimensional cell growth	8	3	0.31	0.00272
GO:0051513	regulation of monopolar cell growth	8	3	0.31	0.00272
GO:0060560	developmental growth involved in morphogenesis	8	3	0.31	0.00272
GO:0000902	cell morphogenesis	9	3	0.35	0.00396
GO:0022603	regulation of anatomical structure morphogenesis	9	3	0.35	0.00396
GO:0022604	regulation of cell morphogenesis	9	3	0.35	0.00396
GO:0009653	anatomical structure morphogenesis	19	4	0.73	0.00523
GO:0031667	response to nutrient levels	10	3	0.38	0.0055
GO:0031669	cellular response to nutrient levels	10	3	0.38	0.0055
GO:0040008	regulation of growth	10	3	0.38	0.0055
GO:0048589	developmental growth	10	3	0.38	0.0055
GO:0048638	regulation of developmental growth	10	3	0.38	0.0055
GO:0006952	defense response	155	13	5.95	0.00664
GO:0008272	sulfate transport	22	4	0.84	0.00902
GO:0072348	sulfur compound transport	22	4	0.84	0.00902
GO:0015698	inorganic anion transport	63	7	2.42	0.01011
GO:0009991	response to extracellular stimulus	13	3	0.5	0.01203
GO:0031668	cellular response to extracellular stimulus	13	3	0.5	0.01203
GO:0071496	cellular response to external stimulus	13	3	0.5	0.01203
GO:0006820	anion transport	115	10	4.41	0.01285
GO:0009267	cellular response to starvation	5	2	0.19	0.01359
GO:0016036	cellular response to phosphate starvation	5	2	0.19	0.01359
GO:0042594	response to starvation	5	2	0.19	0.01359
GO:0050793	regulation of developmental process	25	4	0.96	0.01425
GO:0051128	regulation of cellular component organization	26	4	1	0.01635
GO:0005984	disaccharide metabolic process	46	5	1.76	0.03061
GO:0009311	oligosaccharide metabolic process	47	5	1.8	0.03321
GO:0009607	response to biotic stimulus	65	6	2.49	0.0379
GO:0010208	pollen wall assembly	1	1	0.04	0.03836
GO:0010584	pollen exine formation	1	1	0.04	0.03836
GO:0010927	cellular component assembly involved in morphogenesis	1	1	0.04	0.03836
GO:0080110	sporopollenin biosynthetic process	1	1	0.04	0.03836
GO:0085029	extracellular matrix assembly	1	1	0.04	0.03836

Fig 6

Expression levels of the 40 downregulated mRNAs (A) and 16 lncRNAs (B).

LncRNA-mRNA co-expression analysis of response to clubroot infection process

In this study, all significant differentially expressed lncRNAs and mRNAs were used to calculate the Pearsoncorrelation coefficients based on their expression level. The top 600 potential lncRNA-mRNA regulated pairs whose Pearson correlation coefficient greater than 0.8 were used to construct the regulatory network (S1 Fig). The 40 clubroot diseases related to mRNAs and 16 lncRNAs targeting these significantly differentially expressed mRNAs were also used to construct the correlation network of lncRNA-mRNA. In total, the resulting lncRNA:mRNA association network had 31 nodes and 19 connections between the 15 mRNAs and 16 lncRNAs (Fig 7, S5 Table). Among these molecules, most of mRNAs and lncRNAs are significantly downregulated.This regulation network indicated that four lncRNAs were predicted to be targets of 2 lncRNAs. BraA07g029760.3C and BraA07g0285503C were both targeted by lncRNA TCONS-00034121. In addition, the other three genes were all targeted by lncRNA TCONS-00049044. These results suggest that the expression profiles of mRNA and lncRNA are significantly correlated.

Fig 7

LncRNA-mRNA correlation network of respose to clubroot infection process.

To elucidate the lncRNA-mRNA co-expression network, we annotated the function of the target genes by comparison with Arabidopsis. The annotation showed that they belonged to defense response proteins (66.67%), protein phosphorylation (13.33%), root hair cell differentiation (13.33%) and regulation of the salicylic acid biosynthetic process (6.67%) (Table 3). KEGG annotation showed that the vast majority of the genes involved in the biosynthesis of secondary metabolism pathways and plant-pathogen interactions.

Table 3

Annotation of 15 mRNAs involved in LncRNA-mRNA co-expression network by comparison with the Arabidopsis genome.

Gene	Gene description	Arabidopsis	Functional annotation
BraA01g015860.3C	U-box domain-containing protein 35	AT4G25160	protein phosphorylation, protein ubiquitination
BraA02g012160.3C	calmodulin-binding protein 60 B-like	AT5G57580	regulation of salicylic acid biosynthetic process
BraA02g015930.3C	U-box domain-containing protein 35	AT5G51270	protein phosphorylation, protein ubiquitination,
BraA02g044230.3C	defensin-like protein 6	AT5G63660	defense response,defense response to fungus,killing of cells of other organisms,
BraA05g006080.3C	nematode resistance protein-like HSPRO2	AT2G40000	defense response to bacterium, incompatible interaction, response to oxidative stress, response to salicylic acid, tryptophan catabolic process to kynurenine
BraA06g039160.3C	universal stress protein PHOS32	AT2G03720	root hair cell differentiation,
BraA07g028550.3C	protein SAR DEFICIENT 1-like	AT1G73805	cellular response to molecule of bacterial origin, defense response to bacterium,defense response to oomycetes, plant-type hypersensitive response, positive regulation of defense response to bacterium, regulation of salicylic acid biosynthetic process, regulation of systemic acquired resistance,regulation of transcription, DNA-templated,response to UV-B, response to bacterium
BraA07g029760.3C	MLP-like protein 31	AT1G70850	defense response
BraA08g026660.3C	MLP-like protein 31	AT1G70830	defense response
BraA09g016310.3C	MLO-like protein 6	AT1G61560	defense response, defense response to fungus,incompatible interaction, response to biotic stimulus
BraA09g023180.3C	MLP-like protein 328	AT2G01520	defense response, response to phenylpropanoid, response to zinc ion, vegetative to reproductive phase transition of meristem
BraA09g023480.3C	defensin-like protein 1	AT2G02130	defense response, defense response to fungus,killing of cells of other organism
BraA09g024130.3C	universal stress protein PHOS32-like	AT2G03720	root hair cell differentiation
BraA09g040010.3C	MLP-like protein 43	AT1G35310	defense response
BraA10g020100.3C	Polyketide cyclase/dehydrase andlipid transport superfamily protein	AT1G70860	defense response

Discussion

In the present study, RNA–seq technology was used to investigate the global lncRNA-mRNA regulatory network between the B. rapa line before and after P. brassicae infection. The results of the differentially expressed analysis showed that the number of significantly differentially expressed mRNAs and lncRNAs were approximately 3 times higher in downregulated than in upregulated number. A total of 5193 and 114 mRNAs and lncRNAs were significantly differentially expressed. These results showed that a more complicated regulatory network exists in the clubroot infected plant. The GO annotation of the potential lncRNA targets showed that most upregulated significantly differentially expressed mRNAs were involved in the regulation of gene expression, and the downregulated significantly expressed mRNAs were closely related to stimulus, response to stress, defense response and response to biotic stimulus. The results of KEGG pathway analysis for the above mentioned lncRNA targets showed that they involved in the plant hormone signal transduction. The same conclusion was also reported in previous research [2, 23, 25, 61]. Jasmonic acid and salicylic acid regulate disease resistance in Arabidopsis [62]. The MAPK signaling pathway plays a pivotal role in the cellular processes such as proliferation, apoptosis, and gene regulation [63]. The metabolism of drug and xenobiotic pathway function is to oxidize small foreign organic molecules, such as toxins or drugs [64]. These results suggest that clubroot resistance and some cellular biological processes may be repressed during pathogen infection. The reaction mechanism that responds to xenobiotics may be activated during pathogen infection. Our results provide a distinct landscape in regard to the molecular mechanisms underlying P. brassicae infection.

Some quantitative trait loci (QTLs) related to clubroot diseases have been identified [8–10, 13, 18, 65, 66]. LncRNAs are a group of endogenous RNAs that function as regulators of gene expression, and may play an important role in several biological processes of plants [24]. LncRNA COLDAIR was reported to be required for establishing stable repressive chromatin at FLOWERING LOCUS [67]. LncRNA ASL can be regulated by ATRRP6L to modulate H3K27me3 levels functions in the autonomous pathway in Arabidopsis [68]. Therefore, we first investigated the lncRNA response to Plasmodiophora brassicae infection in Pakchoi and attempted to identify genes regulated by lncRNAs. The markers of QTL intervals that have been identified were mapped to the genome to examine the position relation of the QTLs and the genes (a total of 15 mRNAs and 16 lncRNAs) that were identified as related to clubroot disease in this study. The results show that lncRNA TCONS_00007793 localizes near the QTL Anju1 region on Chromosome A02 [11], two lncRNAs (TCONS_00007004, TCONS_00007046) localize near the QTL Rcr8, which was identified on Chromosome A02 [18], lncRNA TCONS_00014032 localizes near the QTL CRd, which was identified on Chromosome A3 [65], lncRNA TCONS_00038153 localizes near the QTL CRs, which were identified on Chromosome A8 [66], lncRNAs (TCONS_00034121 and TCONS_00036594) localizes near the QTL qBrCR38-1, identified by the bulked segregant analysis (BSA) method [69] on Chromosome A07, lncRNA TCONS_00041523 localized near the QTL qBrCR38-2, which has been identified on Chromosome A08 in the same experience. These lncRNAs associated with the QTL regions maybe have the function of regulating gene expression [70].

We investigated the expression patterns of lncRNAs and mRNAs and constructed a lncRNA-mRNA regulatory network for P. brassicae infected Pakchoi and control. This network can provide a global view of all possible lncRNA-coding gene expression associations based on high-through RNA-seq data. The functional annotation shows that these lncRNAs might exhibit coordinating roles towards transcriptional regulation of the defense responsive genes. KEGG annotation shows that these genes, targeted by lncRNAs, are involved in the biosynthesis of secondary metabolism pathways which are essential for many physiological processes in plants, including pathogen invasion [71]. Although many lncRNAs have been found, their biological functions remain unclear. Further research on the specific role(s) of these lncRNAs will provide additional information regarding their detailed roles in pathogen defense.

Primers used in this study.

(DOC)

Click here for additional data file.

Information of 40 significantly differentially expressed mRNAs in Fig 2A.

(XLSX)

Click here for additional data file.

Information of mRNAs in Fig 5.

(XLSX)

Click here for additional data file.

Information of the 40 mRNAs in Fig 6.

(XLSX)

Click here for additional data file.

Information of the co-expression network in Fig 7.

(XLSX)

Click here for additional data file.

LncRNA-mRNA correlation network of the top 600 potential lncRNA-mRNA regulated pairs based on all significant differentially expressed LncRNAs and mRNAs.

The arrow and circle nodes denote lncRNA and mRNA, respectively. Each gray edge denotes a potential target relationship between a gene or lncRNA.

(PDF)

Click here for additional data file.

28 Aug 2019

PONE-D-19-15638

Intergrating long noncoding RNAs and mRNAs expression profiles of response to Plasmodiophora brassicae infection in Pakchoi (Brassica campestris ssp. chinensis Makino)

PLOS ONE

Dear Dr. Yuying Zhu,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

The topic of the manuscript is very interesting, however both the reviewers raised many criticisms that have to be addressed and substantial changes are needed to make the manuscript suitable for publication.

We would appreciate receiving your revised manuscript by Oct 12 2019 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Serena Aceto, Ph.D.

Academic Editor

PLOS ONE

Journal Requirements:

1. When submitting your revision, we need you to address these additional requirements.

Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. We note that you are reporting an analysis of a microarray, next-generation sequencing, or deep sequencing data set. PLOS requires that authors comply with field-specific standards for preparation, recording, and deposition of data in repositories appropriate to their field. Please upload these data to a stable, public repository (such as ArrayExpress, Gene Expression Omnibus (GEO), DNA Data Bank of Japan (DDBJ), NCBI GenBank, NCBI Sequence Read Archive, or EMBL Nucleotide Sequence Database (ENA)). In your revised cover letter, please provide the relevant accession numbers that may be used to access these data. For a full list of recommended repositories, see http://journals.plos.org/plosone/s/data-availability#loc-omics or http://journals.plos.org/plosone/s/data-availability#loc-sequencing.

3. Please amend the manuscript submission data (via Edit Submission) to include author Dandan Xi.

4. Please amend your authorship list in your manuscript file to include author Wen Zhai.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: No

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: I Don't Know

Reviewer #2: N/A

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: No

Reviewer #2: No

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: This manuscript describes the changes in the expression of protein coding genes (mRNA) as well as those encoding long non-coding (lnc)RNAs in Brassica rapa in response to infection with the clubroot causing pathogen, Plasmodiophora brassicae. It attempts to integrate the findings with respect to lncRNA expression and mRNA expression.

Experiments have been designed and conducted appropriately although there are some significant issues. Please see major comments below.

The manuscript has been generally written well but there are still many issues--see minor comment for a non-exhaustive list of errors. It is impossible for me to list all the errors and the authors should ensure that the submission is free from these types of mistakes.

Major comments:

1. The samples were collected 6 weeks after sowing when full-blown root galls had developed. I do not understand the logic of investigating molecular responses after the development of galls. It would have been more revealing at earlier stages of pathogenesis. In my opinion, earlier stages in the disease should be sampled and mRNA and lncRNA profiles investigated.

2.Were there any pathogen mRNA and lncRNA detected? This is such a late stage that this may have been possible. IF so, it should be reported and if not, a statement should be made to that effect.

3. Line 376, are you talking about QTLs associated with resistance? However, in this case you used susceptible lines. This part of the discussion is confusing and should be clarified.

4. The Discussion section is extremely weak, or non-existent. It has be significantly strengthened to discuss their findings and integrating them with the information that is available in the literature. Currently this section is a sum total of 31 lines or just about one page. The most important part of these types of descriptive articles is their discussion about possible biological relevance of their findings. This should be rewritten in its entirety.

5. Why does the title of the manuscript say Brassica campestris? I thought that Pakchoi is Brassica rapa? Please be accurate and consistent.

Minor comments:

1. Some typographical and grammatical errors in the abstract include, line 11 “researches”; line 12 “works”; line 13 “….profiles of response to…”.

2. Line 28, “…a kind of soil-borne disease…”. It is a soil-borne disease, NOT “a kind-of soil-borne disease”.

3. Line 36, not “clubroot resistance crops” but “clubroot resistant crops”.

4. Line 32, Brassica rapa????

5. Lines 53 and 60, specify the Brassica species! In general, please specify the Brassica (or any plant) species when you are talking about the results from another study.

6. Line 69—do not begin a sentence with “And”.

7. Line 79, wrong reference format.

8. Line 122, NOT “reversely” transcribed.

9. Line 181, should be “Functional” enrichment.

10. Line 281, should be “Functional” annotation.

11. Line 345, “Pearson” correlation coefficient NOT “person”.

Reviewer #2: I have carefully read the manuscript PONE-D-19-15638 “Integrating long noncoding RNAs and mRNAs expression profiles of response to Plasmodiophora brassicae infection in Pakchoi (Brassica campestris ssp. Chinensis Makino)” by Zhu et al. The topic of the manuscript is interesting and the integration of lncRNAs to increase our understanding of clubroot disease is timely and novel. However, there are some shortcomings of the manuscript in its present form most notably the absence of an acceptable discussion of the results. More detailed comments can be found below.

The English of the manuscript is understandable, but I would strongly advice the authors to consult someone proficient in scientific English for language editing as there are several grammatical errors in the MS. There are many sloppy errors (inconsistent spelling, spelling mistakes, etc) throughout the text. I have pointed out some issues in the detailed comments below, but this list is not exhaustive.

Overall the manuscript would add important findings to our understanding of clubroot disease, but I strongly recommend that the authors considerably revise the manuscript to highlight and describe the very interesting findings they collected.

Detailed comments:

Abstract:

L 11: delete “Although lots of researches have been conducted during past decades” – grammer errors and IMHO this is not important in the abstract.

L17: specify which type of enrichment analyses

L19: define which type of interaction relationship, list the most important groups of interactions to provide more info for the reader on the biological relevance of those rather than just providing numbers.

L21: change “15 clubroot disease related” into “15 P. brassicae mRNAs”, maybe provide more detail on these RNAs.

Introduction

The introduction informative about the aims and guides the reader to the topic of the manuscript. Some phrases are difficult to understand, especially if the reader is not familiar with the original literature cited, therefore I would advise the authors to carefully edit the introduction.

L29: no need to write P. brassicae in brackets.

L60: “between CR and clubroot-resistant (CS) lines” this sentence is not clear, should this be clubroot susceptible (CS)?

L60: The sentence starting with “It was also found….” is unclear, please rephrase and explain what you mean with “updated SA function”.

Material and Methods

L102ff: Sample collection: the authors state that they use race 7 of P. brassice (also check for consistent spelling!). Was this a single spore isolate or a field population showing the characteristics of race 7 using the Williams differential? For all further experiments it will make a huge difference if the experiments were conducted using a single spore isolate or a field population.

L105ff: when were the plants inoculated and with approximately which amount of P. brassicae spores?

L107: please change “6 weeks after sowing” by the days post inoculation as this will provide more information on the disease progression

L129ff: what does “clean data” and “dirty reads” refer to? Please rephrase this and avoid laboratory jargon. Having a flow chart of the procedure including versions of the software in the Supplement is usually very helpful for such methods sections (also for LncRNA).

L172: state the type of mastermix used

L181ff: This section is not fully clear, please provide more information on the statistical analyses (R packages, settings, assumptions etc used in the analyses)

L188: The sentence “the p value…” is incomplete, please rephrase

L193ff: this section is not clear.

L193: lncRNAs - please make sure the correct spelling throughout the MS

Results

L217: use control or C instead of CK, this abbreviation is not intuitive.

L223: use “the number of reads after qality filtering (or quality control) were….” Instead of “clean reads”

L225: the sentence startin with Q30 is incomplete

L293: “these findings indicate that mRNAs were participating in the defence of clubroot” – This sentence is not correct. There is not evidence that the mRNAs themselves are involved in clubroot defence, rather it can be assumed that the genes/proteins the mRNA codes for are involved in the reaction. Please rephrase.

L301: The phrase “previous researches” is odd, rephrase to “previous studies” or similar.

L307: why were the GO terms restricted to “biological process” and why were molecular function and cellular component omitted? Is there a reason for this choice?

L346: is it really “person correlation coefficient” or should it be Pearson’s

Tables and Figures:

Overall the figure and table capitations could be a bit more informative and descriptive. Not every image can be easily understood, therefore some comments on specific issues:

Tables: Please describe what the different values stand for. I suspect that “annotated” means the number of transcripts that were assigned this GO term in the full dataset – or does this refer to the number of transcripts that were upregulated in clubroot tissue and assigned a certain GO term? Significant – is the number of these GO terms that were significant in which respect – significantly up/downregulated? What does “expected” mean and where does this come from? Classic Fisher refers to what – the significance of the GO-term, the up/downregulated GO-term, the transcripts?

Fig 2, 5, 6 Heatmap figures: can you provide any biological information to the genes other than the BRA accession? Maybe adding some sort of functional annotation (GO term, gene name, function of the gene,….). Which values are displayed? FPKM values, and if yes were these normalised? DEGs – but then which values were compared to give the values?

Fig3A: what do the numbers on the x-axis refer to? Please describe in the figure capitation or in the image.

Fig7: This figure is not very informative in the present form and also its not a network but a series of correlated genes. Please add information and annotations (which genes are we looking at?) to the figure or convert it into a table (which would provide more information on the individual correlations)

Fig S1: please provide a description of the figure. It is nearly impossible to understand this figure the way it is currently presented.

Discussion:

Unfortunately, the discussion feels very incomplete especially as the authors present a number of fascinating results. The authors fail to discuss what the findings of the correlation of lncRNAs and mRNAs presented mean for the biology of clubroot disease. The concept of lncRNAs is employed to clubroot for the first time, so there are plenty of factors that can be discussed and described here. Also there is no comparison to other transcriptomics studies of which there are plenty on a multitude of clubroot hots, resistant and susceptible interactions, on the intraplant variation etc. Please use this pool of references to discuss the results in a broader context.

L357ff: the first paragraph of the discussion is mostly results. Please move the description of the KEGG analyses into the results section, where only GO terms are described currently. Many of these processes have been identified in previous transcriptomic studies of clubroot, please cite those studies and compare their results to the ones generated in this study.

L376ff: This information is interesting, but most of the data are not yet available. Therefore the validity of more than half of the second paragraph of the discussion cannot be assessed.

Data availability

Data are available.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

16 Oct 2019

Reviewer #1

1. This manuscript describes the changes in the expression of protein coding genes (mRNA) as well as those encoding long non-coding (lnc)RNAs in Brassica rapa in response to infection with the clubroot causing pathogen, Plasmodiophorabrassicae. It attempts to integrate the findings with respect to lncRNA expression and mRNA expression.Experiments have been designed and conducted appropriately although there are some significant issues. Please see major comments below.The manuscript has been generally written well but there are still many issues--see minor comment for a non-exhaustive list of errors. It is impossible for me to list all the errors and the authors should ensure that the submission is free from these types of mistakes.

Response: We appreciated reviewer’s interest in our work. We have corrected these issues in the revised manuscript according to the reviewers’ comments.

2. The samples were collected 6 weeks after sowing when full-blown root galls had developed. I do not understand the logic of investigating molecular responses after the development of galls. It would have been more revealing at earlier stages of pathogenesis. In my opinion, earlier stages in the disease should be sampled and mRNA and lncRNA profiles investigated.

Response: All plants were sown in a pot containing 5×106 spores per gram of dry soil.In 6 weeks after sowing, the plant has grown to 4or5 leaves. In this development stage, infected root showed enlargement. Therefore, our study sampled in the late stage of disease. Please refer to line111-115.

3.Were there any pathogen mRNA and lncRNA detected? This is such a late stage that this may have been possible. IF so, it should be reported and if not, a statement should be made to that effect.

Response: Yes, we actually find pathogen mRNA in our samples. About 3.68% reads of CS22A(infected root) blast to Plasmodiophorabrassicae genome.

4. Line 376, are you talking about QTLs associated with resistance? However, in this case you used susceptible lines. This part of the discussion is confusing and should be clarified.

Response: This is a very good point. Some quantitative trait loci (QTLs) which are related to clubroot diseases have been identified. The lncRNAs are a group of endogenous RNAs that function as regulators of gene expression and may play an important role in several biological processes of plants. So we want to find some negative or positive correlation of lncRNAs and mRNAs. Most of the 15 mRNAs are belonged to defense response proteins by compared with Arabidopsis. However, little is known about the 16 lncRNAs. We agree that further research on these lncRNAs will provide additional information about their detailed roles in pathogen defense.

5. The Discussion section is extremely weak, or non-existent. It has be significantly strengthened to discuss their findings and integrating them with the information that is available in the literature. Currently this section is a sum total of 31 lines or just about one page. The most important part of these types of descriptive articles is their discussion about possible biological relevance of their findings. This should be rewritten in its entirety.

Response: The Discussion sectionhas been thoroughly rewritten and the confusing sentence has been corrected. Please refer to line 418-476.

6. Why does the title of the manuscript say Brassica campestris? I thought that Pakchoi is Brassica rapa? Please be accurate and consistent.

Response: Pakchoi is also named non-heading Chinese cabbage. There have been a lot of studies in the past that have named Brassica campestris, such as Du et al.2008, Ma et al. 2010, Zhu et al. 2017, Fan et al 2019.

Du S and Y. Zhang, et al. (2008)."Regulation of nitrate reductase by nitric oxide in Chinese cabbage pakchoi (Brassica chinensis L.)." Plant Cell Environ 31 (2): 195-204.

Ma, J. and X. Hou, et al. "Cloning and Characterization of the BcTuR3 Gene Related to Resistance to Turnip Mosaic Virus (TuMV) from Non-heading Chinese Cabbage." Plant Molecular Biology Reporter. 2010, 28 (4): 588-596.

Fan, X., Xue, F., Song, B., et al. Effects of Blue and Red Light On Growth And Nitrate Metabolism In Pakchoi. Open Chemistry, 2019, 17(1):456-464.

Zhu H, Li X, Zhai W, et al. Effects of low light on photosynthetic properties, antioxidant enzyme activity, and anthocyanin accumulation in purple pak-choi(Brassica campestris ssp. Chinensis Makino).Plos One, 2017, 12(6):e0179305.

7. Some typographical and grammatical errors in the abstract include, line 11 “researches”; line 12 “works”; line 13 “….profiles of response to…”.

Response: We have corrected this problem in the revised manuscript. Please refer to line 11-12.

8. Line 28, “…a kind of soil-borne disease…”. It is a soil-borne disease, NOT “a kind-of soil-borne disease”.

Response: This sentence has been corrected. Please refer to line 32.

9. Line 36, not “clubroot resistance crops” but “clubroot resistant crops”.

Response:This sentence has been corrected. Please refer to line 41.

10. Line 32, Brassica rapa????

Response: Pakchoi is categorized as Brassica campestris ssp. Chinensis Makino, whose genome is similar to Chinese cabbage (Brassica rapa. L).

11. Lines 53 and 60, specify the Brassica species! In general, please specify the Brassica (or any plant) species when you are talking about the results from another study.

Response:This issue has been corrected. We have specified the brassica species. Please refer to line 56-77.

12. Line 69—do not begin a sentence with “And”.

Response: We have corrected this problem. Please refer to line 58.

13. Line 79, wrong reference format.

Response:We have corrected reference format. Please refer to line 70.

14. Line 122, NOT “reversely” transcribed.

Response: This issue has been corrected. Please refer to line 128.

15. Line 181, should be “Functional” enrichment.

Response: This issue has been corrected. Please refer to line 191.

16. Line 281, should be “Functional” annotation.

Response: This issue has been corrected. Please refer to line 298.

17. Line 345, “Pearson” correlation coefficient NOT “person”.

Response: We have corrected it in the revised manuscript. Please refer to line 389.

Reviewer #2:

1. I have carefully read the manuscript PONE-D-19-15638 “Integrating long noncoding RNAs and mRNAs expression profiles of response to Plasmodiophorabrassicae infection in Pakchoi (Brassica campestris ssp. Chinensis Makino)” by Zhu et al. The topic of the manuscript is interesting and the integration of lncRNAs to increase our understanding of clubroot disease is timely and novel. However, there are some shortcomings of the manuscript in its present form most notably the absence of an acceptable discussion of the results. More detailed comments can be found below.The English of the manuscript is understandable, but I would strongly advice the authors to consult someone proficient in scientific English for language editing as there are several grammatical errors in the MS. There are many sloppy errors (inconsistent spelling, spelling mistakes, etc) throughout the text. I have pointed out some issues in the detailed comments below, but this list is not exhaustive.

Overall the manuscript would add important findings to our understanding of clubroot disease, but I strongly recommend that the authors considerably revise the manuscript to highlight and describe the very interesting findings they collected.

Response:We appreciated reviewer’s interest in our work. The manuscript has been thoroughly rewritten according to the reviewers’ comments.

2.L 11: delete “Although lots of researches have been conducted during past decades” – grammer errors and IMHO this is not important in the abstract.

Response: This sentence has been deleted.

3.L17: specify which type of enrichment analyses

Response: We have corrected it in the revised manuscript. Please refer to line 16.

4.L19: define which type of interaction relationship, list the most important groups of interactions to provide more info for the reader on the biological relevance of those rather than just providing numbers.

Response: This issue has been corrected and we have also provide more detail on the 15 RNAs in the Results section as the reviewer suggested. Please refer to line 385-401 and S5 table.

5.L21: change “15 clubroot disease related” into “15 P. brassicae mRNAs”, maybe provide more detail on these RNAs.

Response: This issue has been corrected and we have also provide more detail on these RNAs in the Results section as the reviewer suggested. Please refer to line 403-409.

6.Introduction.The introduction informative about the aims and guides the reader to the topic of the manuscript. Some phrases are difficult to understand, especially if the reader is not familiar with the original literature cited, therefore I would advise the authors to carefully edit the introduction.

Response: We agree and have revised the section of Introduction accordingly. Please refer to line 32-100.

7.L29: no need to write P. brassicae in brackets.

Response: We have corrected this problem in the revised version of the manuscript. Please refer to line 34.

8.L60: “between CR and clubroot-resistant (CS) lines” this sentence is not clear, should this be clubroot susceptible (CS)?

Response:We have corrected this problem in the revised version of the manuscript. Please refer to line 65.

9.L60: The sentence starting with “It was also found….” is unclear, please rephrase and explain what you mean with “updated SA function”.

Response: We have corrected this problem in the revised version of the manuscript. Please refer to line 66.

10.L102ff: Sample collection: the authors state that they use race 7 of P. brassice (also check for consistent spelling!). Was this a single spore isolate or a field population showing the characteristics of race 7 using the Williams differential? For all further experiments it will make a huge difference if the experiments were conducted using a single spore isolate or a field population.

Response:This was corrected in the Materials and Methods section. Until 2017, 39 counties and 9 towns of Shanghai had a breakout of clubroot disease and the affected area had reached 2500 hm2. The race 7 of P. brassice was characrilized by field population come from disease nurseries in Qingpu district of Shanghai. Some other scientific institutions have also identified the same result, such as East China University of Science and Technology and Chinese Academy of Agricultural Sciences.

11.L105ff: when were the plants inoculated and with approximately which amount of P. brassicae spores?

Response:We have added this information in the Materials and Methods section of the revised manuscript. Please refer to line 111-115.

12.L107: please change “6 weeks after sowing” by the days post inoculation as this will provide more information on the disease progression

Response:We have corrected it in the revised manuscript. Please refer to line 115.

13.L129ff: what does “clean data” and “dirty reads” refer to? Please rephrase this and avoid laboratory jargon. Having a flow chart of the procedure including versions of the software in the Supplement is usually very helpful for such methods sections (also for LncRNA).

Response:We have included this information in the Materials and Methods section of the revised manuscript. Please refer to line 135-140.

14.L172: state the type of mastermix used

Response:We have corrected this problem. Please refer to line 181.

15.L181ff: This section is not fully clear, please provide more information on the statistical analyses (R packages, settings, assumptions etc used in the analyses)

Response:We have included this information in the Materials and Methods section of the revised manuscript. Please refer to line 134-173.

16.L188: The sentence “the p value…” is incomplete, please rephrase

Response: This sentence has been corrected. Please refer to line 198.

17.L193ff: this section is not clear.

Response:Please refer to line204.

18.L193: lncRNAs - please make sure the correct spelling throughout the MS

Response:This sentence has been corrected. Please refer to line 204.

19.L217: use control or C instead of CK, this abbreviation is not intuitive.

Response: We have corrected this problem in the revised version of the manuscript. Please refer to line 231.

20.L223: use “the number of reads after qality filtering (or quality control) were….” Instead of “clean reads”

Response:This sentence has been removed from the revised manuscript. Please refer to line 236.

21.L225: the sentence startin with Q30 is incomplete

Response:We have corrected this problem in the revised version of the manuscript. Please refer to line 239.

22.L293: “these findings indicate that mRNAs were participating in the defence of clubroot” – This sentence is not correct. There is not evidence that the mRNAs themselves are involved in clubrootdefence, rather it can be assumed that the genes/proteins the mRNA codes for are involved in the reaction. Please rephrase.

Response: We have corrected this problem in the revised version. Please refer to line 333.

23.L301: The phrase “previous researches” is odd, rephrase to “previous studies” or similar.

Response:This issue has been corrected. Please refer to line 341.

24.L307: why were the GO terms restricted to “biological process” and why were molecular function and cellular component omitted? Is there a reason for this choice?

Response: We have corrected this problem in the revised version of the manuscript. Please refer to line 347-370.

25.L346: is it really “person correlation coefficient” or should it be Pearson’s

Response: This issue has been corrected. Please refer to line 389.

26. Overall the figure and table capitations could be a bit more informative and descriptive. Not every image can be easily understood, therefore some comments on specific issues: Tables: Please describe what the different values stand for. I suspect that “annotated” means the number of transcripts that were assigned this GO term in the full dataset – or does this refer to the number of transcripts that were upregulated in clubroot tissue and assigned a certain GO term? Significant – is the number of these GO terms that were significant in which respect – significantly up/downregulated? What does “expected” mean and where does this come from? Classic Fisher refers to what – the significance of the GO-term, the up/downregulated GO-term, the transcripts?

Response:We have added more information and description on the figures and tables according to the reviewers’ comments.

27.Fig 2, 5, 6 Heatmap figures: can you provide any biological information to the genes other than the BRA accession? Maybe adding some sort of functional annotation (GO term, gene name, function of the gene,….). Which values are displayed? FPKM values, and if yes were these normalised? DEGs – but then which values were compared to give the values?

Response:We have corrected Figures show the nomalised significant differentially expressed mRNA and lncRNA.

28.Fig3A: what do the numbers on the x-axis refer to? Please describe in the figure capitation or in the image.

Answer: we have added the description on the x-axis in the figure capitation.

29.Fig7: This figure is not very informative in the present form and also its not a network but a series of correlated genes. Please add information and annotations (which genes are we looking at?) to the figure or convert it into a table (which would provide more information on the individual correlations).

Response:We have corrected Figure 7 and added the Supplemental Table3 containing detailed data used to create this figure.

30.Fig S1: please provide a description of the figure. It is nearly impossible to understand this figure the way it is currently presented.

Response: We have added the description of the figure in the figure capitation.

31.Unfortunately, the discussion feels very incomplete especially as the authors present a number of fascinating results. The authors fail to discuss what the findings of the correlation of lncRNAs and mRNAs presented mean for the biology of clubroot disease. The concept of lncRNAs is employed to clubroot for the first time, so there are plenty of factors that can be discussed and described here. Also there is no comparison to other transcriptomics studies of which there are plenty on a multitude of clubroothots, resistant and susceptible interactions, on the intraplant variation etc. Please use this pool of references to discuss the results in a broader context.

Response:We agree that this section, as it was written, was confusing and problematic. We have thoroughly revised this part of the Dicussion and compared our results with other transcriptomics studies to increase the understandiong of clubroot disease. Although lots of lncRNAs have been found, their biological functions remain unclear. We believe that several more studies will be needed to elucidate these issues. .

32.L357ff: the first paragraph of the discussion is mostly results. Please move the description of the KEGG analyses into the results section, where only GO terms are described currently. Many of these processes have been identified in previous transcriptomic studies of clubroot, please cite those studies and compare their results to the ones generated in this study.

Response:We have corrected these issues in the revised manuscript according to the reviewers’ comments.

33.L376ff: This information is interesting, but most of the data are not yet available. Therefore the validity of more than half of the second paragraph of the discussion cannot be assessed.

Response: Some quantitative trait loci (QTLs) which are related to clubroot diseases have been identified. The lncRNAs are a group of endogenous RNAs that function as regulators of gene expression and may play an important role in several biological processes of plants. So we want to find some negative or positive correlation of lncRNAs and mRNAs. Most of the 15 mRNAs are belonged to defense response proteins by compared with Arabidopsis. However, little is known about the 16 lncRNAs. We agree that further research on these lncRNAs will provide additional information about their detailed roles in pathogen defense.

Submitted filename: Zhu-Plosone-Response to Reviewers.docx

Click here for additional data file.

25 Oct 2019

Integrating long noncoding RNAs and mRNAs expression profiles of response to Plasmodiophorabrassicaeinfection in Pakchoi (Brassica campestris ssp. chinensisMakino)

PONE-D-19-15638R1

Dear Dr. Yuying Zhu,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Serena Aceto, Ph.D.

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: I Don't Know

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors have addressed my comments.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

21 Nov 2019

PONE-D-19-15638R1

Integrating long noncoding RNAs and mRNAs expression profiles of response to Plasmodiophora brassicae infection in Pakchoi (Brassica campestris ssp. chinensis Makino)

Dear Dr. Zhu:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr Serena Aceto

Academic Editor

PLOS ONE