Abstract Huanglongbing (HLB), also known as citrus greening, is a devastating disease that affects citrus crops. Citrus australis (Australian round lime), a wild citrus species, has been reported to exhibit some resistance to HLB. In this study, we investigated the biochemical and molecular responses to HLB by several F1 hybrids between the Marisol clementine crossed with the Australian round lime (MRL), with a focus on identifying the mechanisms underlying HLB tolerance. A selected HLB-tolerant hybrid (MRL2-12) showed fewer HLB symptoms and maintained a healthy canopy, whereas the other C. australis hybrids exhibited typical HLB symptoms. The MRL2-12 hybrid exhibited the highest chlorophyll content and the least starch accumulation, both of which are important markers for HLB tolerance. Based on the Candidatus Liberibacter asiaticus (CaLas) titer and biochemical analysis, the gene expression patterns of a selected susceptible hybrid (MRL2-11) and the tolerant MRL2-12 hybrid were further analyzed using RNA-seq data to investigate plant defense responses in the context of HLB. The transcriptomic data from the MRL2-11 and MRL2-12 hybrids revealed different responses to HLB, with a set of differentially expressed genes between the tolerant C. australis hybrid and the susceptible hybrid, which were both grown under the same field conditions. These results revealed that the expression of genes related to cellular defense and pathogenesis-related defense mechanisms was significantly upregulated in the tolerant MRL2-12 hybrid compared with the MRL2-11 hybrid. MRL2-12 showed upregulated expression of pattern recognition receptors (PRRs), receptor-like kinases (RLKs), calcium-dependent protein kinases (CDPKs), and cysteine protease proteins, indicating effective defense mechanisms. Comparative genomic analysis identified significant polymorphic variants in MRL hybrids, indicating a genetically diverse background. These findings suggest that early, coordinated activation of immune signaling and physical defense mechanisms, such as cell wall fortification, plays a critical role in HLB tolerance in the C. australis citrus hybrids. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-025-11942-x. Keywords: Citrus australis, Citrus greening, Australian round lime, Hybrid purity, HLB tolerance, NBS-LRR, Transcriptome Background Huanglongbing (HLB), also known as citrus greening, is the most devastating citrus disease. HLB is caused by the gram-negative, phloem-limited, vector-transmitted bacterial species Candidatus Liberibacter asiaticus (CaLas). This bacterium is spread primarily by the Asian citrus psyllid (Diaphorina citri). HLB symptoms typically include leaf mottling, blotchy chlorosis, yellow shoots, and lopsided fruit with aborted seeds [[32]1]. HLB symptoms are considered to occur because of a series of physiological, cellular, biochemical, and molecular disorders in the citrus host [[33]2]. Most commercial citrus cultivars are susceptible to HLB. However, Citrus australis and other native Australian species, such as C. australasica, C. glauca, and C. inodora, along with their hybrids, have been reported to have varying degrees of tolerance to HLB [[34]3, [35]4]. Due to their ability to enhance a tree's tolerance to HLB, citrus breeders are interested in using these cultivars in citrus improvement programs to introgress their disease tolerance traits. []. Identifying and characterizing the factors responsible for conferring tolerance to HLB is crucial for breeding resistant citrus plants. Transcriptomics provide valuable insights into the molecular responses of plants to a particular condition. Comparative transcriptomics of HLB-susceptible and HLB-tolerant hybrids can provide a platform to elucidate the molecular mechanisms underlying HLB tolerance and identify robust candidate genes for developing commercial HLB-tolerant citrus hybrids. Transcriptomic analysis of hybrids of Poncirus trifoliata and C. sunki, which were classified as susceptible, tolerant, and resistant to HLB, revealed a hypothetical model for understanding the genetic mechanisms involved in HLB tolerance. Genetic tolerance was associated with the activation of receptor-related genes, transcription factors, and cell wall fortification processes [[36]2]. Numerous other comparative transcriptome analyses have also been performed on citrus trees infected with CaLas [[37]2, [38]5, [39]6, [40]7, [41]8, [42]9]. These studies provide an in-depth understanding of the responses of tolerant cultivars to infection. These findings suggest that some cellular and defense responses are common among all tolerant cultivars, whereas others are specific to individual cultivars. With no effective cure for HLB currently, breeding and identifying HLB-tolerant citrus varieties has emerged as a promising strategy for controlling this devastating disease. Transcriptomic profiling is an effective tool for understanding the host responses of citrus hybrids developed for HLB tolerance compared to those of susceptible cultivars, assisting in developing HLB-tolerant citrus varieties. We conducted this study to validate the HLB tolerance of five Marisol Clementine × C. australis (MRL) hybrids. The CaLas titer of the five hybrids was initially assessed, followed by physiological and biochemical analyses to categorize the hybrids as susceptible or tolerant to HLB. The susceptible and tolerant hybrids were then subjected to transcriptomic profiling. Specifically, the transcriptomes of the C. australis hybrids MRL2-11 and MRL2-12 were analyzed to understand the genetic mechanism involved in HLB tolerance. The results revealed that key genes associated with various biological processes were activated in tolerant hybrids to suppress CaLas infection and reduce disease severity. Comparative genomic analysis identified polymorphic variants, confirmed parentage, and showcased the genetic diversity among the individuals. This study elucidates the essential roles of specific genes and pathways in HLB tolerance, providing insights that can be leveraged in breeding strategies and biotechnological approaches, including CRISPR/Cas gene editing. Methods Plant material Five MRL hybrids were used for this study, as outlined in Supplementary Table 1. F[1] hybrid trees were produced by crossing the Marisol clementine with the C. australis cv. DPI-205-2. Trees budded onto Kuharske rootstock were planted in the field under endemic HLB conditions, and their response to HLB was monitored for 10 years. Leaf samples from MRL hybrids were collected for periodic assessment of CaLas titer in the spring and fall. Additional leaf samples were collected for biochemical analysis and transcriptomic profiling in March (spring) when infection peaked in the field. CaLas diagnosis by qPCR Leaf samples were collected biannually (spring and fall) from the five field grown hybrid for CaLas assessment. Genomic DNA was isolated from the petioles of hybrid leaves via a GeneJET Plant Genomic DNA Purification Kit (Thermo Fischer Scientific, Waltham, MA, USA). The genomic DNA concentration was normalized to 20 ng/µL, and qPCR was performed on a QuantStudio™ 3 Real-Time qPCR machine (Thermo Fisher Scientific) with the TaqMan™ Master Mix and CQUL primers (Supplementary Table 2). CaLas genomic DNA was amplified by using primers for the CaLas rpIJ/rpIL ribosomal protein-encoding gene [[43]10]. HLB detection results are displayed as Ct values. Physiological and biochemical analysis Leaf samples were collected to analyze the chlorophyll, starch, phenol, and Malondialdehyde (MDA) contents of the MRL hybrids. Chlorophyll content estimates were generated by calculating the chlorophyll-a, chlorophyll-b, total chlorophyll, and carotenoid contents following the protocol described by Sumanta et al. [[44]11], with slight modifications. Three technical replicates were used for the study. Leaf tissue (100 mg) was homogenized in a tissue lyser II (MINI-BEADBEATER 96 by BIOSPECPRODUCTS) for 2 min at 30 Hz. The homogenized tissue was extracted by adding 1 mL of diethyl ether (DEE) and centrifuged at 10,000×g at 4 °C for 15 min. The supernatant was collected in a fresh Eppendorf tube. The absorbance was measured by adding 2.5 mL of DEE and 250 µL of extracted supernatant using a UV spectrometer at 660.6 nm and 642.2 nm, respectively. The chlorophyll-a, chlorophyll-b, total chlorophyll, and carotenoid contents were calculated using the formulas provided by Sumanta et al. [[45]11]. The starch content of the C. australis hybrids was determined using the protocol described by Rosales and Burns [[46]12], with slight modifications. Leaf tissue (100 mg) was homogenized with 700 µL of distilled water. The samples were vortexed and incubated at the boiling point for 10 min. Then, the tubes were cooled on ice and centrifuged at 6000 rpm for 2 min. A 300 µL aliquot of the supernatant was transferred to a new Eppendorf tube, and 900 µL of 100% ethanol was added before vortexing. The pellet was dissolved in 1 ml of distilled water and added 50 µL of KI: I[2] (8 mM:50 mM). The absorbance was recorded at 594 nm using a spectrophotometer. The starch content was calculated using rice starch as the standard. The total phenolic content was determined by following the method described by Singleton and Rossi [[47]13]. Leaf tissue (100 mg) was extracted in 1 ml of 80% ethanol and centrifuged at 10,000 rpm for 20 min. The supernatant was transferred to a clean glass bottle and evaporated overnight at 25 °C until completely dried. The dried residue was then dissolved in 5 mL of distilled water and used as the plant extract for analysis. For the assay, 200 µL of the extract was diluted in 3 mL of distilled water. Then, 0.3 mL of Folin reagent was added and mixed for 3 min, followed by 2 mL of a 20% Na₂CO₃ solution. The samples were incubated for 60 min in the dark, and the absorbance was measured at 650 nm using a spectrophotometer. A standard curve was prepared using ethanolic gallic acid solutions at 25 to 300 µg/mL concentrations. Lipid peroxidation was evaluated by calculating the MDA content in the leaf samples following the Heath and Packer [[48]14] method. The samples (100 mg) were homogenized and extracted in 0.5 ml of 0.1% (w/v) trichloroacetic acid (TCA). The samples were centrifuged at 14,000 rpm at 4 °C for 10 min. Furthermore, 0.5 ml of the supernatant was transferred to a new Eppendorf tube, and 0.5% (w/v) thiobarbituric acid (1 ml) (TBA) in 20% TCA was added. The samples were mixed and incubated at 95 °C in a water bath for 20 min. The tubes were immediately cooled on ice to stop the reaction, followed by centrifugation at 10,000×g at 4 °C for 5 min. The absorbance was measured at 532 nm and 600 nm. The MDA content was calculated using an extinction coefficient of 155 mM^− 1cm^− 1 and MDA levels are expressed as nmol of MDA per mg of fresh weight. Transcriptome profiling RNA extraction and transcriptome analysis Based on the CaLas titer and biochemical results, MRL2-11 and MRL2-12 were further selected for studying gene expression patterns via RNA sequencing. Leaves from MRL2-11 and MRL2-12 were processed for RNA extraction using the TRIzol method (TRI Reagent^®; Invitrogen, Waltham, MA, USA), following the manufacturer’s instructions. The quality and quantity of the extracted RNA were evaluated using a NanoDrop™ One/OneC Microvolume UV‒Vis Spectrophotometer (Thermo Scientific). RNA samples with an RNA integrity number > 6.5 were treated with oligo (dT) beads to enrich mRNAs with a poly-A tail. The mRNA was fragmented into smaller pieces, and the fragments were synthesized into first-strand cDNA using random primers. The resulting cDNA was subjected to end repair and 3’ adenylation. Adapters were ligated to the 3’ adenylated cDNA fragments and further treated with uracil-DNA-glycosylase to digest the U-labeled second-strand template, followed by PCR amplification. Library quality was assessed, and sequencing was performed on the DNBSEQ platform (BGI, Cambridge, MA, USA). Read mapping and DEG analysis The raw sequence reads in FASTQ format were processed for adapter removal and filtered to eliminate contaminants, yielding valid data using SOAPnuke software developed by BGI [[49]15]. Reads shorter than 150 bp or with N contents less than 0.1% were discarded. Additionally, low-quality reads were removed, defined as those with a quality score below 20. The clean reads were aligned to the Citrus clementina reference genome (Citrus_clementina_85681. JGI.Citrus_clementina_v1.0.v2201) using HISAT and further aligned to reference genes using Bowtie2. Differentially expressed genes (DEGs) were identified using DESeq2 by comparing MRL2-11 as a control and MRL2-12 as a treated sample. Counts were normalized, and DEGs with a p-value ≥ 0.05 were filtered out. The DEGs with log[2](fold-change) ≥ | 1 | were selected for further analysis. Significant DEGs were evaluated via MapMan analysis [[50]16], Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using Enrichment, and Gene Ontology (GO) enrichment analysis using KEGGNet. The DEGs were further categorized based on their biological functions. DEG validation via qPCR To validate the DEGs obtained from the RNA-seq data, 12 genes were selected for qPCR analysis. The primers used were synthesized by Integrated DNA Technologies (Coralville, IA, USA), as shown in Supplementary Table 3. First-strand cDNA synthesis was carried out at 1 µg/µL total RNA with a cDNA synthesis kit from Takara Bio (Shiga, Japan). The qPCR mixture included 9 µL of SYBR^® Green PowerUp™ PCR Master Mix (Applied Biosystems, Foster City, CA, USA), 1 µL of DNA (25 ng/µL), and gene-specific primers (0.25 µM), resulting in a total volume of 10 µL, following the manufacturer’s guidelines. qPCR was conducted on a QuantStudio 3 Real-Time PCR System (Applied Biosystems) in a 96-well PCR plate. Citrus β-actin was the reference gene [[51]17], with each sample tested in triplicate. Comparative genomic analysis of MRL2-11 and MRL2-12 hybrids Leaf samples were collected from the MRL2-11 and MRL2-12 trees for total DNA extraction. The DNA extraction was carried out following the protocol described by Inglis et al. [[52]18]. The quality and quantity of the DNA were assessed using agarose gel electrophoresis and a UV-Vis spectrophotometer. DNA sequencing was performed on the Illumina NovaSeq X platform, generating 250 bp paired-end reads at the University of Florida’s ICBR, Gainesville. Raw sequencing reads underwent quality control, where low-quality bases and adapter sequences were removed using Trimmomatic (v. 0.39). The sequencing quality was further evaluated with FastQC (v. 0.11.7). Genome alignment was conducted using Bowtie2, with C. australis as the reference genome. The aligned reads were further processed with Samtools to convert SAM files to BAM format, followed by sorting, indexing, and removing duplicate reads. Variant calling for MRL hybrids was performed using Freebayes (V 1.3.2) with default parameters to identify single-nucleotide polymorphisms (SNPs) and insertions and deletions (INDELs). Variant filtering was performed using VCFtools with the following parameters: VCFtools –minQ 50 –maf 0.1 –max-missing 0.8 –minDP 10 –maxDP 30 –mac 2 –recode. The filtered SNP dataset was subsequently utilized for heterozygosity and phylogenetic analysis. The genomic distribution of polymorphic variants, including SNPs and INDELs, was analyzed and annotated using SnpEff (V 5.1d). Based on annotation results, SNPs were categorized into intergenic and genic regions. Candidate genes containing polymorphic SNPs were analyzed for GO and KEGG using the Citrus Pan-Genome Breeding Database [[53]19]. The phylogenetic analysis assessed the genetic diversity between the MRL hybrids and their parental backgrounds. FASTA sequences for the parental background, including Marisol Clementine lime and C. australis, were obtained from our ongoing, unpublished project. Variant calling and filtering for MRL hybrids and their parental backgrounds were conducted using previously established parameters. The resulting filtered SNPs were used to create a dendrogram using the neighbor-joining method in Tassel V 5.2.93 [[54]20]. Homozygosity by Descent (HBD) proportions were estimated by calculating runs of homozygosity (ROH). The primarily filtered VCF file was converted to PLINK-compatible binary format. ROH analysis was performed using the -homozyg parameter to identify regions exhibiting homozygosity. The resulting BED file was processed to calculate runs of homozygosity and determine the proportion of segments having homozygosity. The ROH data were analyzed in Python using pandas and matplotlib libraries. The HBD segments were normalized based on the genome size (350 Mb) and categorized into individual classes, with no classes beyond R = 14. The lowest value (R = 2) indicated a recent ancestor, while the R > 8192 suggested a distant ancestor. The r classes were further visualized using stacked bar plots to illustrate the variation in HBD proportions among the parents and the MRL2-11 and MRL2-12 hybrids. Screening of the purity of MRL2-11 and MRL2-12 hybrids using simple sequence repeat (SSR) molecular markers Genomic DNA of Marisol clementine, C. australis, and the hybrids MRL2−11 and MRL2−12 were normalized to 25 ng/µl and used for the hybrid purity analysis. SSR primers synthesized by Integrated DNA Technologies (IDT) were utilized for hybrid testing. The SSR primers were tested on the parental genotypes to test for parental polymorphism. Polymorphic markers were evaluated in the parental genotypes along with their respective hybrids to validate hybrid purity. The polymorphic primer sequences are listed in Supplementary Table 4. The PCR conditions included an initial denaturation step (95 °C for 5 min) followed by 30 cycles of denaturation (95 °C for 45 s), annealing (50–60 °C for 45 s; depending on the temperature of the primers), and extension (72 °C for 45 s) and a final extension step at 72 °C for 7 min. The PCR mixture was prepared for a total volume of 10 µl per reaction. The reaction volume comprised 1 µl of genomic DNA normalized to a 25 ng/µl concentration, 0.75 µl of each forward and reverse primer (each at a 10 µM concentration), 4.2 µl of master mix (Promega GoTaq™ G2 Green Master Mix 2X) and 3.5 µl of nuclease-free water to make the final volume of the reaction mixture 10 µl. PCR was performed on an MJ Research (PTC−200) Peltier Thermal Cycler. The amplified product was visualized with an AXYGEN gel documentation system after resolving on a 3.5% agarose gel via agarose gel electrophoresis. Microscopic analysis of Callose deposition and phloem and xylem examination in infected tissues from MRL2-11 and MRL2-12 To evaluate and compare callose deposition in the hybrids, infected tissues from the MRL2-11 and MRL2-12 hybrids were assessed using the protocol described by Ferrara et al. [[55]21]. For callose deposition analysis, stem samples were collected from infected MRL2-11 and MRL2-12 plants by peeling the bark after making a thin cross-section to expose the vascular tissues. Immediately after collection, the tissues were fixed in 85% ethanol and incubated overnight for staining. The samples were rehydrated for 1 h in 0.01% Tween-20 solution and then transferred to 0.01% aniline blue staining solution for 1 h. After staining, images were captured on a Leica SP8 laser-scanning confocal microscope with imaging parameters based on those described in previously reported methods [[56]22]. The images were further processed for noise reduction in Fiji, and a pixel and object Classifier trained with Ilastik output was used to quantify the callose deposits following the approach described by Welker et al. [[57]23]. Four technical replicates were used for this study. Petioles were excised from MRL2-11 and MRL2-12 leaves and placed into glass vials containing 3% (v/v) glutaraldehyde in 0.1 M potassium phosphate buffer (pH 7.2). The vials were refrigerated at 3 °C overnight for primary fixation, then washed three times with phosphate buffer. Samples were fixed in 1% osmium tetroxide at room temperature for 4 h, washed again, and dehydrated in an acetone gradient (10–100% in 10% increments every 10 min), followed by two washes with 100% acetone. They were incubated overnight in 50% Spurr’s resin-acetone solution, then transferred to 70% Spurr’s resin-acetone for 8 h. Finally, samples were placed in 100% Spurr’s resin, embedded in plastic molds with fresh resin, and baked at 70 °C overnight. The resin blocks were removed and cooled to room temperature. Using a Reichert-Jung Ultracut E microtome, semithin Sect. (0.99 μm) were cut, collected with a loop, and placed on water droplets on a glass slide. A small amount of toluene was added to flatten the sections, and the slide was warmed for an hour at 65 °C. Sections were stained with methylene blue/azure A for 30 s, rinsed, counterstained with basic fuchsin for 30 s, and dried at 65 °C. They were mounted with EUKITT Neo Special and a glass coverslip. Images were captured using an Olympus BX61 microscope (10x resolution) with an Omax A35140U camera. Phloem and xylem areas were measured using ImageJ [[58]24]. Statistical analysis The biochemical data was analyzed using one-way ANOVA, with four replicates per plant. The post hoc Tukey–Kramer HSD test was performed to determine whether comparisons among plants reached a significance level of P < 0.05. Differential gene expression was assessed using the threshold cycle (Ct) method. Gene expression was normalized to the Ct value of the reference gene, producing a ∆Ct value for each sample replicate. Relative gene expression was then calculated using the 2^−ΔΔCt method [[59]25]. Statistical comparisons of gene expression profiles based on qPCR data were performed via Student’s t-test with a significance level of P < 0.05 in SigmaPlot 11.0. Results HLB assessment in the C. australis (MRL) hybrids Periodic assessments (spring and fall) were performed on the MRL hybrids via qPCR analysis. After the two-year assessment, the MRL2-12 hybrid consistently exhibited higher Ct values above 33, which were greater than those of the other MRL hybrids (Fig. [60]1). The CaLas titer in the leaf petioles of MRL2-12 was borderline positive, with Ct values of 33 and 37 in spring and fall, whereas the other MRL hybrids, MRL2-3, MRL2-11 and MRL2-14, exhibited Ct values ranging from 26 to 28, indicating that they are HLB susceptible. MRL2-1 exhibited Ct values of 33 and 30 during these periodic assessments and was considered moderately tolerant. Overall, based on the CaLas titer, compared with other MRL hybrids, MRL2-12 was categorized as tolerant, and MRL2-1 was categorized as moderately tolerant to HLB. Therefore, we further investigated the biochemical differences among the MRL hybrids to determine how HLB affects biochemical processes. Fig. 1. [61]Fig. 1 [62]Open in a new tab CaLas titer in five Marisol clementine × C. australis (MRL) hybrids in spring and fall according to qPCR analysis. The level of significance was analyzed using a paired t test, with the P value indicating a significant difference between the two phenotypes (*** P < 0.001; NS, not significant) Physiological and biochemical analysis Biochemical analysis revealed significant differences among the hybrids. The highest chlorophyll content was detected in MRL2-12 compared with the other four hybrids. The MRL2-12 hybrid exhibited the highest chlorophyll-a (15.27 µg/mL), chlorophyll-b (4.93 µg/mL), and carotenoid (5.29 µg/mL) contents (Fig. [63]2A-C). In contrast, the other C. australis hybrids exhibited similar ranges of chlorophyll contents, with the lowest chlorophyll-a content observed in the MRL2-1 hybrid. Compared with all the other C. australis hybrids, MRL2-12 exhibited significant differences in total chlorophyll content (Fig. [64]2D). Fig. 2. [65]Fig. 2 [66]Open in a new tab Biochemical analysis of the leaves of five Marisol clementine × C. australis (MRL) hybrids. A- Chlorophyll a (Chl-a), B- Chlorophyll b (Chl-b), C- Total chlorophyll, D- Carotenoid content, and E- Starch content. The level of significance was analyzed by Tukey’s–Kramer honestly significant difference test, and different letters indicate significant differences at the level of 0.05. The error bar represents the standard deviation (SD) All the hybrids exhibited similar starch accumulation except for MRL2-12 (Fig. [67]2E). MRL2-11 had the highest starch accumulation (2.58 µg/mL per mg), whereas MRL2-12 had the lowest starch accumulation (0.88 µg/mL per mg). The total phenolic and MDA contents were also estimated but did not significantly differ among the five hybrids (Supplementary Fig. 1). Statistical analysis revealed the observed differences to be nonsignificant. Biochemical analysis revealed significant differences among the MRL hybrids. The MRL hybrids were further categorized based on the CaLas titer and biochemical responses: MRL2-12 was tolerant, whereas MRL2-1, MRL2-3, MRL2-11, and MRL2-14 were susceptible. MRL2-11 was the most vulnerable among the susceptible hybrids, as indicated by its low Ct value and high starch accumulation. In contrast, MRL2-12 exhibited markers of HLB tolerance, including increased chlorophyll content and reduced starch accumulation. Comparative transcriptomic analysis of MRL hybrids Given their contrasting biochemical profiles, we selected MRL2-11 and MRL2-12 for transcriptomic profiling to investigate the molecular mechanism underlying HLB tolerance and susceptibility. In the comparative analysis, MRL2-11 served as the control, while MRL2-12 was treated as the experimental condition for differential gene expression analysis. The transcriptome was examined to elucidate the molecular dynamics and biological factors underlying the increased tolerance of MRL2-12 to HLB compared with that of the MRL2-11 hybrid. RNA sequencing was conducted using the DNBSEQ platform, resulting in approximately 44 million paired-end reads, each with a length of 150 bases and a quality threshold of Q30 > 90 (Supplementary Table 5). The average raw reads were approximately 44.38 million for MRL2-11 and 43.62 million for MRL2-12. These reads were further subjected to cleaning and filtering processes. The cleaned reads were aligned to the C. clementina genome, achieving average mapping ratios of 92.2% for MRL2-11 and 94.16% for MRL2-12. Additionally, alignment with reference genes revealed average total and uniquely mapped gene percentages of 72.11% and 67.37% for MRL2-11 and 73.52% and 68.45% for MRL2-12, respectively (Supplementary Table 5). The tolerant hybrid MRL2-12 was considered the test sample, whereas the susceptible hybrid MRL2-11 served as a control for differential gene expression analysis. Comparative analysis of MRL2-11 and MRL2-12 revealed 2890 differentially expressed genes (DEGs), comprising 1808 upregulated and 1082 downregulated genes defined by a log2-fold change ≥ | 1 |. A scatter plot of the data was demonstrated to visualize the proportion of upregulated and downregulated genes in MRL2-12 compared with MRL2-11 (Fig. [68]3). Fig. 3. [69]Fig. 3 [70]Open in a new tab Scatter plot and heatmap showing expression profiles of differentially expressed genes (DEGs) in Marisol clementine × C. australis (MRL) hybrids. A- Scatter plot of differentially expressed genes (DEGs) in MRL2-11 and MRL2-12. DEGs presented as red points are upregulated, DEGs presented as blue points are downregulated, and genes presented as gray points did not show expression changes. B- Heatmap showing all DEGs in MRL2-11 and MRL2-12. The columns indicate replicates of individual hybrids, and the rows represent individual DEGs. Red represents upregulated DEGs, and blue represents downregulated DEGs Functional annotation of DEGs GO enrichment analysis and KEGG analysis were performed to understand the functional categories of these DEGs. GO enrichment analysis categorized DEGs according to three categories: biological process, cellular component, and molecular function (Fig. [71]4A). The top five most significantly enriched biological process GO terms were defense response (145) and signal transduction (77), followed by response to salt stress (73), cell wall organization (64), and carbohydrate metabolic process (58). Most of the biological terms enriched in this dataset from MRL2-11 and MRL2-12 were associated with plant cellular responses, biotic stress, plant defense and hypersensitive responses to pathogen infection. The most significantly enriched molecular function GO terms were ATP binding (510), ADP binding (113), iron ion binding (83), and carbohydrate binding (72). The largest number of DEGs were enriched in cellular components, including membrane (825), followed by plasma membrane (548), cytoplasm (414), and cytosol (391). The other cellular component GO terms were also enriched in the extracellular region and cell wall. Fig. 4. [72]Fig. 4 [73]Open in a new tab Functional annotation of differentially expressed genes (DEGs) in Marisol clementine × C. australis (MRL) hybrids. A- Gene ontology enrichment in Marisol clementine × C. australis (MRL) hybrids DEGs by GO term classification. GO annotations categorized into three main categories as BP-biological process, CC-cellular component, and MF-Molecular functions. The top 10 terms of all three categories were represented with the adjusted p-value < 0.05. B- KEGG pathway enrichment of differentially expressed genes in Marisol clementine × C. australis (MRL) hybrids KEGG pathway enrichment analysis also revealed significant enrichment of DEGs related to protein processing in the endoplasmic reticulum, biosynthesis of cofactors, starch, and sucrose metabolism (Fig. [74]4B). Important KEGG pathways related to the response to pathogen infection, such as biosynthesis of various plant secondary metabolites, isoflavonoid biosynthesis, glycerophospholipid and glycerolipid metabolism, and ABC transporters, were also enriched. The highest KEGG pathway-rich ratio was 0.23, which was associated with protein processing in the endoplasmic reticulum. The upregulated and downregulated DEGs were further analyzed and categorized based on their biological functions. MapMan analysis revealed significant differences in response to biotic stress, with the most important pathogen recognition, signaling, defense genes, and key transcription factors upregulated in MRL2-12 compared with MRL2-11. A detailed categorization of the genes related to pathogen defense responses is presented below. DEGs related to the defense response The upregulated and downregulated DEGs were visualized to identify genes related to defense responses. These identified pathways provided insights into the mechanism underlying the increased HLB tolerance of the MRL2-12 hybrid compared with MRL2-11. Differences in hormone signaling pathways Most transcripts related to brassinosteroid were significantly downregulated in MRL2-12, and only three transcripts were upregulated. Among these transcripts, the Bbox zinc finger (zfB_box) is involved in the negative regulation of the brassinosteroid-mediated signaling pathway. In rice, brassinosteroid suppresses SA-mediated immunity [[75]26]. Significant downregulation of brassinosteroid-related transcripts may inhibit HLB in tolerant citrus cultivars [[76]8]. MapMan analysis of the RNA-seq data revealed significant upregulation of auxin-related transcripts, whereas ABA-related transcripts were moderately expressed in this dataset (Fig. [77]5). The auxin-related genes identified in this dataset included an auxin efflux carrier family protein, auxin response factor 3, auxin-responsive protein, auxin transporter-like protein 2, IAA acid amido synthetase GH3.5, IAA GH3, and cytochrome_B561 (Table [78]1). Auxin is a key hormone in pathogenesis and plant defense and promotes the expression of expansin in many crops, such as tomato, rice, and soybean [[79]27, [80]28, [81]29, [82]30]. This dataset revealed increased expression of expansin, as described in the above sections. Fig. 5. [83]Fig. 5 [84]Open in a new tab MapMan analysis of differentially expressed genes involved in biotic stress. This study revealed DEGs involved in biotic stress in MRL2-12 compared with MRL2-11 in response to HLB infection. Blue represents upregulated DEGs, and red represents downregulated DEGs Table 1. Differential expression of genes involved in defense and signaling networks in the MRL2-12 hybrid compared with the MRL2-11 hybrid in response to HLB infection Gene Id Description diffexp_log2fc Defense response Ciclev10010370m Legume lectin domain 10.62 Ciclev10009968m Isoaspartyl peptidase/L-asparaginase 1 9.42 Ciclev10022835m Hydrophobic seed protein 8.96 Ciclev10027594m Disease resistance protein RFL1-related 8.71 Ciclev10017274m Probable lipid transfer (LTP_2) 7.99 Ciclev10002288m Pathogenesis-related protein 5 6.11 Ciclev10033090m NPR1/NIM1-like defense protein C 5.94 Ciclev10032973m MLP-like protein 423-related 5.86 Ciclev10003446m MLO family (MLO) 5.72 Ciclev10006102m Pathogenesis-related protein Bet v I family 5.70 Ciclev10004811m Mate efflux family protein 5.28 Ciclev10024379m Disease resistance protein-related 4.76 Ciclev10007167m Dirigent protein 16-related 4.70 Ciclev10018317m Disease resistance-like protein/LRR domain-containing protein 3.64 Ciclev10006212m Potato inhibitor I family (potato inhibit) 3.64 Ciclev10020781m Kdeltailed cysteine endopeptidase CEP1-related 3.30 Ciclev10016745m Cotton fiber expressed protein (DUF761) 3.03 Ciclev10014260m Lipoxygenase 2.62 Ciclev10026624m Cap (cysteine-rich secretory proteins, antigen 5, and pathogenesis-related 1 protein) 3.61 Ciclev10032523m Cysteine-rich repeat secretory protein 55 2.70 Ciclev10009569m Cysteine-rich secretory protein-related 3.08 Ciclev10001885m Pathogenesis-related thaumatin family 2.34 Ciclev10031364m Aspartyl protease family protein 2.33 Ciclev10030550m TIR-NBS-LRR Class disease resistance protein -1.22 Ciclev10027774m Respiratory burst oxidase homolog protein -1.39 Ciclev10025008m Disease resistance protein RFL1-related -1.81 Ciclev10006106m Pathogenesis-related protein Bet v I family -1.96 Ciclev10022100m Thaumatin family -2.21 Ciclev10022863m Protein NIM-interacting 1 -2.28 Ciclev10019279m Protein eceriferum 1 -6.11 Ciclev10030379m Probable lipid transfer (LTP 2) -6.57 Signal transduction Ciclev10024007m G protein-coupled seven transmembrane receptors 1.33 Ciclev10025184m Protein hothead 7.99 Ciclev10019884m Protein IQ domain 21 3.73 Ciclev10016797m RAC-like GTP-binding protein ARAC2 2.68 Ciclev10027581m Receptor-like protein 13 3.20 Ciclev10012256m Calcium/calmodulin-dependent protein kinase -1.08 Ciclev10010314m Calcium-dependent protein kinase 29 -2.35 Ciclev10027761m Extra large guanine nucleotide-binding protein 1 -1.32 Ciclev10033583m Mitogen-activated protein kinase (MAPK) kinase MKK3/MKK6 -4.01 Ciclev10024591m Mitogen-activated protein kinase kinase 7-related -8.62 Brassinosteroid-mediated signaling pathway Ciclev10031341m Betaamyrin 11oxidase/CYP88D6 2.77 Ciclev10001255m C2H2 like zinc finger protein 1.25 Ciclev10004935m Cytokinin 7 betaglucosyltransferase/Uridine diphosphoglucosezeatin 7 glucosyltransferase -6.49 Ciclev10021866m Brassinosteroid metabolic pathway protein ben1 -1.44 Ciclev10031848m Cytosolic sulfotransferase 12 -5.89 Ciclev10011890m Protein kinase family protein -1.09 Response to gibberellin Ciclev10017244m Extensin, proline-rich protein 4.50 Ciclev10024024m Fbox protein GID2 (GID2, SLY1) 2.50 Ciclev10015702m Gibberellin 2betadioxygenase 7-related 1.80 Ciclev10013200m Gibberellin regulated protein 12-related 2.05 Ciclev10029695m Gibberellin-regulated protein 7 5.99 Ciclev10023706m Entcopalyl diphosphate synthase -3.51 Ciclev10017993m Gibberellin A(4) carboxyl methyltransferase -1.34 Ciclev10002984m Gibberellin-regulated protein 1 -1.45 Ciclev10004586m GRAS domain family (GRAS) -2.58 Response to jasmonic acid Ciclev10006614m Copper amine oxidase 7.30 Ciclev10011156m Cullins 3.36 Ciclev10004275m Inosine uridine preferring nucleoside hydrolase family protein 4.00 Ciclev10014209m Linoleate 13 S-lipoxygenase 5.43 Ciclev10018092m Strictosidine synthase-related 2.85 Ciclev10004511m Amine oxidase -1.14 Ciclev10007424m Glutamate receptor 2.1-related -5.37 Ciclev10013735m Ribulose bisphosphate carboxylase -3.85 Response to calcium-mediated signaling Ciclev10017476m Cyclic nucleotide-gated ion channel 1 8.05 Ciclev10027257m Cyclic nucleotide-gated ion channel 7-related 2.16 Ciclev10028938m Cyclin B23-related 3.38 Ciclev10029092m Cyclin 4.50 Ciclev10026759m Protein RALF-like 26-related 4.71 Ciclev10010358m Synaptotagmin 4 3.52 [85]Open in a new tab Transcripts involved in the response to ethylene, including tetratricopeptide repeat protein (TPR), ARGOS-like protein-related, and basic endochitinase B, were also upregulated. Significant changes were observed in the profiles of jasmonic acid (JA)- and gibberellic acid (GA)-related transcripts, but interestingly, the levels of SA-related genes were not significantly affected in this dataset. The induction of genes involved in the GA-mediated signaling pathway was identified, and the downregulation of GA biosynthetic genes was observed. The expression of GRAS domain family-related genes was significantly downregulated, indicating their involvement in the negative regulation of the GA-mediated signaling pathway. The expression of transcripts involved in the response to JA was highly upregulated in this dataset. The expression of the transcript encoding linoleate 13 S lipoxygenase, a key enzyme in JA synthesis, was significantly upregulated (log2 FC = 5.4) in MRL2-12. However, more transcripts related to glutamate receptors 2.1, 2.5, and 3.1 were significantly downregulated in the MRL2-12 hybrid (Table [86]1). This dataset indicates significant changes in the expression profiles of transcripts in the auxin, GA, JA, and ethylene pathways. These alterations likely play crucial roles in the CaLas‒citrus interaction by affecting plant physiology and contributing to the development of HLB symptoms. The expression of transcripts related to Ca^2+-mediated signaling, including cyclic nucleotide-gated ion channel 1 (CNGC), which had a log2 FC of 8.04, CNGC14, and CNGC7, was significantly increased. CNGCs are ion channels regulated by cytosolic signaling molecules, and insect-mediated injury leads to the activation of Ca²⁺ CaM-dependent phosphorylation. This triggers a rapid JA burst and activates plant defenses against pathogens [[87]31]. Ca²⁺ also acts as a secondary messenger in signal transduction pathways. After the initiation of CaLas infection, an increase in cytosolic Ca²⁺ levels can trigger a cascade of signaling events that activate plant defense mechanisms. Differences in signaling factor and signaling receptor expression between the hybrids Plant signaling receptors play crucial roles in the recognition of external stimuli, such as pathogen-associated molecular patterns (PAMPs) and elicitors, leading to the activation of immune responses against pathogens. MapMan analysis revealed significant upregulation of DEGs involved in pathogen recognition, the respiratory burst, and defense response signaling in response to HLB infection (Fig. [88]5). This dataset revealed differential expression of many transcripts encoding receptor-like proteins (RLPs) and receptor-like kinases (RLKs), including G-type and L-type lectin domain-containing receptor-like kinases, leucine-rich repeat (LRR) kinases, serine-threonine protein kinases, and wall-associated kinases (WAKs), in response to pathogen infection. Nearly an equal number of LRR kinases were upregulated and downregulated in the MRL2-12 hybrid. Most protein kinase-related transcripts exhibited moderate expression in this dataset, with a few uniquely upregulated and downregulated transcripts. Compared with that in the MRL2-11 hybrid, the expression of the transcripts encoding cysteine-rich receptor-like protein kinase 27, RLP54, RLP46, and PR5-like receptor kinase was significantly downregulated in MRL2-12. Some of the uniquely expressed RLPs/RLKs in the MRL2-12 hybrid included RLP 13, RLP 55, receptor-like protein kinase TMK1, receptor-like protein kinase HSL1, Ras suppressor protein, L-type lectin domain-containing receptor kinase 1, archaeal ATPase, and WAK receptor 1 (Supplementary Table 6). Notably, RLP13 and WAK1 exhibited relatively high transcript abundances and expression levels in the MRL2-12 hybrid. The significantly increased expression of these receptors may contribute to the increased tolerance to CaLas of MRL2-12. Plant hormone receptors, such as phytohormone receptors, are also directly associated with signaling pathways, which trigger defense responses against pathogens [[89]32]. This receptor activity leads to significant upregulation of associated signaling cascades, typically involving MAPK5, MAPKK6, and MAPKKK-SSK2. MAPKs regulate the gene expression of antimicrobial compounds, reinforce the cell wall, and activate other defense-related processes. Our data revealed that different cyclin family-related transcripts and cyclin-dependent protein kinases were uniquely and significantly upregulated in the MRL2-12 hybrid compared with the MRL2-11 hybrid. A greater number of transcripts related to cyclins, which act as regulatory subunits that bind to and activate cyclin-dependent kinases (CDKs), were significantly upregulated in this dataset. The significantly increased expression of cyclins and CDKs could influence the plant response to disease due to the roles of these proteins in cell cycle regulation, stress responses, and interactions with hormone pathways in the context of HLB-induced stress. These findings also indicate that CDKs increase plant stress tolerance by controlling the phosphorylation of SR-splicing factors [[90]33]. Differences in defense-related transcript activity between the hybrids Defense-related genes are directly involved in the production of compounds that can inhibit a pathogen and limit its growth. Many disease tolerance genes were significantly upregulated in the MRL2-12 hybrid in response to HLB infection. Specifically, six transcripts of the disease tolerance protein RFL1 were highly expressed, with the highest log2FC value (8.7), whereas one transcript was downregulated. This dataset revealed key pathogenesis-related (PR) proteins that were significantly upregulated, such as NPR1-like C-terminal, PR 5, the PR protein Bet v I family-PR10, and the PR thaumatin family, in the MRL2-12 hybrid. Conversely, two transcripts of each PR protein, Bet v I family-PR10, the PR thaumatin family, and one transcript of a negative regulator of NPR1, i.e., NIM1-interacting protein, were also downregulated in the MRL2-12 hybrid in response to HLB infection. The transcript encoding EDS1, a transmembrane signaling receptor involved in SAR responses, was upregulated in the MRL2-12 hybrid [[91]34]. In addition, other important defense-related genes that were significantly and uniquely upregulated in the MRL2-12 hybrid included dirigent proteins 16 and 19, DUF761, AIG1 domain-containing protein, ATP binding cassette subfamily F member 3, late embryogenesis abundant hydroxyproline-rich glycoprotein, MLO11, MLO12, MLO13, MLP-like protein 423-related, multidrug resistance protein, gem-like protein 4, legume lectin domain, D-mannose binding lectin (B_lectin), AAA domain, and hydrophobic seed protein. The legume lectin domain transcripts exhibited the highest expression among the defense-related genes, with a log2FC value of 10.62, indicating their probable involvement in plant defense against insect pests and bacteria [[92]35]. Three transcripts of hydrophobic seed protein (log2 FC-8.95) encode a putative lipid transfer protein, which is also involved in the priming of the SAR and ISR responses, particularly in the spread of cell-to-cell signals during pathogen attack [[93]36]. Some defense-related genes, such as probable lipid transfer protein (LTP 2), potato inhibitor I family (potato inhibit), mate efflux family protein, and Ca^2+-independent phospholipase A2, were moderately expressed in MRL2-11 and MRL2-12. Genes encoding lipid transfer proteins play a role in inhibiting bacterial growth [[94]37], and they were differentially expressed in this dataset. The expression of these genes revealed the activation of the defense response in both hybrids in response to CaLas infection. The dataset revealed that four transcripts of cysteine-rich secretory proteins, antigen 5, and pathogenesis-related 1 protein (CAP) superfamily proteins, were upregulated in MRL2-12; these proteins are involved in immune defense and pathogen response pathways [[95]38]. The transcript encoding isoaspartyl peptidase (AP) was highly expressed, with a log2FC value of 9.4 in the MRL2-12 hybrid. Plant aspartic proteases are proteolytic enzymes crucial for plant defense mechanisms during plant‒pathogen interactions [[96]39]. The transcript encoding AP was abundant and significantly upregulated in this dataset in the MRL2-12 hybrid. The transcripts encoding cysteine proteases of the C1 family and the papain-like cysteine protease CEP1 were upregulated in the MRL2-12 hybrid compared with MRL2-11. C1 cysteine proteases can degrade pathogen proteins, inhibit pathogen growth, and limit disease progression [[97]40]. Transcripts involved in cellular responses were differentially regulated between the hybrids Many transcripts related to cell wall modification were differentially expressed in the MRL2-12 hybrid. Transcripts encoding xyloglucan endotransglucosylase (XTR) 10, XTR 25, XTR 33, and xyloglucan glycosyltransferases 4 and 5 were significantly upregulated in MRL2-12. Most cell wall modification-related DEGs, including cellulose synthases, expansins, and pectinesterases, were upregulated in MRL2-12. Pectin hydrolysis in response to bacterial infection [[98]41] was observed, and many genes related to pectin synthesis and degradation were identified in this dataset (Table [99]2). Another important transcript, plant invertase/pectin methylesterase inhibitor and pectin esterase inhibitors 12, 22, and 34, was significantly upregulated in MRL2-12. Pectin methyl esterase plays a role in pectin modification, and its inhibitor strengthens the cell wall against HLB infection, a basal cellular defense response [[100]42, [101]43]. Expansin helps modify the cell wall structure, a crucial physical barrier against pathogen invasion. Remodeling of the cell wall by expansin can increase a plant’s ability to prevent pathogen entry [[102]44]. The dataset revealed several expansin transcripts, including A13, A6, A8, A1 and B1. Overall, many upregulated DEGs related to cell wall biogenesis, primary and secondary cell wall modifications, and organization were identified in this dataset. The transcripts involved in these cellular defenses include protein trichome birefringence-like protein, pectate lyase, laccase, germin-like protein, cobra-like extracellular glycosyl-phosphatidyl inositol-anchored family protein, cinnamoyl coa reductase-like protein, galacturonosyltransferase, and fasciclin-like arabinogalactan protein. Very few genes involved in cellular defense responses were downregulated, such as transcripts encoding alpha/beta hydrolase, inorganic pyrophosphatase-like protein, and casparian strip membrane protein 5. It was also observed that genes encoding chitinase commonly exhibited altered expression, whereas two chitin elicitor receptor kinase 1 transcripts were downregulated in this dataset. Table 2. Differentially expressed genes involved in cell wall modification and stress responses in the MRL2-12 hybrid compared with the MRL2-11 hybrid in response to HLB infection Gene Id Description diffexp_log2fc Cellulose synthesis Ciclev10000288m Cellulose synthase 3.12 Ciclev10018639m Cellulose synthase A catalytic subunit 7 4.39 Ciclev10018574m Cellulose synthase-like protein D5-related 2.67 Cell wall modification Ciclev10015643m alpha/beta hydrolase fold (Abhydrolase 1) -3.67 Ciclev10033634m Casparian strip membrane protein 5 -6.46 Ciclev10002604m Cell wall/vacuolar inhibitor of fructosidase 2-related 2.38 Ciclev10015683m Cinnamoyl COA reductase-like protein 1.39 Ciclev10020292m Cobra-like protein 4 4.62 Ciclev10016660m Expansin B1-related 2.92 Ciclev10029088m Expansin-like A1-related 2.12 Ciclev10006347m Extensin proline-rich protein 6.24 Ciclev10017676m Fasciclin-like arabinogalactan protein 20-related 5.58 Ciclev10013506m GDSL esterase/lipase EXL1 8.04 Ciclev10011798m GDSL/SGNHlike AcylEsterase family found in Pmr5 and Cas1p 3.75 Ciclev10011358m Laccase2 7.01 Ciclev10011635m Pectate lyase 12-related 3.69 Ciclev10007806m Pectin methylesterase 2-related 6.09 Ciclev10025336m Plant protein of unknown function (DUF247) 4.43 Ciclev10005627m Pollen allergen (Pollen_allerg_1) 4.54 Ciclev10021897m Pollen allergen (Pollen_allerg_1) -5.36 Ciclev10031403m Protein Eskimo 1 3.32 Ciclev10011779m Protein trichome birefringence-like 3 6.06 Ciclev10003963m Transducin/wd40 domain-containing protein -6.42 Ciclev10002730m Transmembrane protein 18 (TMEM18) -1.08 Ciclev10023335m UDP-N-acetylglucosamine 1carboxyvinyltransferase 5.52 Ciclev10032164m Xyloglucan endotransglucosylase/hydrolase protein 33 5.31 Ciclev10014531m Xyloglucan glycosyltransferase 4 2.42 Phloem Ciclev10003574m Fbox domain/Phloem protein 2 (PP2) 5.63 Ciclev10033113m Phloem protein 2 (PP2) 2.69 Ciclev10025490m UDP glycosyl transferase 72E1 2.38 Ciclev10014352m Anthranilate phosphoribosyl transferase-like 2.70 Ciclev10023804m Delta(4)-3-Ketosteroid 5β-Reductase 4.17 Ciclev10005557m Fbox protein PP2 A14 2.04 Ciclev10002074m Fbox protein PP2 B13 1.62 Ciclev10019128m Sieve element occlusion N-terminus (SEO_N) 2.70 Ciclev10025332m Plant protein of unknown function (DUF936) 3.68 Ciclev10025189m Protein NRT1/PTR family 1.1-related 1.68 Ciclev10021426m PP2 A13 1.64 Ciclev10015879m Dof domain, zinc finger (zfDof) -1.07 Ciclev10002204m Fbox protein PP2 B13 -6.74 Ciclev10002065m Fbox protein PP2 B1 -4.46 Ciclev10001839m Homogentisate phytyltransferase/HPT -1.56 Ciclev10031318m Protein DA1-related 2 -1.11 Stress response Ciclev10003792m Glutathione S transferase GST 7.62 Ciclev10026130m Peroxidase/Lactoperoxidase 7.56 Ciclev10001877m Peroxidase 64 22.07 Ciclev10009079m IQ domain 9 protein 5.59 Ciclev10022871m Glutathione peroxidase 8-related -1.13 Ciclev10002678m Glutathione s transferase -4.89 Ciclev10031883m IQ calmodulin-binding motif (IQ) -2.65 Ciclev10012554m L-ascorbate peroxidase 2, cytosolic -3.29 Ciclev10021010m Peroxidase/Lactoperoxidase -3.24 Proteases and inhibitors Ciclev10022156m Trypsin and protease inhibitor 11.07 Ciclev10018236m Cysteine protease family C1-related 4.52 Ciclev10003531m Cysteine proteinase-like protein-related 1.32 [103]Open in a new tab Transcripts encoding the phloem proteins PP2-B13 and PP2-B1 were strongly downregulated, whereas those encoding PP2-A14, PP2-A13, and PP2-A2 were significantly upregulated in response to HLB infection. PP2-B13 was moderately regulated in this dataset. In addition to these genes, transcripts related to phloem formation and transport, such as DUF936, sieve element occlusion N-terminus (SEO N), protein NRT1/PTR family 1.1 and 5.8, and testosterone 5-beta reductase, were also upregulated in this dataset. In contrast, protein DA1-related 2, which plays a role in phloem development in response to HLB infection, was downregulated. The gene encoding transmembrane protein 18, which plays a role in callose deposition in the cell wall, was downregulated in MRL2-12. This dataset revealed altered expression of phloem genes and the absence of callose-induced transcript expression. This result suggested that MRL2-12 modulated the expression of phloem genes in response to HLB infection without overdeposition of callose, indicating that it might not cause any phloem function disorders. Differences in biotic stress responses between the hybrids Transcripts encoding peroxidases such as peroxidases 64, 55, 48, 42, 22, and 21 were significantly upregulated in the MRL2-12 hybrid in response to HLB infection (Fig. [104]5). Among these genes, peroxidase 64 exhibited the highest expression, with a log2FC value of 22.07. The accumulation of peroxidases in citrus indicates their involvement in plant defense activation [[105]45]. In citrus, redox-related genes have been linked to HLB tolerance in both young and mature leaves [[106]46]. Genes encoding glutathione S-transferase were consistently modulated in the MRL2-12 hybrid. Other stress-related genes, including L-ascorbate peroxidase 2, superoxide dismutase, lactoperoxidase, calmodulin-binding family protein, glutathione peroxidase 8, and PCD4, were downregulated in this dataset. Conversely, major stress-related genes, including aquaporin PIP14, CML30, calmodulin-like protein 3, LACT3, peroxisomal sarcosine oxidase, glutaredoxin, and glutathione hydrolase, were upregulated in the MRL2-12 hybrid compared with the MRL2-11 hybrid (Table [107]2). Transcripts encoding heat shock proteins were differentially expressed in this dataset in response to HLB infection. Six transcripts encoding trypsin and protease inhibitors (Kunitz legume), responsive to peroxidase, were significantly upregulated in MRL2-12, with the highest fold change in expression being 11.07. Transcription factors The plant defense response requires large-scale transcriptional reprogramming in response to pathogen attack. This explains why this dataset’s transcriptional expression of transcription factors, including WRKY, MYB, AP2 domain, BHLH, ERF, and NAC domain, was differentially modulated in response to pathogen infection. Most of the 66 transcription factors and related domains were upregulated, and 39 were downregulated in MRL2-12. The data revealed many transcripts related to ethylene-responsive transcription factors, general transcription factor 2-related zinc finger proteins, GATA transcription factors, SHN shines, DNA binding/transcription factors, the transcription factor MYC1, the transcription factor XBOX, and the MYB-like DNA binding protein MYB, which were uniquely and strongly upregulated in MRL2-12. The data revealed only two upregulated transcripts of WRKY45, whereas more transcripts related to the WRKY domain were downregulated in the MRL2-12 hybrid. The transcription factor MYB plays crucial roles in disease tolerance and other biological processes, and this dataset identified a greater number of transcripts related to the MYB domain. More transcripts associated with MYB and NAC domain transcription factors were significantly upregulated, as shown in Supplementary Table 7. Most transcription factor domains were upregulated in both hybrids, whereas some were uniquely upregulated in MRL2-12, potentially regulating genes responsive to HLB infection. The uniquely upregulated transcription factor domains in the MRL2-12 hybrid included general transcription factor 2, CDP1, DCP1, TT2, MYC1, SHN1, and ADA2. Additionally, ethylene-responsive factors (ERFs) involved in ethylene signaling, such as ERF13, ERF16, ERF35, and ERF SHINE2 transcription factor domains, were upregulated, whereas the CRF5 and ERF096 transcripts related to ERF transcription factors were downregulated. Differences in secondary metabolites between the hybrids Several transcripts associated with biological processes related to secondary metabolism were identified in the MRL2-12 hybrid. These DEGs are involved in the biosynthesis of metabolites, including terpenoids, oxylipins, isoprenoids, flavonoids, and phenylpropanoids. This dataset revealed moderate alterations in flavonoid 3’ monooxygenase and O-methyl transferase (OMT) expression in the MRL2-12 hybrid. Many OMT transcripts involved in flavonoid and phenylpropanoid metabolism were considerably upregulated, with the highest log2 FC value of 9.83 among all OMT transcripts. Phenylpropanoid and flavonoid biosynthesis were more highly enriched in the MRL2-12 hybrid than in the MRL2-11 hybrid. Transcripts involved in terpenoid and isoprenoid biosynthetic processes, including nerolidol synthase, beta-amyrin synthase, geraniol 8-hydroxylase, all-trans-nonaprenyl diphosphate synthase, geranylgeranyl diphosphate synthase, and dehydrogenase, were uniquely upregulated in the MRL2-12 hybrid. An increase in secondary metabolite levels can be triggered following CaLas infection, which might contribute to HLB tolerance in the MRL2-12 hybrid. Other genes involved in secondary metabolism-related processes, including serine carboxypeptidase-like 1, shikimate o-hydroxycinnamoyl transferase, phenylalanine N-monooxygenase, flavonol synthase, acyl-coa desaturase, and norcoclaurine synthase, were significantly upregulated in this dataset. In contrast, the downregulated genes were mainly associated with UDP-glucosyl transferase 73B2 and shikimate kinase. The transcripts encoding glucosyl transferase and cytochrome P450 were moderately expressed in the MRL2-12 hybrid. These genes are involved in plant secondary metabolism, and responses to biotic and abiotic stress might play a role in defense responses. Differences in carbohydrate metabolism, transporters, and miscellaneous functions between the hybrids Biological functions such as carbohydrate metabolism and the transport of substances were strongly affected by HLB infection. Genes related to these biological functions were differentially expressed in the MRL2-12 hybrid. Some major transcripts encoded alpha/betahydrolases, o-glycosyl hydrolase family 17, beta-galactosidases 2, 3, 10, and 45, malate dehydrogenase, and pectate lyase 5. A few transcripts, including beta-galactosidase 1 and NAD-dependent epimerase/dehydratase, were modulated in both hybrids. The primary downregulated transcripts in the MRL2-12 hybrid included starch synthase 1, AAGR1, alpha-amylase 2, and DJ1. Genes involved in starch synthesis and degradation were downregulated in the MRL2-12 hybrid. Transcripts related to transport functions, including ABC transporter c, amino acid transporter, oligopeptide transporter, potassium transporter, and sugar transport protein, were differentially expressed in the MRL2-12 hybrid. The transporter family genes that were most affected were related to the transport of ions, amino acids, and sugars. Some of the genes uniquely upregulated in this dataset included a zinc/iron transporter, boron transporter 4, auxin transporter protein 1, and phosphate transporter PHO1. Ubiquitin-mediated protein degradation plays a crucial role in plant-pathogen interactions [[108]47]. This dataset identified fifteen transcripts involved in the ubiquitin-proteasome system that were upregulated, including those encoding E3 ubiquitin ligases, the E3 ubiquitin-protein ligase ATL41, XBAT31, the ubiquitin-conjugating enzyme E2, and ubiquitinyl hydrolase 1. One transcript of the ubiquitin-conjugating enzyme E2 was downregulated in the MRL2-12 hybrid. Validation of DEGs via qPCR The RNA-seq data were validated via qPCR analysis. Twelve DEGs involved in defense in response to HLB were selected for validation. Among them, six upregulated genes in MRL2-12, including peroxidase 64, PP2-A14, PR10, the BTB/POZ domain, pectin esterase inhibitor (PI) 22 and Kunitz trypsin inhibitor (KTI), were confirmed by qPCR. Additionally, six downregulated DEGs, PP2-B13, PP2-B1, starch synthase, NIMIN1, chitinase, and DUF3420, were also validated. The qPCR results of these 12 DEGs were consistent with the RNA-seq data (Fig. [109]6), confirming the reliability of the gene expression patterns identified in the transcriptomic data. Fig. 6. [110]Fig. 6 [111]Open in a new tab Differentially expressed gene validation by qPCR. The relative levels of mRNA transcript accumulation for six upregulated (A) and six downregulated (B) transcripts in comparisons of the Marisol clementine × C. australis (MRL) hybrids MRL2−11 and MRL2−12, as determined by qPCR (2^−ΔΔCt). The level of significance of differences between samples is indicated as P ≤ 0.05; this notation indicates that the calculated p value was less than or equal to 0.05 Comparative genomic analysis of MRL2-11 and MRL2-12 hybrids The whole genome sequencing of MRL2-11 and MRL2-12 hybrids produced an average of 52,990,939 and 41,363,982 reads, respectively, with a PHRED quality score exceeding 37. After removing low-quality sequences, clean reads were mapped to the reference genome and used for variant calling. Variant calling identified 19,746,782 variants, including 18,608,524 SNPs and 1,138,258 INDELs. This extensive variation highlights the genetic diversity present in MRL2-111 and MRL2-12 hybrids. After stringent filtering, 16,310 SNPs and 3,060 INDELs remained, with an average depth of 23.48 and 18.45, respectively. After filtering out the 16,310 SNPs, 5086 heterozygous SNPs were detected, resulting in a heterozygosity rate of 0.1559. (Supplementary Table 8). Annotation and effect of SNPs/INDELs The filtered SNPs and INDELs were further analyzed to identify polymorphic and unique markers. Among the filtered variants, 13,767 SNPs and 2,686 INDELs were identified as unique, polymorphic, or informative, displaying a distinct alternate allele between the MRL2-11 and MRL2-12. Unique heterozygous SNPs were annotated using the C. australis database with SnpEff to understand their plausible effects. The polymorphic variants consisted of 6,536 SNPs and 1,599 INDELs located in transcript and protein-coding regions (Fig. [112]7A). The distribution of these variants and their effect on genomic regions are also illustrated in Fig. [113]7B. The highest distribution of variants was observed in the intergenic region, followed by the upstream gene variant category in SNPs and INDELs, which affects gene regulation. Out of 6,536 SNPs found in the transcript region, 227 SNPs were identified as non-synonymous variants. The most significant INDELs were classified as frameshift mutations, stop-gained variants, and splicing-related changes, with 34 such INDELs observed. Fig. 7. [114]Fig. 7 [115]Open in a new tab SNP classification and its functional annotation. A- The number of polymorphic variants distributed in genomic and transcript regions is presented in a bar chart. B- The number of SNPs effect in the transcript region in each class is presented in a bar chart C- The number of INDELs effect in the transcript region in each class is presented in a bar chart D- GO analysis of SNP-containing genes in transcript region: The analysis is divided into three categories -biological process, cellular component, and molecular function. E- KEGG pathway enrichment: KEGG pathway analysis of SNP-containing genes distribution from transcript region is presented in a pie chart Gene ontology pathway analysis was performed on genes affected by heterozygous SNPs in the transcript region using the Citrus Pan-Genome to Breeding Database (CPBD) to assess their functional annotation (p < 0.05; Fig. [116]7C). GO enrichment analysis revealed SNP-containing genes were primarily enriched in molecular function, with “protein binding” showing the highest enrichment (146), followed by “ATP binding” and “protein kinase activity.” In the biological process, the most enriched terms were “protein phosphorylation” (59), followed by transmembrane transport (41). In the cellular component category, the most enriched terms were “membrane” (57), and “integral component of membrane” (40). Several key biological processes related to plant defense were significantly enriched, including “defense responses”, “ethylene-activated signaling pathway” and “transmembrane transport”. KEGG pathway analysis further revealed significant enrichment of candidate genes containing heterozygous SNPs in pathways associated with metabolic processes, biosynthesis of secondary metabolites, plant hormone signal transduction, and plant-pathogen interactions (Fig. [117]7D). These pathways play a crucial role in plant defense mechanisms against pathogen infections. Phylogenetic analysis The phylogenetic tree was analyzed to examine the relationship among the MRL hybrids and their parental lines, Clementine and C. australis (Supplementary Fig. 2). As expected, the MRL2-11 and MRL2-12 hybrids formed a branch in the phylogeny, originating from the Clementine and C. australis parental lines. The tree shows that MRL2-11 and MRL2-12 are closely related, while C. australis is slightly more diverged. Clementine is the most distinct among the samples, showing a significant genetic distance from the others. Homozygosity by descent (HBD) Homozygosity by descent (HBD) analysis was performed to investigate regions inherited from the common parent. The HBD analysis revealed the distribution of ROH (Runs of Homozygosity) size classes among Clementine, C. australis, and their hybrids MRL2-11 and MRL2-12 (Fig. [118]8). The hybrids, MRL2-11 and MRL2-12, share significant homozygous regions, inherited from their parents, with variations in the proportion of different ROH classes. Clementine exhibits a high proportion of long ROH segments, indicating a more inbred genetic background, whereas C. australis shows a more diverse ROH distribution. The ROH classification suggests that MRL2-12 has a higher proportion of shared homozygous segments with Clementine. Since this is the F[1] generation, ROH from the susceptible parent in the tolerant hybrid could be explained by segregation distortion [[119]48]. Fig. 8. [120]Fig. 8 [121]Open in a new tab Variation in runs of homozygosity (ROH) among parental lines and their hybrids. ROH classes represent segments inherited from the parental lines to their hybrids, ranging from the ancient class (R = 8192) to the recent class (R = 2). Microscopic analysis of callose deposition and phloem and xylem examination in infected tissues from MRL2-11 and MRL2-12 Microscopic examination of callose deposition was performed for the MRL2-11 and MRL2-12 hybrids. Processed Ilastik output images were used for callose quantification. However, image analysis revealed a higher accumulation of callose in the MRL2-11 hybrid than in the MRL2-12 hybrid (Fig. [122]9A and B). Callose deposition was observed in both MRL2-11 and MRL2-12 hybrids, with visibly higher accumulation in the MRL2-11 hybrid compared to MRL2-12 (Fig. [123]9C). However, statistical analysis indicated that the difference was non-significant (P < 0.05). Fig. 9. [124]Fig. 9 [125]Open in a new tab Callose deposition of the infected stem phloem and morphological analysis of infected petioles of Marisol clementine × C. australis (MRL) hybrids. A- MRL2-11 B- MRL2-12 hybrids. These are the unaltered images of a high-callose section of MRL hybrids stem phloem, stained with aniline blue and fully segmented image, with all objects classified correctly as callose, background and tissue artifact. C - The mean callose formation count for MRL2-11 and MRL2-12. Brightfield images of petioles of MRL2-11 (D), MRL2-12 (E). Phloem/xylem Ratio of MRL hybrids (F). The level of significance of variation between samples is indicated as P ≤ 0.05. This notation indicates that the calculated p-value is ≤ 0.05; NS- not significant Morphological differences between phloem and xylem were observed in MRL2-11 and MRL2-12 hybrids. A notable difference was observed in the phloem and xylem thickness in the tolerant hybrid MRL2-12 compared to the MRL2-11 hybrid (Fig. [126]9D and E). The mean phloem/ xylem ratio (Pa/Xa) was significantly different between MRL hybrids, with a higher ratio in the tolerant hybrid (Fig. [127]9F). This suggests that MRL2-12 regenerates phloem more effectively under HLB pressure. In contrast, the susceptible hybrid exhibited a lower ratio, indicating more severe phloem degradation without effective regeneration. Hybrid purity testing of MRL2-11 and MRL2-12 SSR markers were tested on the parental lime genotypes (Marisol clementine and Australian Round Lime) and their hybrids (MRL2-11 and MRL2-12) to determine hybrid purity. Parental genotype screening identified a total of 13 polymorphic primers. The polymorphic primers were further validated in the hybrids. SSR primers exhibiting male- and female-specific amplicons in the hybrids were used to validate hybrid purity. A total of 65% of amplicons were polymorphic between the parental genotypes. Five SSR markers that validated the amplicons for both parents were used to determine hybrid purity. The SSR markers used to validate hybrid purity are illustrated in Supplementary Fig. 3. Most SSR markers had amplicon sizes ranging from 100 to 250 bp. Discussion C. australis is a citrus species from Australia that is tolerant to HLB, and its genetics can be used to improve disease tolerance in other citrus cultivars through hybridization [[128]49]. Identifying and characterizing the factors responsible for HLB tolerance is crucial for developing resistant citrus varieties. Effective screening of tolerant hybrids primarily involves field trials, molecular analyses, and biochemical analyses, which are instrumental in advancing breeding efforts for HLB tolerance. This study aimed to identify tolerant C. australis hybrids by studying their physiological, biochemical, and molecular elements. The gene expression patterns of two hybrids with similar genetics -an HLB susceptible hybrid (MRL2-11) and a tolerant hybrid (MRL2-12) - were analyzed using RNA-seq to investigate plant defense responses in the context of HLB. Transcriptomic data revealed genes differentially expressed between the tolerant C. australis hybrid and the susceptible sibling, which were both grown under the same field conditions. This dataset revealed biological processes and pathways activated or suppressed in the HLB-tolerant C. australis hybrid MRL2-12, leading to basal tolerance against HLB. We confirmed the hybrid purity of MRL2-11 and MRL2-12 by using SSR markers to determine the authenticity of the genetic material [[129]50]. Comparative genomic analysis provided valuable insights into the polymorphic variants in MRL hybrids. ROH analysis indicated that the tolerant hybrid MRL2-12 shared a higher proportion of homozygous segments with Clementine than with the tolerant parent. However, despite MRL2-12 exhibiting tolerance, its homozygous HBD segments appear to have been primarily inherited from the susceptible parent. Previous studies have also demonstrated tolerance in Clementine-derived hybrids, such as the Sugar Belle hybrid (Clementine × Minneola) [[130]51, [131]52] or the USDA’s FF-5-51-2 (Clementine × Orlando) [[132]53], indicating that Clementine also plays a role in the HLB tolerance observed in its progeny. The phenotypic dominance of the tolerant trait may be influenced by epistatic interactions between genes from both parents while retaining homozygous regions from the susceptible parent. Silva et al. [[133]54] support this concept, aligning with observations from this study on the tolerance hybrid MRL2-12, which retained regions from the susceptible parent. Pathogenesis-related defense responses Common gene modulation patterns between the two hybrids were observed, such as altered expression of genes involved in pathogen recognition, reception, and cell wall modification, which are the initial responses generated against pathogens. For example, NBS-LRR genes were the most abundant genes among genes related to signaling receptors, likely because CaLas exist within phloem sieve elements and NBS-LRRs are involved in the direct recognition of CaLas proteins. A slightly greater number of transcripts related to NBS-LRR genes were expressed in the MRL2-11 hybrid than in the MRL2-12 hybrid, suggesting that CaLas might trigger more severe immune responses in the susceptible (MRL2-11) hybrid than in the tolerant (MRL2-12) hybrid. Similar results were observed in genome-wide association studies of HLB-susceptible and HLB-tolerant citrus accessions, where NBS-LRR genes were significantly upregulated in HLB-susceptible citrus accessions [[134]6]. In addition to NBS-LRR, transcripts encoding PRRs, RLPs, RLKs, and CDPKs were significantly upregulated in the tolerant MRL2-12 hybrid in our dataset. These defense patterns against CaLas infection play crucial roles in HLB tolerance. RLK genes play pivotal roles in plant defense mechanisms against pathogen infection [[135]55]. CDPKs play essential roles in diverse plant immune responses, including regulation of the oxidative burst, hormone-mediated signal transduction, and defense gene expression [[136]56]. ATCNGC1, a member of CNGC group 1, is involved in Ca^2+ uptake in plants and root growth [[137]57]. CNGC1 acts as a Ca^2+-permeable channel and is involved in Ca^2+ oscillations and receptor-mediated signaling during plant defense [[138]58]. In particular, Ca^2+ influx has been shown to activate the oxidative burst after elicitor treatment, which causes pathogen- or elicitor-induced cell death in cowpea, tobacco, and soybean [[139]59, [140]60, [141]61, [142]62]. CNGCs induce a hypersensitive response and act as positive mediators in various pathogen tolerance responses involving PR1 protein accumulation and SA signaling pathways [[143]63, [144]64]. Significant upregulation of CNGC and its isoform was also observed in HLB tolerant AtNPR1 transgenic citrus lines, and this gene is also important for PAMP-triggered immunity during pathogen infection in Arabidopsis plants [[145]17, [146]65]. Activating defense response patterns at early stages in the tolerant MRL2-12 hybrid might be responsible for developing rapid and strong responses that limit CaLas growth after infection. Antimicrobial peptides (AMPs) and cysteine-rich peptides are major components of innate immunity in plants [[147]38]. CAP proteins were reported to be abundant in finger lime and HLB-tolerant Australian desert lime [[148]9, [149]66]. This dataset also revealed significant upregulation of CAP transcripts in C. australis hybrids, specifically in MRL2-12. Most genes involved in hormonal crosstalk were expected to be commonly regulated in the two hybrids. The upregulation of several genes involved in JA-mediated signaling confirmed the typical defense responses induced by HLB infection. The pathogenesis of CaLas and its modulation of hormone-mediated defense responses can be blocked to achieve tolerance in the host [[150]67]. Plant proteases are crucial in pathogen perception, defense priming, signaling, and PCD [[151]68]. Papain-like cysteine proteases (PLCPs) positively regulate immunity and defense hormone accumulation. PLCPs are involved in immune responses in phloem-enriched tissue and are targeted by a CaLas effector, SDE1, in citrus [[152]69]. The elevated expression of genes encoding PLCPs suggests their involvement in HLB tolerance in the MRL2-12 hybrid. The first evidence of the role of plant aspartic proteases in biotic stress was reported in tomato leaves, which degraded pathogenesis-related (PR) proteins secreted in response to pathogen attack, preventing their overaccumulation [[153]70]. Aspartic protease is believed to cleave PR-1b, which releases a peptide, CAP-derived peptide 1 (CAPE1), which triggers the expression of genes involved in innate immunity, stress, and defense responses, and SAR [[154]71, [155]39]. CAPE1 is proposed as a novel damage-associated molecular pattern (DAMP) linked to the JA and SA pathways and SAR activation [[156]39]. Reactive oxygen species (ROS) are crucial for plant defense against pathogen attacks. In this dataset, the levels of nine transcripts encoding different peroxidases (peroxidases 64, 55, 48, 42, 22, and 21) were elevated in the transcriptomic analysis (Fig. [157]5). Similar results were obtained for citrus during a proteomic analysis, which revealed abundant peroxidase activity [[158]45]. In melon, redox-related proteins, including peroxidase, are localized in the phloem and significantly increase in abundance during pathogen infection [[159]72]. Similarly, in grapefruit, the levels of ROS-regulating enzymes are elevated at the protein level following CaLas infection [[160]73]. Kunitz-type protease inhibitor, which is located in the cell wall and responsive to peroxidase, was significantly upregulated in this dataset. It also has an inhibitory effect on insects and pathogens [[161]74]. These findings suggest that tolerance is also associated with the induction of ROS-generated genes [[162]46], indicating the importance of the redox balance process during CaLas infection [[163]45]. Cell wall modification defense responses Plant cell walls are composed of cellulose, hemicellulose, and pectins, which play a defensive role against pathogen infection. Significant pectin breakdown is linked to tissue decay, leaf abscission, and swelling of the middle lamella, which has been observed around sieve elements [[164]75, [165]76] during the early stages of CaLas infection. Pectin degradation occurs due to infection with bacteria and fungi, with pectin methylesterase inhibitors (PMEIs) playing a defensive role by regulating cell wall susceptibility to microbial endopolygalacturonases [[166]34]. The expression of PMEI was upregulated in the tolerant hybrid, which may be linked to an interaction with CaLas infection. Xyloglucan endotransglycosylase (XET) is involved in the synthesis of xyloglucan, the main hemicellulose in the cell wall. Transcripts encoding XET are upregulated in C. sinensis leaves but downregulated in stems infected with CaLas [[167]77]. XET is important during fruit ripening, a process that CaLas disrupts. The expression of genes involved in cellulose synthesis, which play essential roles in regulating cell wall structure, was also upregulated in this dataset [[168]78]. COBRA-like genes play crucial roles in cell wall maintenance and biosynthesis. They regulate the orientation of cell wall expansion and modulate cellulose deposition in the root cell wall [[169]79, [170]80]. The tolerant MRL2-12 hybrid exhibited activation and expression of enzymes that strengthen the cell wall, likely reinforcing physical barriers to limit CaLas invasion. Similar to our findings, a previous study revealed that the plant defense response to the Citrus Tristeza Virus (CTV) involved the activation of enzymes required for cell wall strengthening [[171]34]. We observed the elevated expression of the lectin domain, which prevents the pathogen from entering the cytoplasm and is involved in plant defense mechanisms. Lectins are particularly important in defending against indirect interactions of bacteria with cell wall carbohydrates or extracellular glycans [[172]81]. Other indirect mechanisms involve lectin binding to bacterial cell wall peptidoglycans, which exclusively interact with muramic acid, N-acetylmuramic acid, and muramyl dipeptide. Peptidoglycans in bacterial cell walls are crucial for their structural integrity and survival, enhancing plant defenses against microbes [[173]82]. The elevated expression of transcripts, particularly those related to lectin domains, Aps, CAP genes, and others, suggests their involvement in HLB tolerance in the MRL2-12 hybrid. Excessive callose deposition in phloem sieve plates leads to blockage of assimilate transport. This blockage of sieve pores contributes to starch accumulation and the development of HLB symptoms [[174]83]. In this study, the relatively low callose deposition and a higher Pa/Xa ratio in the MRL2-12 hybrid indicate a more tolerant response to HLB. In contrast, the susceptible hybrid MRL2-11 exhibited higher callose deposition. A previous study highlighted that the phloem of Valencia sweet orange produced more callose in response to HLB when compared to a tolerant Finger lime. The study also suggests that the increased callose production in susceptible plants is primarily a mechanical response to HLB, as it blocks the transport of nutrients and other essential molecules but doesn’t prevent the spread of CaLas [[175]9]. Other defense responses in response to HLB infection Carbohydrate metabolism and transporter genes were the most affected biological functions in the hybrid plants. Delayed carbohydrate metabolism responses as have also been observed in Mexican lime might be related to the mechanism of tolerance, avoiding the early phloem blockage caused by callose and starch deposition [[176]5], which might be the reason for the significant upregulation of phloem proteins and the corresponding downregulation of callose synthase in the MRL2-12 hybrid. Regarding nutrient uptake, manganese (Mn) and zinc (Zn) are essential nutrients for plant defense, especially in plants affected by HLB. These nutrients play a critical role in protecting plants from damage caused by ROS [[177]76]. Zn deficiency can disrupt the activity of membrane-bound NADPH oxidase, which impairs the detoxification of H[2]O[2]. This disruption can lead to HLB symptoms such as chlorosis, necrosis, and stunted growth [[178]84]. This dataset revealed the upregulation of the Zn transporter, suggesting that this protein may play a role in alleviating Zn deficiency or scarcity by providing additional zinc to the affected areas of the plant. This response appears to contribute to the tolerance of plants to HLB. Conclusions This study evaluated hybrids between the Marisol clementine and C. australis (MRL), highlighting their diverse biochemical and molecular responses to HLB infection. This study provides insights into the mechanisms of HLB tolerance in C. australis hybrids, as C. australis has been identified as an HLB-tolerant species. Among the MRL hybrids investigated, MRL2-12 exhibited high chlorophyll content and relatively low starch accumulation compared to the other hybrids. This hybrid was further studied to identify molecular determinants of tolerance. Transcriptomic analysis revealed key gene modulation patterns associated with pathogen recognition, signaling, and defense responses in citrus hybrids during CaLas infection. The greater expression of NBS-LRR genes in the susceptible hybrid (MRL2-11) suggests a strong but potentially ineffective immune response. In contrast, the tolerant hybrid (MRL2-12) exhibited upregulation of PRRs, RLKs, CDPKs, and CAP proteins, indicating a more effective defense mechanism. The involvement of plant proteases, JA- and SA-mediated hormonal crosstalk, and ROS-regulating enzymes further highlights the complexity of HLB tolerance. Additionally, the tolerant hybrid exhibited a higher Pa/Xa ratio and regulated callose deposition, along with the upregulation of genes involved in cell wall reinforcement, likely strengthening physical barriers against CaLas invasion. Comparative genomic analysis revealed significant polymorphic variants in MRL hybrids, suggesting a genetically diverse background. These findings suggest that early and coordinated activation of defense mechanisms, including cell wall fortification and immune signaling, contributes to enhanced HLB tolerance. Supplementary Information Below is the link to the electronic supplementary material. [179]Supplementary Material 1^ (526.5KB, pdf) [180]Supplementary Material 2^ (238.3KB, pdf) Abbreviations CaLas Candidatus Liberibacter asiaticus CNGC Cyclic nucleotide-gated ion channel DEGs Differentially expressed genes GO Gene ontology HLB Huanglongbing INDEL Insertion or deletion LRR Leucine-rich repeat LTP Lipid transfer protein PAMP Pathogen-associated molecular patterns PMEIs Pectin methylesterase inhibitors PR Pathogenesis-related RLKs Receptor-like kinases RLPs Receptor-like proteins ROH Runs of homozygosity ROS Reactive oxygen species SNP single nucleotide polymorphism WAKs Wall-associated kinases Author contributions SR - wrote the manuscript, conducted experiments and analysis; LM- performed RNA extraction; JD-performed hybridity validation; SW-microscopic analysis; MD-generated the hybrids, designed the study, obtained funding, reviewed the manuscript, and supervised the project. Funding This work was supported by the USDA-ECDRE project award no. 2021-70029-36055 from the U.S. Department of Agriculture’s National Institute of Food and Agriculture. Data availability The datasets presented in this study are available in the SRA repository with bioproject accession number PRJNA1245434. Declarations Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. References