Nicotiana tabacum L. cv. K326 and N. benthamiana plants were grown in the artificial climate chamber that maintained at 25°C (day/night), 16 h/8 h (light/dark) cycles and 65% relative humidity. Potato virus Y (PVY-LN, GenBank ID: JQ971975) was isolated and purified by our laboratory and propagated on tobacco. Tobacco leaves were sprayed twice with H₂O or MgSO₄ with concentration of 240 mg·L⁻¹ at 4–5 leaf stage at 72 and 24 h before phosphate-buffered saline (PBS) solution or PVY inoculation, respectively (PBS solution + H₂O or P + H; PBS solution + magnesium or P + Mg; PVY + H₂O or PVY + H; PVY + magnesium or PVY + Mg). Crude extracts from 1 g of PVY-infected tobacco leaf tissues homogenized with 0.01 mol·L⁻¹ PBS (pH = 7.2) was mechanically inoculated on the surface of tobacco leaves. The inoculated leaves of tobacco at 1 day post inoculation (dpi) and the systemic leaves of tobacco at 3 dpi, 5 dpi, 7 dpi and 9 dpi were harvested for measurement of virus accumulations. The inoculated leaves at 1 dpi and the systemic leaves of tobacco at 3 and 9 dpi were harvested for RNA-seq analysis. Each treatment was performed for three biological replicates with at least nine plants.

About 1 μg of total RNA from each sample was used as input for RNA-seq. The libraries were generated using a NEB Next Ultra RNA Library Prep Kit. The libraries were sequenced on an Illumina HiSeq TM2500 (Biomarker Technologies Co. Ltd., Beijing, China). The clean reads were mapped to the reference genome of N. tabacum.¹ The relative gene expression levels were normalized as fragments per kilobase of transcript per million mapped reads (FPKM). We set the threshold of false discovery rate (FDR) < 0.05 and |log2 fold change| ≥ 1 as DEGs. Gene Ontology (GO) enrichment of DEGs were analyzed by a GOseq R packages based Wallenius non-central hyper-geometric distribution (Young et al., 2011). We used KOBAS software to perform KEGG analyses of DEGs (Mao et al., 2005). The sequencing data were deposited in the SRA database at NCBI with the accession number PRJNA903693.

Weighted gene co-expression network analysis (WGCNA) was used to construct gene co-expression networks (Langfelder et al., 2009). Highly co-expressed gene modules were obtained using the WGCNA v3.1.1 package in R language (Langfelder et al., 2009). A gene expression adjacency matrix was constructed to analyze the network topology with an unsigned type of topological overlap matrix (TOM), a power β of 6, a minModuleSize of 15, and minimum height for merging modules of 0.08036.

We used a previously reported tobacco rattle virus (TRV) vector for virus-induced gene silencing (VIGS) assays (Bachan and Dinesh-Kumar, 2012). The constructions of TRV-based vectors were performed according to Guo et al. (2022). The primers used are shown in Supplementary Table 1. The TRV1 and TRV2 plasmids were transformed individually into Agrobacterium tumefaciens strain GV3101. Agrobacterium cultures carrying TRV1 vector and TRV2 vector were mixed equally in volume with each final concentration of OD₆₀₀ = 0.5, which were then infiltrated into the fifth and sixth leaves of eight-leaf N. benthamiana plants. The upper two non-infiltrated N. benthamiana leaves were mechanically inoculated with PVY or PBS solution after 10 days post infiltration with TRV. The upper two systemically infected leaves were collected at 10 days post PVY inoculation.

Total RNA of tobacco K326 and N. benthamiana leaf tissues were extracted using TRIzol reagent (TIANGEN, Beijing, China). The first-strand cDNA was synthesized using 2 μg of total RNA. The real time quantitative PCR (RT-qPCR) was performed as previously reported (Guo et al., 2022). The expression level of N. tabacum Ntubc2 (AB026056.1) and N. benthamiana NbActin (AY179605.1) was used as an internal control, respectively. The specific primers used in RT-qPCR detection are listed in Supplementary Table 1. All the experiments were performed with at least three independent biological replicates.

