Molecular responses of arbuscular mycorrhizal fungi in tolerating root rot of trifoliate orange

Arbuscular mycorrhizal fungi (AMF) enhance plant disease resistance, while the underlying mechanisms in the molecular levels are not yet known. In this study, five-leaf-old trifoliate orange seedlings were inoculated with Funneliformis mosseae for 14 weeks and subsequently were infected by a citrus root rot pathogen Phytophthora parasitica by 7 days. The transcriptome results by Illumina HiSeq 4000 revealed that the percentage of Q30 bases reached 92.99% or above, and 29696 unigenes were annotated in a total of 63531 unigenes. 654 and 103 differentially expressed genes (DEGs) were respectively annotated in AMF-inoculated versus non-AMF-inoculated plants under non-infection and infection with P. parasitica, respectively, whilst these DEGs were related to defense mechanisms, signal transduction mechanisms and secondary metabolites biosynthesis. Forty-two genes were functionally annotated as the putative 'defense mechanism', whilst AMF inoculation induced 1 gene down-regulated and 3 genes up-regulated under P. parasitica infection. AMF inoculation stimulated more genes linked to signal transduction mechanism down-regulated than non-AMF plants. Eight genes were involved in secondary metabolites biosynthesis in AMF versus non-AMF seedlings under P. parasitica-infection conditions. Such transcriptome database provided total information in the molecular levels regarding mycorrhizal roles in tolerating Phytophthora parasitica infection.


Introduction
Citrus trees are susceptible to a large number of diseases, including root rot. The main pathogen of citrus root rot in China is Phytophthora parasitica (Tian et al., 2018). Soil arbuscular mycorrhizal fungi (AMF) can build mutual symbionts with most plants (including citrus), viz. arbuscular mycorrhizas (AMs) (He et al., , 2020Wu et al., 2019;Zhang et al., 2020). Many studies showed a crucial function of AMF on increased plant disease resistance (Xie et al., 2019;Zhang et al., 2019;Gao et al., 2020). Ozgonen et al. (2010) inoculated AMF on Arachis hypogaea plants to study its effects on stem rot caused by Sclerotium rolfsii, and found that all the selected AMF species, including Funneliformis etunicatum, F. clarum, F. caledonium, and F. fasciculatum collectively reduced the incidence of stem rot. Inoculation with F. mosseae and Rhizophagus irregularis also mitigated symptoms of root rot in Aphanomyces euteiches-infected pea (Bodker et al., 1998;Slezack et al., 1999). As a result, it seems that AMF has the capacity to mitigate root rot, whereas the underlying mechanisms are not yet known.
Transcriptome analysis with high-throughput sequencing (RNA-Seq), a new molecular biology technique, is designed to elucidate the molecular responses to bacteria-disease occurrence and fungi-disease occurrence (Gao et al., 2016). RNA-Seq technology also identifies responsive genes of pathogen infection from large-scale transcripts of plants, and further analyzes disease resistance mechanisms or pathogenic mechanisms through gene function analysis. Yap et al. (2005) found wound-induced proteins to activate the phosphorylation cascade in mycorrhiza-inoculated-Medicago truncatula by RNA-Seq. Dao et al. (2011) found a chalcone synthase from the transcriptome data of M. truncatula inoculated with Glomus versiforme, which could enhance plant disease resistance. Lambais and Mehdy (2010) found that in transcriptome of Phaseolus vulgaris plants, inoculation with Rhizophagus irregularis up-regulated chitinase and β-1,3-glucanase expression in differentially expressed genes (DEGs) to collectively hydrolyze cell walls of pathogenic fungi (Esquerré-Tugayé et al., 2000;Hu et al., 2017). These results conclude that RNA-Seq could reveal the potential mechanisms at the molecular level in disease resistance. However, information regarding AMF-enhanced resistance of root rot in molecular levels is scarce.
The objective of the present study was to establish the transcriptome of trifoliate orange (Poncirus trifoliata L. Raf.) roots after inoculated with an arbuscular mycorrhizal fungus Funneliformis mosseae and a citrus root rot pathogen Phytophthora parasitica, and to identify DEGs and their metabolic pathways after AMF inoculations.

