Tissue-Specific Transcriptomic Responses and Viral Accumulation in Lily Cultivars Infected with Cucumber Mosaic Virus

1 Introduction

Cucumber mosaic virus (CMV) is among the most widespread plant pathogens worldwide, infecting over 1200 plant species across more than 100 families [1,2,3]. Beyond its impact on food crops, CMV also threatens ornamental horticulture, where outbreaks have been reported in lilies (Lilium spp.), chrysanthemums, tulips, and many other floricultural plants [4,5,6]. In these hosts, CMV causes mosaic and mottling of leaves, chlorosis, vein clearing, deformation of leaves and flowers, and general growth stunting. In lilies specifically, infections lead to severe foliar mosaics, distortion of new leaves, necrotic spotting, and deterioration of bulb quality [7]. Such symptoms reduce the aesthetic value and marketability of both cut flowers and bulbs, with reported yield losses ranging from 20% under typical conditions to over 50% in severe outbreaks [5].

Plant viruses like CMV extensively reprogram host gene expression, altering cellular pathways to favor infection or counteract defenses [8]. For example, Arabidopsis RCY1 confers resistance to a yellow strain of CMV, and this resistance is accompanied by elevated defense gene expression and salicylic acid (SA) accumulation [9]. In contrast, mutations in essential host factors like eIF4E can confer CMV resistance by impeding viral protein translation [10]. To capture such host responses, high-throughput transcriptomic tools (especially RNA sequencing) are now widely used to profile global expression changes during infection and identify virus-responsive genes and pathways. RNA-seq studies in model plants have illuminated key CMV-responsive pathways. For instance, in tobacco, hundreds of genes are differentially expressed upon CMV infection, notably those linked to reactive oxygen species, SA signaling, and plant-pathogen interactions, which are strongly activated in a resistant cultivar [11]. Similarly, transcriptomic analysis in Arabidopsis revealed broad modulation of defense- and hormone-related genes by CMV [12]. Notably, such approaches are now extending to ornamental crops: transcriptome profiling in CMV-infected chrysanthemum identified stress- and hormone-related genes modulated by the virus [4], and an RNA-seq study in lily uncovered a NAC transcription factor associated with CMV resistance [13]. These findings underscore the broad utility of RNA-seq in elucidating plant-virus interactions. In virus-inoculated leaves (the primary source tissues), plants often exhibit a sharp depression in photosynthetic activity alongside heightened immune signaling. Many viruses actively induce a “sink” metabolism in infected source leaves, repressing photosynthesis while diverting resources toward defense and viral replication [14].

Correspondingly, strong defense responses are triggered in these leaves (e.g., pathogenesis-related proteins and hormone-mediated pathways), as the plant attempts to contain the infection. In contrast, lily bulbs serve as nutrient-storage tissues and sinks for photoassimilates, their infection has more prolonged and pernicious consequences. Because lilies are propagated clonally through bulbs, any virus reaching the bulb can persist and carry over to the next generation [15]. Indeed, systemic virus invasion of sink tissues leads to deterioration of bulb quality and yield losses (often >50% in severe cases) [5]. This makes infection of the bulb a major economic concern, as it not only stunts long-term plant vigor but also threatens the production of healthy, marketable bulbs and flowers. Despite these implications, to our knowledge no study has yet systematically compared virus-induced transcriptomic changes between source leaves and storage tissues in lily. Previous RNA-seq analyses in ornamental hosts have focused only on leaves [13], leaving a significant knowledge gap regarding how sink tissues like bulbs respond at the molecular level to virus stress.

Reports in lilies indicate variation in CMV impact across species or hybrids, and some germplasm shows notable tolerance, for example, Lilium regale has been characterized as virus-resistant and yielded CMV-responsive defense regulators via RNA-seq [13,16]. However, formal comparative studies at the cultivar level within commercial Lilium remain scarce. Most assessments in lilies have emphasized symptom observation or presence/absence diagnostics (ELISA/RT-PCR, multiplex RT-PCR) rather than molecular dissection of tolerance mechanisms [5,15,17]. This evidence gap leaves breeding programs with few robust molecular markers for CMV tolerance, slowing progress toward resistant cultivars. Moreover, defense signaling that is excessively activated can aggravate disease through cell-death pathways, whereas moderated responses may support tolerance, underscoring the need for integrated genotype- and tissue-level analyses of CMV responses in lily.

In this study, we integrate comparative transcriptome profiling with quantitative measurements of viral accumulation. Two Lilium cultivars with contrasting CMV responses were analyzed, and gene expression was profiled in source leaves and storage bulbs at a later stage of infection alongside mock controls. In parallel, leaf virus titers were tracked by ELISA, enabling linkage between leaf viral dynamics and the tissue-resolved transcriptomes. This combined approach jointly compares tissue- and cultivar-specific molecular responses to CMV and connects them to virus accumulation, yielding candidate defense pathways and genes that can guide molecular marker development and breeding strategies for CMV-tolerant lilies.