Total proteins of tobacco K326 and N. benthamiana were extracted using a Plant Protein Extraction Kit (Solarbio, Shanghai, China). The proteins separated by 12% SDS-PAGE electrophoresis and were transferred to 0.20 μm polyvinylidene fluoride (PVDF) membranes (Sangon Biotech, Shanghai, China), which were then incubated in blocking solution for 1 h (Solarbio, Shanghai, China). PVY CP antibody was used at a dilution of 1: 1000 (Youlong, Shanghai, China). Beta-actin antibody was used at a dilution of 1: 5000 (Proteintech, Chicago, USA). Secondary antibody was used at a dilution of 1: 10000 (ABclonal, Wuhan, China). After incubation with primary and secondary antibodies, membranes were washed twice for 15 min with 1 × TTBS. Finally, the membranes were transferred into ECL solution (Millipore, Billerica, USA) to detect signals by Tanon Chemiluminescence Gel Imager (Tanon, Shanghai, China).

IBM SPSS Statistics 25.0 software (IBM Inc., Armonk, USA) was used for data analysis. The differences among groups were analyzed through two-tailed t test and one-way analysis of variance (Duncan).

To verify the effect of Mg on tobacco anti-PVY infection, we set up four treatments (PBS solution + H₂O or P + H; PBS solution + magnesium or P + Mg; PVY + H₂O or PVY + H; PVY + magnesium or PVY + Mg) at five different time points. The results showed that milder PVY symptoms occurred on PVY + Mg plants compared with that on PVY + H plants (Figures 1A,B). At 9 dpi, the PVY + H plants showed severe stems necrosis and leaves chlorosis, while those PVY + Mg plants showed only slight stem browning (Figures 1A,B). RT-qPCR and Western blot were used to detect the accumulations of genomic RNAs and CP proteins of PVY in the upper two tobacco leaves, and the results were consistent with the severity of our observed symptoms (Figures 1C,D). At 7 dpi, the accumulations of PVY genomic RNAs and CP proteins were decreased by 48 and 64% in PVY + Mg plants, compared with that in PVY + H plants (Figures 1C,D). These results indicated the positive effect of Mg on inhibiting PVY infection in tobacco.

In order to explore the molecular mechanism of Mg regulation on tobacco resistance to PVY infection, we performed four different treatments on tobacco plants at 1, 3, 9 dpi for Illumina RNA-seq, resulting in a total of 36 libraries. The correlation analysis of three replications of 36 samples showed that Pearson’s correlation coefficients were between 0.586 and 0.995, which proved that the correlation between each group of biological replicates was high (Supplementary Figure 1). Each library contained ≥6.11 Gb of clean data with Q30 quality scores ≥92.98%, and CG content percentage between 43.36 and 45.11% (Supplementary Table 2). Mapped the sequencing reads to reference genome of N. tabacum cv. K326, the comparison efficiency ranged from 78.63 to 96.23% (Supplementary Table 3).

To further elucidate the transcriptomic variations of tobacco leaves resistant to PVY infection induced by Mg application, we conducted pairwise comparison of different treatments (i.e., PVY + Mg vs. P + Mg, PVY + Mg vs. PVY + H, P + Mg vs. P + H, PVY + H vs. P + H and PVY + Mg vs. P + H) at 1, 3, 9 dpi. Based on the standard, a total of 11,970 DEGs were identified, of which 429, 27, 287, 78, and 146 DEGs were found in PVY + Mg vs. P + Mg, PVY + Mg vs. PVY + H, P + Mg vs. P + H, PVY + H vs. P + H and PVY + Mg vs. P + H at 1 dpi, respectively (Figures 2A,B). At 3 and 9 dpi, 163 and 92 DEGs were found in PVY + Mg vs. PVY + H (Figures 2A,B). The expressions of a total of 9,708 genes were changed in the PVY + H vs. P + H comparison at 9 dpi, significantly more than that at 1 and 3 dpi (Figures 2A,B). In P + Mg vs. P + H comparison, we found that the DEGs were mainly concentrated in 1 and 3 dpi (Figures 2A,B). In addition, more DEGs in PVY + Mg vs. PVY + H comparison at 3 dpi than that at 1 and 9 dpi (Figures 2A,B). These results indicated that the changes of gene expression in tobacco plants induced by PVY infection were increased with time, while the regulation of Mg on tobacco gene expression was mainly concentrated at 3 dpi during PVY infection.