Experimental design
The experiment consisted of a completely randomized block design with the inoculation with or without Funneliformis mosseae and the infection with or without Phytophthora parasitica. The four treatments were (i) the inoculation without F. mosseae and P. parasitica (-F. m-P. p), (ii) the inoculation with F. mosseae and without P. parasitica (+F. m-P. p), (iii) the inoculation with P. parasitica and without F. mosseae (-F. m+P. p), and (iv) the inoculation with F. mosseae and P. parasitica (+F. m+P. p). Each of four treatments was replicated five times.
Three 5-leaf-old trifoliate orange seedlings grown in autoclaved (0.11 MPa, 121 ºC, 2h) sands were planted into 1.6-L plastic pots supplied with 1.5 kg of autoclaved substrates of soil and sand (5:1, v/v). At the same time, mycorrhizal inoculums (90 g per pot) of F. mosseae were applied into the rhizosphere of seedlings as the AMF treatment. In addition, an equal number of autoclaved mycorrhizal inoculums plus 2 mL filtrate (25 µm filter) of the inoculum were mixed with growth substrates as the non-AMF treatment. All the seedlings were grown in a greenhouse for 14 weeks with 721 to 967 µmol/m 2 /s photon flux, 25/19 ºC average day/night temperature, and 75-85% relative humidity.
The pathogenic fungus P. parasitica was freely provided by the Citrus Research Institute, Chinese Academy of Agricultural Sciences. Such pathogenic infection was conducted out according to Li et al. (2014) with slight modification. Before the infection, the P. parasitica was cultured on potato dextrose agar (PDA) at 28 ºC for one week. The 70% alcohol solution was utilized to sterilize the root neck of trifoliate orange seedlings for 10 s. After rinsed with sterilized water, a sterilized needle was used to make a wound, and 5-mm-diameter mycelial plug of P. parasitica was placed onto the wound. The non-P. parasitica infection treatment received sterilized PDA with same diameter. Subsequently, moistly sterilized absorbent cotton covered on the wound. After 7 days of the pathogen infection, the treated seedlings were harvested.

Illumina sequencing and data processing
The extraction of total RNA in the roots was done using an EASY spin Plus plant RNA kit (Aidlab Biotecnolohies CO. Ltd, China). DNAase (Takara Bio. Inc, Japan) was utilized to remove genomic DNA. After checked RNA yield, purity, and integrity, a total of the 12 RNA samples (three replicates (seedlings)/treatment) with 10 μg total RNA per replicate were sent for RNA-Seq by means of Illimina Genome Analyzer at Biomaker (Beijing, China) in 2017. After purified and concentrated polyadenylated mRNAs with dT-conjugated magnetic beads, directional RNA-Seq library was prepared. After reversed transcription from RNA to cDNA, PCR products with 200-500 bp were purified and quantified for sequencing. Whereafter, Illumina HiSeq 4000 platform was used to sequence the cDNA libraries of each treated root. The raw data were gathered by the sequencer. The raw sequence reads were uploaded into NCBI, in which the SRA number is SRR9665367.
The reads containing two N were eliminated, the adaptor was changed according to adapter information, and then the low-quality reads were trimmed to obtain clean data. Trinity software was used to assemble the clean data and acquire unigenes libraries of each sample.

Identification and annotation of DEGs
To identify DEGs, the edgeR software was utilized. Statistical significance indexes were established based on the fold change (|log2FC| ≥ 2) and p-value (p ≤ 0.01). BLAST software (version 2.2.26) was used to compare DEGs with databases of COG, eggNOG4.5, GO, KOG, KEGG, NR, and Swiss-Prot to obtain annotation information. BLAST parameters e-value not greater than 1×e -5 and HMMER parameters e-value not greater than 1×e -10 were selected.

Results
Evaluation and splicing of sequencing data A total of 85.80 Gb of clean data was obtained. The clean data of each samples reached 6.32 Gb, and the Q30 bases percentage was 93% or above ( Table 1). As a result, the sequencing data in this study had good quality, which could meet the need for subsequent analysis. A total of 63531 unigenes were obtained by splicing, with a total length of 53714816 nt and N50 length of 1715 nt (Table 2). Among them, 15449 unigenes were longer than 1000 bp, and the average length of unigenes was 845.49 bp. In short, a high-quality transcriptome assembly database was obtained.

Screening of DEGs
Depending on the relative expression levels of the regulated samples, the DEGs were defined as upregulated and down-regulated genes. Compared with -F. m-P. p treatment, -F. m+P. p treatment induced 99 genes up-regulated and 111 genes down-regulated in roots, while +F. m-P. p treatment induced 630 genes upregulated and 662 genes down-regulated in roots (Table 3). Compared with -F. m+P. p treatment, +F. m+P. p treatment induced 36 genes up-regulated and 151 genes down-regulated in roots, and+F. m-P. p treatment modulated 665 genes up-regulated and 578 genes down-regulated. The abbreviations are the same as in Table 1.

