The R. solani AG-3 TB (YC-9) strain was isolated from severely infected tobacco plants in Kuandian County, Dandong City, Liaoning province of China (Wu et al., 2012). The YC-9 strain was activated in potato dextrose agar medium at 28°C for 3 days in the dark. The potato dextrose agar (PDA) medium with R. solani AG-3 TB (6 mm diameter) were inoculated on the acupuncture point and the cotton with sterile water was used for moisturizing. A time course study was performed by acupuncture inoculating R. solani YC-9 on the 5th and 6th leaf of tobacco variety Yunyan 87 (one of the commonly cultivated susceptible variety) at the 9th leaf stage, and harvested at 0, 6, 12, 24, 36, 48, 72 hpi. The center of the acupuncture part was taken, which was drilled with a 1.5 cm diameter punch, then the inoculated leaves were collected and frozen with liquid nitrogen. Four leaves were inoculated per tobacco and each leaf was inoculated four acupuncture points. A total of 105 tobacco plants were inoculated at 0, 6, 12, 24, 36, 48, 72 hpi, among which, 15 tobacco plants were inoculated at each time point and five of them were measured once as one biological replicate. Tobacco leaves inoculated with PDA medium serve at 0, 6, 12, 24, 36, 48, 72 hpi as mock treatment, the method of sample collection was the same as above. At the same time, the fungi were cultivated at PDA (0, 6, 12, 24, 36, 48, and 72 h), which were taken as the fungal control group for subsequent analysis.

Total RNAs were extracted from the fungus inoculated leaf tissues at each time point using TRIzol Reagent (Invitrogen cat. NO.15596026). All the RNA samples were treated with DNase prior to mRNA isolation and sequencing, then the quality was determined using NanodropTM One Cspectrophotometer. And the 1.5% agarose gel electrophoresis using to determine the RNA integrity and using the Qubit 3.0 to quantify the final qualified RNAs. The total RNAs were subjected to stranded RNA sequencing library preparation. Each sample mentioned above was measured one time as one biological replicate. The generation sequencing library was constructed followed the Illumina’s recommendations. Oligo (dT) was used to purify poly (A)-containing mRNA from total RNA. Then the purified mRNA was fragmented and reverse transcribed to cDNAs. The short fragments were connected with adapters at both ends. Thereafter, the adaptor-ligated cDNA was performed using AxyPrep Mag PCR clean-up (Axygen) and recovered the fragments of ∼360 bp. The products were purified and enriched by PCR (11 cycles), and generated the indexed double-stranded cDNA library. The cDNA libraries were analyzed by Agilent 2100 Bioanalyzer and quantified by a Qubit 3.0 Fluorometer (Invitrogen, Carlsbad, CA, USA). Subsequently, the libraries were sequenced by paired-end sequencing under the platform of an Illumina HiSeq 6000 (SeqHealth Co., Ltd, Wuhan, China).

For transcriptomic analysis of R. solani AG-3 TB, raw sequencing data was first filtered by fastp (version 0.23.0) (Chen S. et al., 2018), low-quality reads were removed and the reads with adaptor sequences were trimmed. Then clean and deduplicated data were mapped to the reference genomes of Nicotiana tabacum from https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/715/135/GCF_000715135.1_Ntab-TN90/using STAR software (version 2.5.3a) with default parameters (Dobin et al., 2013) to remove the host transcripts. Then unmapped reads were de novo assembled by Trinity with the default parameters (Grabherr et al., 2011). Sequencing reads were mapped back to the assembled transcripts for assessing the quality of the transcriptome assembly using the Bowtie2 (Langmead and Salzberg, 2012). The longest transcripts of the same genes were screened as the unigenes for annotation and DEG analysis. For functional annotations of the unigenes, the protein databases Nr (NCBI non-redundant protein database), UniProt (universal protein database), Pfam (homologous protein family), eggNog (orthologous groups of genes), GO (Gene Onotology), and KEGG (Kyoto encyclopedia of genes and genomes) were used to infer the amino acid sequences.

The reads per kilobase per million mapped reads (RPKM) were used to compare the levels of differentially expressed genes (DEGs). Each sample in different time points (0, 6, 12, 24, 36, 48, and 72 hpi) was compared with the control. The tool of the EdgeR package (version 3.12.1) was utilized to identify the expression of DEGs (Robinson et al., 2010; McCarthy et al., 2012). The p value cut-off of 0.05 and a fold-change cut-off of 2 was used to determine the statistical significance of gene expression differences. The gene ontology (GO) analysis and Kyoto encyclopedia of genes and genomes (KEGG) enrichment for DEGs were analyzed by KOBAS software (version: 2.1.1) with a p value cut-off of 0.05 to determine the statistically significant enrichment (Wu et al., 2006).