2 Materials and Methods

2.1 Plant Materials and Virus Inoculation

In vitro-propagated plantlets of Lilium ‘Cancun’ and ‘Connecticut King’ were supplied by the tissue culture facility, National Institute of Horticultural and Herbal Science (Republic of Korea). The Cucumber mosaic virus (CMV) Teaan strain was propagated in Nicotiana tabacum (Xanthi-nc) and mechanically inoculated onto test plants at the three-leaf stage. Virus-infected leaf sap was extracted in 0.1 M sodium phosphate buffer (pH 6.9) and applied to Lilium leaves using carborundum dust to facilitate infection. Each plant was inoculated twice to ensure successful viral transmission, while mock-inoculated plants received only buffer without the virus. Symptom development was monitored, and leaf and bulb samples were collected at 0 and 28 days post-inoculation (dpi) for transcriptomic analysis. The 0 dpi samples represent baseline tissues before systemic infection was established and were used as a reference state in which CMV was not yet detectable by DAS-ELISA. The 28 dpi time point was used as a late systemic stage, when systemic mosaic symptoms and growth reduction were clearly visible in the inoculated plants in this experiment. An independent cohort grown and inoculated under identical conditions was used for ELISA sampling at 0, 14, and 24 dpi.

2.2 RNA Extraction and Sequencing

Leaf and bulb samples from ‘Cancun’ and ‘Connecticut King’ were collected at 0 and 28 days post-inoculation (dpi). Total RNA was extracted using the Qiagen RNeasy Plant Mini Kit (Qiagen, Germany) according to the manufacturer’s protocol. RNA quality and concentration were assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific) and an Agilent 2100 Bioanalyzer (Agilent Technologies). RNA integrity was determined using the RNA Integrity Number (RIN), and only samples with RIN values ≥7.0 were used for sequencing. Libraries were prepared with the Illumina TruSeq RNA Sample Preparation Kit (Illumina, USA), and sequencing was performed on an Illumina HiSeq 2000 platform, generating 100 bp paired-end reads. The resulting raw sequence reads were subjected to quality control filtering before assembly.

2.3 De Novo Transcriptome Assembly, Mapping, and Exploratory Differential Expression

The Illumina adapter sequence was removed from the sequencing data using Cutadapt, and read preprocessing was performed with the SolexaQA package, employing DynamicTrim and LengthSort [18]. De novo transcriptome assembly was performed using Velvet (v.1.1.06) and Oases (v.0.2.01) software based on the de Bruijn graph algorithm [19,20]. Trimming resulted in reads with a mean length of 88.57 bp across all samples and a minimum length of 25 bp. The final assembled transcripts, obtained by optimizing the k-mer and removing duplicate sequences, were used for downstream analyses [21]. Cleaned reads from each sample were then mapped to the assembled sequences using Bowtie2 [22], and a gene-by-sample count matrix was generated from the alignment output.

Since this study is based on a legacy dataset in which only one RNA-seq library was available per condition (cultivar × tissue × time), count data were normalized primarily to correct for differences in sequencing depth using the DESeq framework [23], and fold changes between conditions were estimated from these normalized counts. Genes were classified as differentially expressed within this dataset using a false discovery rate (FDR) threshold of 0.05 with the Benjamini-Hochberg procedure as a heuristic filter. In the absence of biological replication, the resulting DEG lists and associated FDR values should be interpreted as exploratory rankings of expression patterns rather than as formally inferential statistics from a replicated design.

2.4 Functional Enrichment Analysis

Assembled transcripts were functionally annotated by comparing them with gene sequences in the euKaryotic Orthologous Groups (KOG), NCBI non-redundant (NR), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases using the BLASTX algorithm with an e-value threshold of 1 × 10−10. Gene Ontology (GO) annotation was performed by assigning functional classifications based on sequence similarity to proteins in the GO database, with an e-value cutoff of 1 × 10−10 and selection of the best hits [24].

For each comparison, up- and down-regulated gene sets defined within this dataset (Section 2.3) were analyzed separately for functional enrichment. Enrichment of GO terms and KEGG pathways was evaluated using a hypergeometric test with Benjamini-Hochberg correction for multiple testing (FDR). In line with the non-replicated design of the RNA-seq experiment, these enrichment analyses are used here as descriptive tools to summarize over-represented functional categories among the exploratory DEG lists, and the associated FDR values serve as heuristic ranking criteria rather than as formally inferential statistics.