COG classification of DEGs
Compared with -F. m-P. p treatment, -F. m+P. p treatment induced a total of 80 DEGs annotated to 19 COG classifications (Figure 2a). Among the 19 COG classifications, "general function prediction" catched most of the DEGs (18, 23%), followed by "transcription" (11, 14%). Besides, "replication, recombination and repair", "posttranslational modification, protein turnover, chaperones", and "carbohydrate transport and metabolism" all catched 10 DEGs. In addition, 9 DEGs were captured by "secondary metabolites biosynthesis, transport and catabolism", indicating that these genes may be involved in the secondary metabolites biosynthesis of citrus. Additionally, 5 DEGs and 7 DEGs were annotated respectively to "signal transduction mechanisms" and "defense mechanisms", suggesting that these genes may play a role in resisting pathogen infection.

Root DEGs related to defense mechanism
After blasted with a serial of databases, there were 42 unigenes functionally annotated as the putative "defense mechanism" of all differential treated groups in at least one database (Table 4). Compared with nonmycorrhizal plants, mycorrhizal plants induced 24 unigenes down-regulated and 3 unigenes up-regulated under the condition of no pathogen infection. The biggest drop was c42319.graph_c0 (-6.85) and the biggest rise was c44347.graph_c0 (1.59). Compared with non-mycorrhizal plants, mycorrhizal plants induced 1 unigene (c41443.graph_c0, -3.51) down-regulated and 3 unigenes up-regulated in the presence of the pathogen infection, whilst the biggest rise was c46155.graph_c0 (1.92). Additionally, P. parasitica infection induced 4 unigenes down-regulated and 19 unigenes up-regulated than non-P. parasitica infection treatment in mycorrhizal seedlings. In non-mycorrhizal seedlings, P. parasitica infection induced 8 unigenes down-regulated and 1 unigene up-regulated, indicating that mycorrhizal seedlings performed more actively than nonmycorrhizal seedlings to response to P. parasitica infection.  Table 1 Root DEGs related to signal transduction mechanism After blasted with a serial of databases, lots of unigenes functionally annotated as the putative "signal transduction mechanism" of all differential treatment groups in at least one database, and the top 10 unigenes with the highest expression in each difference treatment group was selected for analysis. The blast results of the 40 unigenes in the Orange Genome Annotation Project database were shown in Table 5. The 40 unigenes were mostly annotated as calcium-binding protein, protein serine/threonine kinase, phosphoinositol kinase, and leucine-rich repeat receptor protein kinase, which were acted as the mediators of signal transmission and play a role in the middle and downstream of signal transduction in the Orange Genome Annotation Project database. In the -F. m-P. p vs -F. m+P. p group, there were 6 unigenes up-regulated and 4 unigenes downregulated in the first 10 expression levels. The most significantly up-regulated gene was c48447.graph_c0 by 1.7 times, and the most significantly down-regulated gene was c43034.graph_c0 by 4.3 times. The annotation information of c47577.graph_c1 and c41969.graph_c1 in the sweet orange database was WRKY19 and MAPK3, respectively, and the expressions of the two genes were collectively down-regulated. In the -F. m-P. p vs +F. m-P. p group, the unigenes in the top 10 expressions were up-regulated for 3 unigenes and downregulated for 7 unigenes. The most significant up-regulation was c48165.graph_c0, up by 3.44 times, and its annotation information in the sweet orange database was WRKY46. The most significant down-regulation was c49730.graph_c0, down by 7.98 times. In the -F. m+P. p vs +F. m+P. p group, the unigenes in the top 10 expressions were up-regulated for 2 unigenes and down-regulated for 8 unigenes. The up-regulated genes were c48593.graph_c1 (1.34) and c44707.graph_c0 (1.29). Among the down-regulated genes, the most significant one was c49730.graph_c0, which was down-regulated by 8.19 times.
In the +F. m-P. p vs +F. m+P. p group, the unigenes in the top 10 expressions were up-regulated for 8 genes and down-regulated for 2 genes. The most significant up-regulation was c35123.graph_c0 increased by 5.02 times. The most significant down-regulation was c28686.graph_c0 decreased by 3.16 times. Interestingly, whether trifoliate orange seedlings were infected by the pathogen or not, mycorrhizal plants induced more down-regulated genes than non-mycorrhizal plants.