The screening of secretory proteins was performed based on the presence or absence of the predicted coding sequence of the signal peptide, transmembrane domain, ω-sites for glycosylphosphatidylinositol (GPI) anchor, transit peptides to mitochondrion and nuclear localization signal (Zheng et al., 2013). SignalP6.0¹ was used to perform signal peptide cleavage site prediction. Transmembrane helices in the proteins were predicted using TMHMM.² GPI-anchored proteins were identified using PredGPI.³ Proteins located in the mitochondria were determined by TargetP.⁴ Nuclear localization signal was predicted using NetNES.⁵ The proteins that comprise signal peptide cleavage sites and nuclear localization signal, helices without transmembrane domain together with the GPI-anchored proteins were retrieved as secreted proteins. Effector candidates were searched from among the predicted coding sequences of the transcriptome contigs using effectorP.⁶ Localization of probable effectors was predicted using the apoplastP.⁷

To analyze the gene expression of the selected DEGs from each time point, Real-time quantitative PCR (qRT-PCR) was performed using a real-time PCR system Q711 (Vazyme Biotechnology, Nanjing, China) according to the manufacturer’s instruction. The quantitative PCR reaction was carried out in a 20 μl volume containing 1 μl of reverse transcription product, 10 μl of ChamQ Universal SYBR qRCR Master Mix, 0.4 μl of each primer (10 μM) and 8.2 μl of dd H₂O. The reaction conditions for RT-qPCR including three steps (Step 1: 95°C, 30s, Reps1; Step 2: 95°C, 10s, 60°C, 30s, Reps40; Step 3: 60–95°C, increment 0.5°C/5s, Reps1). To verify gene expression of R. solani AG-3 TB during the growth and invasion stage, total cDNA was extracted from total RNA by time course (0, 6, 12, 24, 36, 48, and 72 hpi). The primers were designed by Primer Premier 5 (Supplementary Table 1).

A time course observation and sampling were conducted to clarify the induction of host symptoms and gene expression of pathogenic fungus R. solani AG-3 TB at different infection stages. The results indicated that tobacco leaves with pathogen did not show observable symptoms at 6 and 12 hpi, while the obvious yellow halo can be observed at the inoculated site at 24 hpi. The symptom aggravated with the appearance of wheel pattern after 72 hpi (Figure 1A). The diameters of the necrotic lesions were 0.0036, 0.0152, 0.1709, 0.2633, 0.3939, 0.6434 cm at 6, 12, 24, 36, 48, 72 hpi, respectively (Figure 1B).

Thereafter, an integrated transcriptomic analysis was conducted to globally reveal the crucial genes and pathways involved in R. solani AG-3 TB infection on N. tabacum at different stages, specifically at 6, 12, 24, 36, 48, and 72h after pathogen inoculation. The assembly statistics of R. solani AG-3 TB infection showed a total of 35,415,039 contigs, lengths of N50 and N90 with 1487 and 321 bp, respectively, with 49.06% GC content, of which the GC content maximum was 87.07% and the GC content minimum was 22.52% (Table 1).

The six time points were designated as the early stage (6–12 hpi), middle stage (24–36 hpi), and late stage (48–72 hpi) after R. solani AG-3 TB inoculation. The DEGs change showed that 37,999 DEGs were detected in the early stage, including 18,317 up-regulated genes and 19,682 down-regulated genes after R. solani AG3 TB inoculation. In contrast, the number of DEGs were 43,371, including 28,133 up-regulated and 15,238 down-regulated genes in the middle stage, which comprised largest amounts of DEGs (Supplementary Figure 1 and Supplementary Table 2). These results suggested that the number change of DEGs have the difference during in various hours post inoculation, which may be related to pathogenic factors secreted of R. solani AG-3 TB.

Functional enrichment analysis is an important way to retrieve some significant DEGs for organisms. The Gene Ontology related to three items including biological processes, molecular function and cellular components was analyzed. The results showed that DEGs involved in biological processes including the fatty acid metabolism process, pyridine-containing compound metabolic process and cellular components including proton-transporting V-type ATPase complex and cytosolic were significantly enriched compared with other items in the early infection stage of R. solani AG-3 TB (6–12 hpi) (Supplementary Figure 2A). At 24 hpi after the fungus inoculation, the critical items such as biological processes involved in cellular carbohydrate biosynthetic process, response to stress were enriched, and cell periphery related to cellular components was enriched (Supplementary Figure 2B). In addition, the MAPK cascades involved in biological processes was significantly enriched in the late stage of infection (48 –72 hpi) (Supplementary Figure 2C).