GO enrichment analysis was conducted using Blast2GO [21], with GO terms classified into Biological Process (BP), Cellular Component (CC), and Molecular Function (MF) categories. FDR values (≤0.05) are reported to highlight categories that appear most over-represented in this dataset. The number of assembled transcripts assigned to each GO term was counted using in-house scripts from SEEDERS Co.

2.5 Data Analysis and Visualization

A gene × sample count matrix was assembled from the RNA-seq output for two cultivars (‘Cancun’, ‘Connecticut King’), two tissues (leaf, bulb), and two time points (0 and 28 dpi), yielding eight samples in total. A separate sample metadata table encoded sample, cultivar, tissue, and time. Raw counts were normalized using the DESeq “median-of-ratios” size-factor method: for each gene, we computed the geometric mean across samples; per-sample size factors were obtained as the median of counts/geometric-mean ratios; and counts were divided by these size factors. A log2 transformation (log2(normalized + 1)) was then applied to stabilize the dispersion. To emphasize structure driven by high-variance transcripts, we selected up to the top 5000 most variable genes based on the variance of the log2-normalized matrix.

PCA analysis PCA was performed on the (genes × samples) matrix transposed to samples × genes using scikit-learn (PCA with two components). The proportion of variance explained by PC1 and PC2 is reported in the axis labels. Points are colored by tissue and shaped by cultivar. For paired conditions (same cultivar × tissue at 0 and 28 dpi), vectors (0 → 28 dpi) were drawn to visualize infection-driven trajectories; time labels (“0 d”, “28 d”) were placed near each point with small offsets to avoid overlap. Pearson correlation coefficients were computed from the log2-normalized expression matrix for all sample pairs and displayed as a heatmap, with samples ordered as in the PCA plot. Because only one RNA-seq library was available per condition, the PCA and correlation heatmap provide a descriptive ordination and summary of sample relationships rather than inferential statistics. Analyses used Python 3 (NumPy, pandas, scikit-learn, SciPy) and Matplotlib.

GO/KEGG enrichment analysis For GO/KEGG enrichment visualization, up- and down-regulated gene sets defined within this dataset (Section 2.3) were first obtained from cultivar-specific contrasts (28 dpi vs. 0 dpi). Genes were classified as Up or Down when |log2 fold-change| ≥ 1.0; FDR values from the DESeq framework were used as a heuristic ranking criterion where available, rather than as formal significance thresholds. For tissue-level summaries, DEG lists were then combined across cultivars (Cancun ∪ CK) within each direction (Up/Down). Enrichment statistics for GO terms and KEGG pathways were computed as described in Section 2.4 (one-sided Fisher’s exact test with Benjamini-Hochberg FDR correction). Terms were ranked by FDR and, secondarily, by overlap size and GeneRatio. To aid interpretation in this exploratory dataset, for each panel (Leaf-Up, Leaf-Down, Bulb-Up, Bulb-Down) up to 10 terms with the lowest FDR were displayed; in cases where very few terms met a strict cutoff, an FDR ≤ 0.10 was used only to select terms for plotting. GeneRatio is defined as k/term_size, where k is the number of genes in the DEG set intersecting the term. Bubble area is proportional to k (overlap), and color encodes −log10(FDR). These plots are intended to summarize functional categories that appear over-represented among the exploratory DEG lists rather than to provide formal hypothesis tests.

Heatmap For the marker-focused heatmap, four log2 fold-change contrasts were used: ‘Cancun’ leaf (28 dpi vs. 0 dpi), ‘Connecticut King’ leaf (28 dpi vs. 0 dpi), ‘Cancun’ bulb (28 dpi vs. 0 dpi), and ‘Connecticut King’ bulb (28 dpi vs. 0 dpi). Candidate genes were identified by keyword matching of functional annotations (TAIR symbol/description, GO, and KEGG orthology) across five categories: defense/salicylic acid (SA), jasmonate/ethylene (JA/ET) signaling, WRKY/NAC transcription factors, reactive oxygen species (ROS)/antioxidant processes, and photosynthesis. Within each category, genes were ranked by across-contrast variance and up to eight top genes were retained; remaining slots (up to a maximum of 40 genes) were filled by the highest-variance genes to ensure coverage of additional responsive candidates. The heatmap was generated from the four log2 fold-change values per gene after row-wise z-scaling and clipping to the range [−2.5, 2.5]. Rows were clustered using hierarchical clustering (average linkage, Euclidean distance). Because these contrasts are derived from a non-replicated design, the heatmap should be interpreted as a descriptive summary of candidate marker behavior rather than a statistically validated expression signature.

2.6 CMV Quantification by DAS-ELISA