To further explore the effects of PVY infection on tobacco under Mg treatment, we selected DEGs in PVY + Mg vs. PVY + H at 3 and 9 dpi for GO and KEGG pathway enrichment analyses, respectively. The GO terms in PVY + Mg vs. PVY + H at 3 dpi were mainly enriched in the biological process (BP) terms ‘metabolic process’ and ‘single-organism process’, cellular component (CC) terms ‘cell part’ and ‘cell’, and the molecular function (MF) terms ‘catalytic activity’ and ‘binding’ (Figure 3A; Supplementary Table 4). The results of GO analyses at 9 dpi were basically consistent with those at 3 dpi (Figure 3A; Supplementary Table 4). Subsequent KEGG enrichment analyses showed that DEGs in PVY + Mg vs. PVY + H at 3 dpi were significantly enriched in ‘starch and sucrose metabolism’ (ko00500), ‘fatty acid elongation’ (ko00062), ‘amino sugar and nucleotide sugar metabolism’ (ko00520) and ‘nitrogen metabolism’ (ko00910), while ‘photosynthesis’ (ko00195), ‘oxidative phosphorylation’ (ko00190) and ‘arginine and proline metabolism’ (ko00330) were the top three pathways at 9 dpi (Figure 3B).

To further analyze the gene regulatory network of Mg regulating resistance to PVY infection in tobacco, we performed WGCNA analysis on all the obtained genes. We finally identified 10 different gene regulatory modules containing DEGs ranging from 19 to 4006 (Figure 4A; Supplementary Table 5). The correlations between modules and modules, and modules and treatments were showed (Figures 4B,C). The results showed a strong correlation between green and magenta modules and these DEGs were highly expressed in PVY + Mg at 3 dpi. Top GO and KEGG analyses showed that these DEGs were mainly enriched in ‘energy production and conversion’, ‘photosynthesis’, ‘plant-pathogen interaction’ and ‘protein processing in endoplasmic reticulum’. These two modules mainly included DEGs related to oxidative phosphorylation (NtAS, NtNQOR), photosystem II (NtPsbC), photosystem I (NtPSI-A2) and some chaperones (NtHSP90, NtHSP20). The turquoise and black modules were highly correlated with PVY + Mg at 9 dpi. These DEGs were mainly enriched in ‘inorganic ion transport and metabolism’, ‘carbohydrate transport and metabolism’ and ‘plant-pathogen interaction and respiratory burst’, which were mainly involved in starch and sucrose metabolism (NtSS, NtTPS), nitrogen metabolism (NtrTl-Cyn, NtNRT), Ca²⁺ signal transduction (NtCML44, NtCML35, NtCML41, NtCML36) and ROS scavenging (NtPrx) (Supplementary Tables 6, 7). Through WGCNA analysis, we found that the DEGs in green and magenta modules were involved in the resistance of tobacco to PVY infection regulated by Mg at 3 dpi, while the DEGs in turquoise and black modules were involved in the resistance of tobacco regulated by Mg at 9 dpi.