The KEGG pathway analysis indicated that the enrichment pathways of R. solani AG-3 TB were diverse in different infection stages. The results showed that amino sugar and nucleotide sugar metabolism, carbon metabolism, biosynthesis of amino acids and other pathway were significantly enriched in the early infection stage (6 –12 hpi) (Figure 2A). In contrast, the ubiquitin mediated proteolysis pathway was the most significantly enriched item in the middle stage (Figure 2B). Specifically, pyrimidine metabolism pathway was only enriched in the middle infection stage (24 hpi). Furthermore, the amino sugar and nucleotide sugar metabolism, and biosynthesis of amino acid pathways were enriched, while the carbon metabolism pathway was not enriched at 72 hpi stage (Figure 2C).

The synthesis of fungal toxin PAA requires five important enzymes, including shikimate kinase, 3-phosphoshikimate 1-carboxyvinyltransferase (EPSP synthase), chorismate synthase, prephenate dehydrogenase and prephenate dehydratase, of which, shikimate is used as the initial material and the phenyl-pyruvate is the precursor of PAA (Cook et al., 2016; Figure 3A). In the DEGs during R. solani AG-3TB infection, the critical enzymes required for PAA synthesis were retrieved, and their expression levels were gradually increased in different infection stages. A total of 10 DEGs were selected and the expression levels of chorismate synthase, prephenate dehydrogenase were increased in the early and middle stage, while those of shikimate kinase, EPSP synthase and prephenate dehydratase were up-regulated in the middle and late stage of infection (Figure 3B). Specifically, the expression of five enzyme genes in PAA synthesis pathway increased to the highest levels in the middle infection (24 hpi).

Prenyltransferases (DMATSs) and polyketide synthases (PKSs) are two important enzymes involved in the production of secondary metabolites (Slot and Rokas, 2010). The number of DMATSs and PKSs were 11 and 2, and those DMATSs and PKSs genes expression continuously increased from 24 to 48 hpi stage (Figure 3C). The expression levels of most DMATSs were increased at 24 hpi and peaked at 48 hpi. As for PKSs, the expression of DN1756_c0_g1_i1 was increased in the early stage, while the expression of DN14539_c0_g1_i2 increased significantly at 48 hpi (Figure 3C).

Carbohydrate-active enzymes (CAZymes) is a large gene family involved in the construction and breakdown of complex carbohydrates and glycoconjugates, and mainly comprise Glycoside hydrolases (GHs), Glycosyl transferases (GTs), Polysaccharide lyases (PLs), Carbohydrate esterases (CEs) and Carbohydrate-binding modules (CBMs) (Cantarel et al., 2009). The transcriptomic results showed that the gene expression levels of 12, 13, 51, 55, 69, and 53 CAZyme families were gradually up-regulated from the 6 to 72 hpi (Figure 4A). In the early infection stage (6–12 hpi), the secretion of CAZymes was low. Then, the number of CAZymes related DEGs rapidly increased at 24 hpi and peaked at 48 hpi (Figure 5). The quantities of GHs and PLs were slightly higher than other components during R. solani AG-3 TB infection. Results of the heatmap indicated that expression levels of GH45, GH18 and GH16 genes were progressively increased in the early infection stage, while those of PL1, CE4 and CBM12 genes were increased in the middle and late infection stages, respectively (Figure 4C).

Cell wall-degrading enzymes (CWDEs) produced by plant pathogenic fungi, especially those without special penetration, can damage the cell wall polymers (Kubicek et al., 2014). Among these fungi, R. solani can produce CWDEs including pectinase, xylanase, laccase, cutinase and cellulase (Zheng et al., 2013). Here, the pectinase, xylanase, laccase, cutinase and cellulase of R. solani AG-3 TB were retrieved from the RNA-seq (Supplementary Figure 3). The result indicated that expression of 7 CWDEs genes (pectinase and cellulose) were increased in the initial stage of infection (6 hpi), while 79 total genes including pectinase, xylanase, laccase, cutinase and cellulase were significantly increased at 72 hpi (Figure 4B). According to the expression levels of CWDEs, xylanase (TRINITY_DN21213_c0_g1_i1) and cutinase (TRINITY_DN2158_c0_g1_i1) were increased at 24 hpi, while those of the PG1 (TRINITY_DN22192_c0_g1_i1) were up-regulated at 6 hpi (Figure 4C).