Plants used for ELISA were grown and inoculated with the CMV-Teaan strain as described in Section 2.1. Leaf and bulb tissues of ‘Cancun’ and ‘Connecticut King’ were collected at 0, 14, and 24 days post-inoculation (dpi). The 0 dpi samples were taken immediately after mock or CMV treatment and represent a pre-systemic baseline stage; at this time point, CMV absorbance values remained at background levels and virus was not detectable above the DAS-ELISA threshold in either cultivar.

For each condition, approximately 100 mg of tissue was homogenized in the kit extraction buffer at a 1:10 (g/v) ratio and clarified by brief centrifugation. Virus concentration was measured by double-antibody sandwich ELISA (DAS-ELISA) using a commercial CMV kit (Agdia, Elkhart, USA), following the manufacturer’s instructions. Aliquots of 100 μL per well were loaded, and absorbance was read at 405 nm (OD405) on a microplate reader. Mock-inoculated tissues and buffer blanks served as negative controls, and a CMV-positive control was included on each plate. Biological replicates were assayed with technical duplicates.

3 Results

3.1 Global Separation of Transcriptomes by Tissue and Time

Principal component analysis (PCA) of the eight RNA-seq libraries revealed a strong tissue effect (leaf vs. bulb) along PC1, which explained 47.7% of the variance (Fig. 1A). For both cultivars, the 0 → 28 dpi vectors progressed largely within the tissue clusters, indicating infection-driven shifts that were smaller than the baseline leaf-bulb difference. PC2 (25.9% variance) further separated samples by time within tissues. Consistently, the sample-sample Pearson correlation heatmap showed high within-tissue similarity and lower cross-tissue similarity (Fig. 1B). Together, these global patterns in this dataset indicate that tissue identity largely structures the transcriptome landscape, with CMV treatment producing tissue-anchored shifts across cultivars. We note that one library per condition was sequenced, so these ordinations are descriptive summaries rather than inferential statistics.

Figure 1: PCA of RNA-seq and sample-sample correlation heatmap. (A) Principal component analysis (PCA) of leaf (green) and bulb (orange) samples from Lilium cultivars ‘Cancun’ (●) and ‘Connecticut King’ (▲) at 0 and 28 days post-inoculation (dpi). Arrows connect the same cultivar × tissue pair from 0 → 28 dpi; labels indicate time points. PC1 and PC2 explain 47.7% and 25.9% of the variance, respectively. RNA-seq was performed with one library per condition; this is a single-replicate, descriptive ordination rather than an inferential analysis. (B) Pearson correlation matrix calculated from log₂-transformed, DESeq-style normalized expression for the same eight samples shown in (A). Higher correlations (yellow) indicate similar global transcriptional profiles; lower correlations (purple/blue) indicate dissimilarity.

3.2 Functional Programs Modulated by CMV Infection Differ between Leaves and Bulbs

Across tissues, the exploratory GO/KEGG analyses suggested distinct functional signatures associated with CMV treatment (Fig. 2 and Fig. 3). In leaves, up-regulated DEGs were associated with GO terms related to carbohydrate catabolism and translation, including polysaccharide catabolic process, ribosome and structural constituent of ribosome, together with stress- and defense-related categories (defense response, response to biotic stimulus, oxidoreductase activity, heme binding). Concordantly, KEGG highlights in this dataset included JA/ET-linked regulators (transcription factor MYC2, jasmonate ZIM domain-containing protein) and primary/translation modules (beta-amylase, replication factor A1, small subunit ribosomal protein S7, histone H4, L-3-cyanoalanine/cysteine synthase, phosphoserine aminotransferase). By contrast, leaf down-regulated sets appeared enriched for cell-wall remodeling and transport together with mobile genetic element/replication features, with GO terms such as xyloglucan:xyloglucosyl transferase activity, cell wall, carbohydrate metabolic process, apoplast, RNA-directed DNA polymerase activity, RNA-templated DNA biosynthetic process, and DNA integration. KEGG mirrored these trends with decreases in xyloglucan:xyloglucosyl transferase (including TCH4), aquaporin PIP, 1-aminocyclopropane-1-carboxylate synthase, HSP20 family protein, EREBP-like factor, calcium-binding protein CML, and ribose-5-phosphate isomerase A.