Carbohydrate metabolism is closely related to other metabolic activities, which provides energy for plant growth and enhances plant stress resistance (Pommerrenig et al., 2018). In starch and sucrose metabolism pathway, we found 83 DEGs at 1, 3 and 9 dpi under different treatments, including five starch synthases (NtSSs), four amylases (NtAMSs), three sucrose-phosphate synthases (NtSPSs), five sucrose synthases (NtSUSs), 15 alpha, alpha-trehalose-phosphate synthases (NtTPSs), 10 UDP-glucuronate 4-epimerases (NtUDPGLEs), three UDP-glucuronic acid decarboxylases (NtUXSs), three 1,4-alpha-glucan-branching enzymes (NtGBEs), one 4-alpha-glucanotransferase (NtGTFB), two cellulose synthase-like proteins (NtGESAs), four galacturonosyltransferases (NtGAUTs), 22 glucosidases (NtGBAs), two endoglucanases (NtEGUs), two beta-D-xylosidases (NtXyls) and two alpha-1,4 glucan phosphorylase L-2 isozymes (NtGPs) (Figure 5A; Supplementary Table 8). At 9 dpi, the expression levels of almost all these DEGs were up-or down-regulated in PVY + H vs. P + H. In PVY + Mg vs. PVY + H, the expression level of NtSPSs, NtGAUT12, NtXyl5 were up-regulated at 3 dpi, and NtGBEs and NtTPSs were down-regulated to some extent at 9 dpi. We identified 13 DEGs associated with fructose and mannose metabolism pathway, of which fructokinase-1 isoform X1 (NtFKs) and mannose-6-phosphate isomerase 1-like (NtMPIs) were up-regulated in PVY + H vs. P + H at 9 dpi, while mannan endo-1,4-beta-mannosidases (NtManAs) were down-regulated (Figure 5A; Supplementary Table 8). In PVY + Mg vs. PVY + H, we found that NtManAs were up-regulated at 3 dpi. In glycolysis / gluconeogenesis pathway, seven DEGs were obtained, including one phosphoglucomutase (NtPGM), four hexokinases (NtHKs) and two triosephosphate isomerases (NtTims) (Figure 5A; Supplementary Table 8). Most of these DEGs were down-regulated in PVY + H vs. P + H at 9 dpi, while unchanged in PVY + Mg vs. PVY + H. There were 28 DEGs related to amino sugar and nucleotide sugar metabolism pathway (Figure 5A; Supplementary Table 8). In PVY + H vs. P + H, UDP-arabinose 4-epimerases (NtUDP-Araps), chitinases (NtCases), most of endochitinases (NtEases) and chitotriosidases (NtCHTs) were up-regulated at 9 dpi, while UDP-D-apiose/UDP-D-xylose synthases (NtUAXs) were down-regulated. These results suggested that PVY infection seriously affected plant carbohydrate metabolism at 9 dpi, while mannose metabolism and sucrose metabolism were regulated by Mg under PVY infection.

Reactive oxygen species (ROS) bursts and allergic reactions induced by Ca²⁺ concentration in intracellular and extracellular are important immune mechanisms in plant response to stresses (Moeder et al., 2019). In this study, we identified 68 DEGs in Ca²⁺ signaling transduction pathway at 9 dpi in PVY + H vs. P + H (Figure 5B; Supplementary Table 8), of which 35 DEGs were up-regulated including 14 NtCAM/CMLs, eight calcium-dependent protein kinases (NtCDPKs), nine calcium-transporting ATPases (NtCa²⁺-pumps), one calmodulin-binding transcription activator (NtCAMTA), one cation/calcium exchanger (NtCCX) and two calmodulin-binding receptor-like cytoplasmic kinases (NtCRCKs), while 13 DEGs were down-regulated including three NtCAM/CMLs, three NtCDPKs, two NtCa²⁺-pumps, three calcium sensing receptors (NtCaSRs), one NtCCX and one CDPK-related kinase (NtCRK), but no DEGs were found at 1 dpi and 3 dpi. In PVY + Mg vs. PVY + H, we found that only NtCDPK29 was up-regulated at 3 dpi, while three NtCAM/CMLs and one NtCa²⁺-ATPase were down-regulated at 9 dpi. These results showed that Ca²⁺ signaling transduction was involved in response to PVY infection and Mg application may inhibit part of Ca²⁺ flow at 9 dpi.