When pathogenic fungi infect plants, the respiratory burst oxidase homologs play an important role to increase their pathogenicity (Ghosh et al., 2014). To clarify the gene change involved in respiratory burst oxidase homologs during the interaction between R. solani AG-3 TB and tobacco, the catalase, glutaredoxin, glutathione peroxidase, glutathione S-transferase, copper/zinc superoxide dismutase, and iron/manganese superoxide dismutase were retrieved for further analysis. The results showed that the expression of respiratory burst oxidase homologs increased exponentially in the middle and late stages of infection (Supplementary Figure 4). Additionally, the expression of catalase (TRINITY_DN19585_c0_g1) was significantly up-regulated in the middle stage of infection (36 hpi), and the expression of glutathione S-transferase (TRINITY_DN14304_c0_g1) was increased in the late stage of infection (48 hpi) (Figure 4C).

The secretion of biologically active proteins is a fundamental infection strategy during the interaction between plants and fungi (Lo Presti et al., 2015). Secretomes of fungi play important roles in infection, colonization and pathogenicity (Xia et al., 2020). We adopted the classical secretion pathway (searching protein domains) as well as the apoplastic and cytoplasmic effectors (Effectrop 2.050 and ApoplastP) to retrieve the secreted proteases during the R. solani AG-3 TB infection. A total of 807 potential secretomes were retrieved, and 124 apoplastic effectors and 236 cytoplasmic effectors were predicted. There are few reports on apoplastic and cytoplasmic effectors of fungus, but some report had revealed that cytoplasmic effectors in phytophthora sojae encoding conservative sequences such as RxLR can cause tissue necrosis of plants (Zhang et al., 2015; Sperschneider et al., 2018). A total of 13 cytoplasmic effectors and 28 apoplastic effectors were retrieved, and most of these effectors can be classified in to serine protease (TRINITY_DN17844_c0_g3_i1), eukaryotic metallothionein (TRINITY_DN22450_c1_g5_i2), polysaccharide deacetylase (TRINITY_DN3377_c0_g1_i1), acetylxylan esterase (TRINI TY_DN7106_c0_g1_i1), extracellular metalloproteinases (TRINITY_DN18673_c0_g3_i1), and deuterolysin metalloprotease (TRINITY_DN19525_c0_g1_i1) families (Supplementary Table 3). The expression analysis of pathogenic protease genes showed that the candidate cytoplasmic effector genes were significantly up-regulated at the late infection stage (48 hpi), while the genes of apoplastic effector were enriched significantly in the middle and late infection stages (24, 48 hpi) (Figure 6).

Small cysteine-rich proteins play functional roles in the molecular interaction between fungi and plant (Stergiopoulos and de Wit, 2009), and generally have a classical structural characteristic with a less than 300aa protein and more than 4% cysteine (Yamamoto et al., 2019). Here, a total of 78 potential small cysteine-rich proteins were retrieved from R. solani AG-3 TB data (Supplementary Table 4). It should be noted that most of the small cysteine-rich proteins are unnamed protein products, which still require further investigation for clarification of their functions.

We chose the DEGs with large difference expression and revealed the potential critical roles (CAZymes, toxins and effectors) during R. solani AG-3 TB infection in different infection stages (6 hpi, 12 hpi, 24 hpi, 36 hpi, 48 hpi, 72 hpi). The expression levels of DN3377 (TRINITY_DN3377_c0_g1_i1), DN14681 (TRINITY_DN14681_c0_g1_i2), DN17797 (TRINI TY_DN17797_c0_g2_i1), DN18673 (TRINITY_DN18673_c0 _g3_i1), DN20404 (TRINITY_DN20404_c0_g2_i1), DN20778 (TRINITY_DN20778_c0_g1_i1), DN29708 (TRINITY_DN2 9708_c0_g1_i1), DN21213 (TRINITY_DN21213_c0_g1_i1), DN18070 (TRINITY_DN18070_c0_g1_i4), DN21085(TRIN ITY_DN21085_c0_g1_i1) were verified by qRT-PCR. The result indicated that CAZymes (DN29708, DN3377, DN20778, DN21213, and DN17797) expression were diverse, but increased levels were found for those five genes in the middle and late infection stages. EPSP synthase (DN14681) and chorismite synthase (DN18070) were two important genes involved in PAA synthesis pathway, and their expression started increasing at 24 hpi. The 36 hpi and 48hpi were important infection stages for secretory protein (DN20404, DN18673, and DN21085). These results indicated that the pathogenic genes expression obviously changed and also proved the reliability of transcriptomics during R. solani AG-3 TB infection (Figure 7).