In bulbs, the dominant up-regulated signal pointed to chromatin/replication and RNA metabolism. Top GO terms in the exploratory enrichment analysis comprised nucleosome, nucleosome assembly, RNA-directed DNA polymerase activity, RNA-templated DNA biosynthetic process, RNA-DNA hybrid ribonuclease activity, and RNA binding, together with polysaccharide catabolic process and beta-amylase activity. KEGG likewise featured histones (histone H3/H4), ribosomal protein L21e, S-adenosylmethionine synthetase, lipoxygenase, EREBP-like factor, pectate lyase, glycerol-3-phosphate acyltransferase, and chalcone synthase, suggesting chromatin-associated remodeling with concurrent lipid/phenylpropanoid and cell-wall-modifying activities in this dataset. Bulb down-regulated DEGs were comparatively broad, with GO terms capturing ATP/nucleotide binding and kinase/phosphorylation processes (ATP binding, nucleotide binding, protein kinase activity, protein phosphorylation), together with cell-wall/carbohydrate metabolism and apoptotic process. KEGG decreases included components of translation and ribosome biogenesis/trafficking (elongation factor 2, midasin), protein homeostasis (heat shock 70 kDa protein 1/6/8), DNA replication/repair (replication factor A1), regulatory proteins (plant MYB, transformation/transcription domain-associated protein), cell-wall enzymes (xyloglucan:xyloglucosyl transferase), specialized metabolism (indol-3-yl-methylglucosinolate hydroxylase), and ion transport (P-type Ca2+ transporter type 2C).

Together, these descriptive patterns suggest that, in this dataset, leaves tend to emphasize immune and hormone-related (JA/ET-linked) programs while dampening selected cell-wall and water-transport modules, whereas bulbs show a coordinated chromatin-replication/RNA-associated program with selective suppression of kinase/transport and protein homeostasis functions. Full GO/KEGG term names and FDR values are provided in Supplementary Table S1, and should be interpreted as exploratory summaries rather than formally tested effects.

Figure 2: GO term enrichment of DEGs at 28 dpi versus 0 dpi. Dot-plots show over-represented Gene Ontology (GO) terms for the union of DEGs across cultivars within each tissue. The x-axis indicates GeneRatio (overlap/term_size). Bubble size encodes the number of overlapping DEGs, and color encodes the −log10(FDR) (BH corrected Fisher’s exact test). For readability, term labels are attached to bubbles (y-axis tick labels hidden) and automatically adjusted to avoid overlap; labels are clamped to remain inside axes. Up to 10 terms/panel are displayed (FDR ≤ 0.05, relaxed to 0.10 if fewer than 10 terms are available). Background: all genes with GO annotation in the dataset.

Figure 3: KEGG pathway enrichment of DEGs at 28 dpi versus 0 dpi. Dot-plots based on KEGG Orthology (KO) annotations. x-axis shows GeneRatio; bubble size indicates overlap; color denotes −log10(FDR). Up to 10 pathways/panel are shown (FDR ≤ 0.05, relaxed to 0.10 when KO coverage is sparse; minimum term size = 5). Labels are placed directly on bubbles and constrained to remain inside the plot area. Background: all genes with KO annotation.

3.3 Tissue-Divergent Response Modules Resolved by Marker Genes

To visualize gene-level trends underlying the global patterns, we compiled a focused panel (≤40 genes) spanning defense/SA, JA/ET signaling, WRKY/NAC transcription factors, ROS/antioxidant enzymes, and photosynthesis. Row-wise z-scaled log2 fold changes (28 dpi vs. 0 dpi) revealed a pronounced tissue bias in this dataset that was consistent across cultivars (Fig. 4).

Figure 4: Heatmap of selected genes (≤40) spanning defense/SA, JA/ET, WRKY/NAC TFs, ROS/antioxidant, and photosynthesis categories. Each column shows the 28 vs. 0 dpi log₂ fold change within a tissue and cultivar (Cancun_leaf, CK_leaf, Cancun_bulb, CK_bulb). Values were row-wise z-scaled and clipped to [−2.5, 2.5]; rows were hierarchically clustered (average linkage, Euclidean distance). Warmer colors indicate higher expression at 28 d relative to 0 d.

In leaves, both ‘Cancun’ and ‘Connecticut King’ showed preferential up-regulation of the JA pathway and a stress-response module. Canonical biosynthetic and signaling components (AOS/CYP74A/DDE2 and OPR1/OPR2) were among the genes with the largest positive fold changes, together with multiple WRKY family transcription factors (e.g., WRKY40/51/22/14) and glutathione-S-transferase-type detoxification markers (e.g., GSTU/ERD). These increases were accompanied by elevated expression of pathogenesis-related markers (PRXR-type entries) and auxin/tryptophan-pathway enzymes (ASA1/BIG/CRM1), a combination that is consistent with mobilization of immunity-linked metabolism in infected foliage within this dataset. The pattern appeared particularly pronounced in CK leaves, where the entire JA/WRKY/GST block shifted toward higher expression.

In the same marker panel, bulbs tended to show the opposite tendency. Compared with 0 dpi, 28 dpi bulbs displayed stronger induction of photosynthesis- and chloroplast-associated transcripts (LHCA/LHCB antenna proteins, CAB family, HCF107) together with chlorophyll-cycle and redox enzymes (CSD2) and several ethylene/NAC regulators (ERF1/ERF8 and NAC083/VNI2). MYC2-related entries were also elevated in bulbs, pointing to a remodeled hormone-associated network in this storage tissue in the context of CMV treatment.