Root DEGs related to secondary metabolites biosynthesis
After blasted, lots of unigenes were functionally annotated to the putative "secondary metabolites biosynthesis" of all different treatment groups in at least one database, whilst the top 10 unigenes with the highest expression in each different treatment group was selected for analysis. The blast results regarding the 38 unigenes in the Orange Genome Annotation Project database were shown in Table 6. In the -F. m-P. p vs -F. m+P. p group, there were one gene (c43852.graph_c0, 1.12) up-regulated and nine genes down-regulated in the unigenes according to the top 10 expression levels. The most significant down-regulation was c29552.graph_c0, down by 5.72 times. In the -F. m-P. p vs +F. m-P. p group, the unigenes in the top 10 expression levels all presented down-regulation pattern, and the most significant down-regulation was c29552.graph_c0, down by 7.29 times. In the -F. m+P. p vs +F. m+P. p group, only 8 unigenes showed differential expression pattern, among which 5 genes were up-regulated and 3 genes were down-regulated. The most significant up-regulation was c41941.graph_c0 increased by 2.1 times. The most significant downregulation was c41443.graph_c0 decreased by 3.51 times. In the +F. m-P. p vs +F. m+P. p group, the unigenes in the top 10 expression levels all presented up-regulation pattern, whilst the most significant up-regulation was c41443.graph_c0 increased by 4.73 times.

Discussion
The construction of gene expression profiles through transcriptome sequencing is an important approach for non-model plants lacking genomic sequence data. Lambais and Mehdy (2010) screened out chitinase-and β-1,3-glucanase-related DEGs from the transcripts of Phaseolus vulgaris inoculated with Rhizophagus irregularis. Chitinase and β-1,3-glucanase, as important secondary metabolites in plants, inhibit the infection of pathogenic fungi and improve the disease resistance of plants (Ebrahim et al., 2011). Ward and Weber (2012) used RNA-Seq technology to sequence the resistant strains of raspberry 'Latham' infected by root rot pathogen (Phytophthora rubi), and found that pathogenesis-related protein genes, as well as genes related to tricarboxylic acid cycle and lignin synthesis pathway all presented up-regulation. Cicatelli et al. (2012) and Vangelisti et al. (2018) also demonstrated the positive effects of secondary metabolites induced by AMF in tolerating abiotic stress through RNA-Seq technology. In this work, DEGs related to secondary metabolites biosynthesis presented up-regulation in mycorrhizal plants when suffered the pathogen infection, while DEGs related to secondary metabolites biosynthesis presented down-regulation in non-mycorrhizal plants when suffered the pathogen infection. This indicated that mycorrhizal plants can actively induce the upregulated expression of the genes related to secondary metabolites biosynthesis to respond to the pathogen infection.
Signal transduction is a key link in plant response to external stimulus (Kaur and Gupta, 2005). In this work, compared with -F. m-P. p vs -F. m+P. p, the number and expression level of DEGs annotated as receptor-like protein kinases (RLKs) in +F. m-P. p vs +F. m+P. p differential gene set were more and higher. This kind of protein kinase is a transmembrane protein kinase located on the plasma membrane surface, which can be activated in large quantities under the action of various stress factors, and participate in the signal transduction process to induce plant disease resistance (Xu et al., 2002). In our work, the expression of gene that annotated as mitogen-activated protein kinase kinase kinase 3 in -F. m-P. p vs -F. m+P. p differential gene set was down-regulated by 2.07 times, which was in accordance with the result of Liu et al. (2017) in mango after infected by Fusarium mangiferae. However, compared with -F. m-P. p vs -F. m+P. p differential treatment group, DEGs in +F. m-P. p vs +F. m+P. p were more up-regulated, suggesting that mycorrhizal plants can regulate the up-regulated expression of more genes related to signal transduction to respond to the pathogen infection. However, mycorrhizal inoculation induced the top 10 DEGs related to signal transduction assumed down-regulation of 8 genes and up-regulation of 2 genes when suffered the pathogen infection in comparison with non-mycorrhizal treatment. This might be because mycorrhizal plants selectively mobilize related signal transduction pathways, which needs to be further studied.
Many of the DEGs were linked to cell physiological processes, metabolic processes, stimulus response, catalytic activity, and transcriptional activity in the GO secondary node. Moreover, expression of more genes was induced by AMF inoculation in non-pathogen-infected plants than in pathogen-infected plants. However, from the number of DEGs between -F. m-P. p vs -F. m+P. p and +F. m-P. p vs +F. m+P. p, it found that mycorrhizal plants induced more gene expression than non-mycorrhizal plants when suffered P. parasitica infection. It suggests that although the pathogen infection significantly inhibited the positive effect of AMF inoculation in plants, mycorrhizal plants still regulate various defense mechanisms inside plants to respond to the pathogen infection. In addition, the expression situation of DEGs in secondary metabolites biosynthesis, signal transduction mechanism and defense mechanism in the COG database in the four differentially treatment groups were similar to that in GO secondary node. This further illustrates the important role of AMF inoculation in the response of P. parasitica infection. The annotation results of DEGs in each functional classification of each treatment group also provided important information and data for digging the key genes associated with the defense mechanism of trifoliate orange after AMF inoculation. All these unigenes may become the important information resource for future genetic research in mycorrhizal roles in enhancing tolerance of trifoliate orange in response to biotic stress.