Nitrogen is an essential nutrient for plant growth, and nitrogen metabolism is one of the important ways for plants to resist stresses (Wang et al., 2014). In this study, we found 19 DEGs involved in nitrogen metabolism (Figure 5C; Supplementary Table 8). At 9 dpi, one nitrite reductase (NtNR) and two carbonic anhydrases (NtCAHs) were down-regulated in PVY + H vs. P + H, while two cyanate hydratases (NtrTI-Cyns) were up-regulated. At 3 dpi, only one NtCAH was down-regulated in PVY + Mg vs. PVY + H.

Plant energy metabolism is the basis of life activities (Siqueira et al., 2018). By analyzing all DEGs at 1, 3 and 9 dpi, we identified a total of 80 DEGs related to oxidative phosphorylation (Figure 5D; Supplementary Table 8). At 9 dpi, most of the enzymes related to oxidative phosphorylation were up-or down-regulated to varying degrees (i. e. NtASs, NtACLYs, NtDLDs, NtETFs, NtFNRs, NtFHs, NtMDHs, NtSnRK1s, NtPPases, NtAACs, NtALDHs, NtCACs, NtAKRs, NtMCPs, NtUCPs, NtCSs, NtNDHs and NtGDPDs) in PVY + H vs. P + H comparison. At 1 and 3 dpi, however, the expression levels of these DEGs were unchanged. In PVY + Mg vs. PVY + H comparison, we found that three NtASs related to mitochondrial oxidative phosphorylation reaction, and one NtNDH related to respiratory chain were up-regulated at 9 dpi. These results indicated that PVY infection interfered with plant energy metabolism, and exogenous Mg could improve the activity of related enzymes to a certain extent, thus promoting the energy metabolism of plants.

We randomly selected six genes to verify their expression levels in four different treatments at three time points (Figure 6). RT-qPCR results showed that the expression levels of LOC104121166 (gene_42165) and NtTHI (gene_395) were down-regulated in PVY + Mg vs. PVY + H, while up-regulated in PVY + H vs. P + H at 1 dpi. The expression levels of NtrbcS (gene_61355) were unchanged in PVY + Mg vs. PVY + H at 1 dpi, while down-regulated in PVY + H vs. P + H. The expression levels of NtPI (gene_80129) were down-regulated in PVY + Mg vs. PVY + H and PVY + H vs. P + H at 3 dpi. The expression levels of NtWRKY40 (gene_48638) were unchanged in four different treatments at 3 dpi. The expression levels of LOC102603354 (PB.15726) were unchanged in PVY + H vs. P + H, while up-regulated in PVY + Mg vs. PVY + H at 9 dpi. These results are basically consistent with our transcriptome sequencing results.

Transcriptome sequencing results showed that NtNR and NtPPases that involved in nitrogen metabolism and energy metabolism were down-regulated in PVY + H vs. P + H at 9 dpi. NtTPS and NtGBE related to carbohydrate transport and metabolism were down-regulated in PVY + Mg vs. PVY + H at 9 dpi, while NtCML36 related to Ca²⁺ signaling transduction was up-regulated. In order to further explore the roles of these five DEGs, we selected their homologous genes in N. benthamiana and performed functional analyses through TRV-based VIGS assays. At 10 dpi, we found that NbGBE-silenced N. benthamiana had lower PVY accumulation and weaker symptoms, while NbNR-and NbPPases-silenced plants accumulated higher PVY and showed more serious leaf curl and chlorosis (Figures 7A–D). The PVY accumulation in NbTPS-silenced plants was significantly reduced compared with that in control plants, and the severe yellowing symptoms of these plants may be caused by gene silencing combined with PVY infection (Figures 7A–D). The virus accumulation in NbCML36-silenced plants was almost the same as that in control plants (Figures 7A–D). The silencing efficiency of target genes ranged from 55 to 85% determined by RT-qPCR (Figure 7E).