The separation of “leaf-defense” versus “bulb-photosynthetic/organellar maintenance” signatures was most pronounced in CK, but the same directionality was evident in ‘Cancun’ as well. Overall, the gene-level map in Fig. 4 is consistent with the enrichment view: in this exploratory analysis, leaves tend to emphasize immune and JA/ET-linked programs whereas bulbs favor organellar and metabolic maintenance with selective hormone-transcription factor engagement (compare with the GO/KEGG summaries in Fig. 2 and Fig. 3).

3.4 Virus Accumulation Dynamics in Leaves Differ by Cultivar

Leaf DAS-ELISA revealed distinct CMV trajectories between the two cultivars (Fig. 5). In ‘Cancun’, mean leaf OD405 values rose monotonically across the series (0 → 14 → 24 dpi), indicating continuous accumulation through the sampling window. In ‘Connecticut King’ (CK), leaf titers displayed a transient maximum at 14 dpi followed by a decline at 24 dpi. These patterns were reproducible across the biological replicates measured (technical duplicates per replicate were concordant), and the trajectories are displayed for each cultivar as separate panels in Fig. 5.

To relate the serology to the transcriptome patterns, we compared the leaf ELISA trajectories with the 28 dpi expression profiles. In CK, the early rise and subsequent decline in ELISA signal are compatible with the stronger JA-/WRKY-associated expression pattern observed at 28 dpi in leaves (Fig. 4) and with the down-shift of wall/transport-related terms in the exploratory enrichment analyses (Fig. 2 and Fig. 3). Conversely, the progressive increase in ‘Cancun’ leaves is compatible with the comparatively weaker shift of these candidate defense markers by 28 dpi in the marker panel. Because RNA-seq was sampled only at 0 and 28 dpi and without replication, these correspondences should be viewed as hypothesis-generating rather than as evidence for time-resolved causality. In addition, ELISA was performed on leaves only, so any inferences about sink/bulb behavior are drawn from the descriptive RNA-seq results rather than from serology.

Figure 5: Cucumber mosaic virus accumulation in leaves and bulbs of Lilium cultivars ‘Cancun’ (left) and ‘Connecticut King’ (right) measured by DAS-ELISA at 0, 14, and 24 days post-inoculation (dpi). Symbols represent mean absorbance at 405 nm (OD₄₀₅) ± SEM of biological replicates, each measured in technical duplicate. Asterisks denote p < 0.05 from two-way ANOVA (factors: organ, time) with Tukey’s HSD on simple effects when the interaction was significant (per cultivar).

4 Discussion

4.1 Global Structure of the Transcriptome

In our dataset, tissue identity was the principal determinant of variance: PC1 (47.7%) separated leaf and bulb samples cleanly in both cultivars, and the infection trajectories (0 → 28 dpi) were within-tissue shifts rather than cross-tissue jumps (Fig. 1A,B). Thus, at 28 dpi the CMV treatment effect was superimposed on, but did not eclipse, the strong baseline differences between a photosynthetic source tissue and a storage sink. This global pattern motivates our tissue-resolved analyses and frames the later enrichment and marker-based summaries.

Placing these ordinations in context, plant viruses must traverse cell-to-cell interfaces and vasculature to achieve systemic infection, and their success is constrained by tissue-specific barriers [25,26]. Our observation that infection produced consistent but tissue-anchored shifts is therefore consistent with a scenario in which pre-existing tissue programs shape infection trajectories and systemic progression remains contingent on movement across source-sink interfaces, rather than overriding those programs.

The emphasis on bulbs as a distinct cluster is biologically and economically relevant for lilies. Bulbs are the vegetatively propagated material, and virus persistence in stored bulbs is a recognized driver of carry-over across production cycles; practical management even prioritizes bulb testing during storage as more reliable than leaf testing at floral stages. These industry observations reinforce why separating leaf and bulb transcriptomes is useful for interpreting disease risk in Lilium [5,27].

In addition, most molecular CMV studies in lily have profiled leaves and identified resistance regulators (e.g., LrNAC35, LreEF1A4) [13,16], but have not jointly contrasted source vs. sink programs within one experiment. Our ordinations therefore provide a tissue-comparative framework for lilies by showing that CMV-associated shifts ride on robust tissue identity, motivating the tissue-specific enrichment and marker analyses that follow. At the same time, because only one RNA-seq library was available per condition, these global patterns should be viewed as descriptive and hypothesis-generating, rather than as statistically validated effects from a replicated design.

4.2 Functional Programs That Diverge between Leaves and Bulbs

In leaves, CMV treatment was associated with up-regulation of defense- and hormone-linked functions together with translation/oxidoreductase terms, and with down-regulation of cell-wall remodeling and water/transport categories (Supplementary Table S1 and Fig. 2A,C). The KEGG panels highlight induction of JA-related components (e.g., MYC/JAZ context) and decreases in xyloglucan modifiers and aquaporins (Fig. 3, Leaf Up/Down). These descriptive patterns are in line with current views that jasmonate signaling can contribute to antiviral outcomes and potentiate silencing/defense layers in infected leaves, and they are compatible with the idea that restricting wall-loosening and membrane transport may constrain local cell-to-cell spread toward the vasculature [28,29]. Recent syntheses on hormone-virus crosstalk and movement through plasmodesmata (PD) are consistent with both points [30,31].

In bulbs, the dominant Up sets in the exploratory enrichment analysis involved nucleosome/chromatin, RNA metabolism, polysaccharide catabolism (β-amylase), and chloroplast/assembly factors, whereas Down sets were broader, touching ATP/kinase functions and elements of translation/proteostasis (Fig. 2B,D; Fig. 3, Bulb Up/Down). Induction of β-amylase and carbohydrate-use modules alongside maintenance of photosynthetic machinery is consistent with sink-side carbon reallocation under stress, where starch turnover supports metabolic homeostasis rather than a full-scale immune escalation [32]. Work on stress/sugar signaling shows that amylases often rise under carbon draw or starvation-like cues, providing a mechanistic link to the bulb enrichment patterns observed here [33,34].

Taken together, the enrichment contrasts in this dataset are consistent with a source-sink split: leaves tend to prioritize immune/JA programs and temper wall/transport processes that would favor spread, whereas bulbs appear to emphasize metabolic maintenance and carbon remobilization with limited up-shift of canonical defense machinery. The proximity of these trends to established mechanisms, WRKY/JA-centered reprogramming in antiviral defense and PD/cell-wall control of movement lends plausibility to our interpretation of the tissue-resolved responses [35,36], but these inferences remain provisional given the non-replicated design.

4.3 Marker-Guided Tissue Responses and Viral Load Dynamics

In the marker panel (Fig. 4), leaves of both cultivars showed coordinated induction of JA-pathway enzymes (AOS/OPR1/2, CYP74A/DDE2) and multiple WRKY transcription factors (e.g., WRKY40/41/51/53), together with GST-type detox markers. Bulbs tended to show the inverse tendency: higher photosynthesis/assembly signals (LHCB/CAB, HCF107) and β-amylase, with comparatively weak induction of canonical defense markers. These tissue-specific signatures are compatible with the ELISA curves (Fig. 5): ‘Connecticut King’ (CK) exhibited a transient leaf peak at 14 dpi followed by decline, whereas ‘Cancun’ showed progressive increases in leaves. Bulb titers were not measured by ELISA; bulb-related inferences are drawn solely from the transcriptome patterns (Fig. 2, Fig. 3 and Fig. 4).

The leaf program aligns with current views that jasmonate signaling can contribute to antiviral defense and potentiate RNA-silencing layers, while WRKY factors often act as positive regulators of antiviral responses. These mechanisms provide a plausible basis for the containment-like pattern observed in CK leaves, but they are not directly demonstrated here [28,29,37]. The concurrent down-trend of wall/transport terms (Fig. 2 and Fig. 3) complements this working model, given that cell-to-cell and systemic movement rely on plasmodesmata and wall remodeling; restricting those processes could reduce spread from source leaves into the vasculature [31,38]. In CMV specifically, the movement protein interacts with cell-wall apoplastic factors such as ascorbate oxidase, underscoring why wall-linked changes in leaves could influence systemic delivery [39].

On the bulb side, enrichment of β-amylase and chloroplast/antenna assembly factors suggests carbon remobilization with maintenance of photosynthetic capacity rather than broad immune escalation. β-amylase up-regulation under stress/carbon draw is well documented in plants, and together with LHCB/CAB it points to a tolerance/maintenance mode that may help preserve storage-tissue function during infection [40,41]. The observation that CK shows declining leaf titers together with a maintenance-skewed bulb program is consistent with limited systemic establishment in our working model, although bulb titers were not directly measured and causality cannot be inferred from these data alone.

Practically, this exploratory marker-titer concordance nominates leaf readouts (WRKY40/41/51/53; AOS/OPR1/2; CYP74A/DDE2) as candidate indicators of early containment capacity, and bulb readouts (β-amylase; LHCB/CAB/HCF107) as candidate indicators of sink-side maintenance. However, these markers should be regarded as hypotheses generated from a non-replicated dataset; their predictive value and robustness will require validation by qPCR and/or replicated RNA-seq across additional cultivars, timepoints, and environmental contexts before practical deployment. Given the clonal propagation of lilies via bulbs and the risk of carry-over of infected planting stock, the tissue-aware marker sets proposed here may nevertheless provide a useful starting point for designing such validation studies [5].

4.4 Tissue-Specific Implications for Breeding and Stock Management

Within the limitations of this dataset, our results suggest two potential levers for breeding and stock management. First, leaves (source tissue) tend to show coordinated induction of JA/WRKY markers with concurrent dampening of wall/transport functions (Fig. 2, Fig. 3 and Fig. 4). In ‘Connecticut King’, this pattern coincides with a transient leaf titer peak followed by decline (Fig. 5), which is compatible with the idea that earlier or stronger leaf responses could contribute to partial containment. For breeding and screening, a compact leaf panel containing WRKY40/41/51/53 and JA-pathway enzymes (AOS/OPR1/2, CYP74A/DDE2) emerges from this study as a promising candidate readout of containment capacity in challenge tests around 7–14 dpi, but its utility remains to be established empirically.

Second, bulbs (sink tissue) emphasize photosynthetic/assembly factors (LHCB/CAB, HCF107) and β-amylase rather than broad immune escalation (Fig. 2, Fig. 3 and Fig. 4). Given that lilies are propagated via bulbs, limiting sink establishment is critical for preventing carry-over. For stock management and quality control, our results support the concept that storage-phase ELISA on bulbs could in principle be paired with a small bulb expression panel (LHCB/CAB/HCF107; β-amylase) to flag lots that maintain tissue function while keeping virus below actionable thresholds. At present, however, this remains a working hypothesis; testing such a scheme would require replicated trials, defined threshold criteria, and evaluation across multiple genetic backgrounds.

Taken together, the tissue split in this exploratory study offers a simple conceptual framework: prioritize genotypes that are able to mount an early leaf JA/WRKY shift and that show limited sink-side establishment under infection pressure, and then ask whether these tendencies are reproducible across environments and genotypes. In practice, favoring CK-like profiles characterized by leaf containment and bulb protection over ‘Cancun’-like patterns of progressive accumulation is an attractive idea, but it still needs to be validated in dedicated breeding and stock-management trials.

Finally, adding early-timepoint qPCR checks on the leaf and bulb marker sets would be an important step toward substantiating the tissue-specific readouts proposed here, and pairing selection with routine bulb ELISA during storage should help to test whether these markers track virus dynamics reliably in real-world propagation pipelines. This recommendation is consistent with current practice in lilies, where post-storage bulb testing is widely considered more reliable than flowering-stage leaf tests, and where serological or RT-PCR diagnostics are standard for stock hygiene [5,17]. In this sense, our study provides a descriptive, tissue-resolved transcriptomic framework that can be used to generate and refine marker candidates for future replicated validation.

5 Conclusion

Across two commercial lilies challenged with CMV, the transcriptome in this legacy dataset separated first by tissue and then by infection state: in leaves, CMV treatment was associated with a shift toward a JA/WRKY-biased defense program and reduced expression of wall/transport-related functions, whereas bulbs emphasized photosynthetic/organellar maintenance and β-amylase-linked carbon use. These tissue-specific tendencies were qualitatively similar across cultivars and were compatible with the virological readouts: ‘Connecticut King’ showed a transient leaf titer peak followed by decline, while ‘Cancun’ accumulated progressively in leaves. Agreement among ordination, exploratory enrichment patterns, marker behavior, and ELISA kinetics is consistent with a simple working model in which stronger or earlier leaf responses may contribute to partial containment and reduced sink-side establishment, although cause-effect relationships cannot be inferred from these data.

From a practical perspective, our findings suggest potential levers for breeding and stock hygiene in Lilium, while remaining provisional. The leaf marker set highlighted here (WRKY40/41/51/53; AOS/OPR1/2; CYP74A/DDE2) and bulb-associated candidates (β-amylase; LHCB/CAB/HCF107) emerge as promising indicators of containment capacity and sink maintenance that warrant further testing, for example by qPCR and/or replicated RNA-seq in additional cultivars and at earlier timepoints (e.g., 7–14 dpi). In the longer term, such markers could be combined with storage-phase bulb testing (ELISA or RT-PCR) as part of tissue-aware quality control, with the goal of favoring CK-like profiles (leaf containment with limited sink establishment) over patterns resembling ‘Cancun’. At present, however, these applications should be viewed as hypotheses to be evaluated in dedicated validation trials. Overall, our tissue-resolved CMV transcriptome analysis provides a descriptive framework and a set of candidate markers that can guide the design of future replicated studies aimed at improving selection and stock hygiene in lily production.