Exploring the Genetic Basis of Heterosis in Calonectria spp. Resistance on Eucalypt

1. Introduction

1.1. Background and Significance of Research on Eucalypt Leaf Blight

Eucalypt leaf blight, a significant global epidemic, is triggered by fungi from the Calonectria genus (Ca.) of the Nectriaceae family [1]. Eucalypt (Eucalyptus, Corymbia and Angophora) industrial plantations in Guangdong, Guangxi, Fujian, and Hainan are generally affected by eucalypt leaf blight [2]. Based on the model estimated by Zhu Jianhua et al. (2011), the average annual loss caused by eucalypt leaf blight in Fujian amounts to 51.8 million RM [3]. Therefore, eucalypt leaf blight has become one of the serious factors affecting the investment return and sustainable development of eucalypt plantations in China. The species diversity of the pathogenic bacteria in the genus Calonectria is rich, and 171 species have been reported [4,5]. Six pathogens pose a serious threat to eucalypt plantations in Guangdong and Fujian [6]. The pathogens responsible for eucalypt leaf blight are diverse and variable. Given the diverse and expansive artificial forest environment, as well as the tall and open crown of the trees, implementing effective biological and abiotic control measures for that is a huge challenging. The only solution is to select new resistant genotypes for specific environments and fundamentally solve eucalypt leaf blight.

1.2. About Eucalyptus pellita

The high-temperature and humidity summer season in South China coincides with the outbreak of eucalypt diseases. The superior disease resistance of the hybrid of E. pellita ensures that its growth remains largely unaffected, and accentuating its rapid growth [7]. Our multi-site testing results on hybrid of E. pellita indicated that its rapid growth rate is not inferior to that of hybrid from E. grandis × E. urophylla, and its disease resistance is significantly higher than that of hybrid of E. grandis × E. urophylla in areas severely affected by diseases such as eucalypt leaf blight. For example, the superior hybrid (EC150) of E. urophylla × E. pellita, at 45 months, had an average DBH of 13.7 cm, average height of 18.4 m, and an average volume of 138.2 m³/ha, it is 52.9% -81.5% higher than the control DH3229 in the western Guangdong.

E. pellita has excellent disease resistance and rapid growth, and prefers environments with strong sunlight, humid heat, and abundant rainfall [8,9,10]. E. pellita and its hybrids are one of the ideal tree species for enriching eucalypt plantations species in southern China [11]. Moreover, E. pellita has good hybrid cross-compatibility and trait-complementary with many other eucalypt species [12,13]. It is an important interspecific hybrid parent of eucalypt in South China. In 2007, we evaluated the early introduction (from Australian and Indonesia) of E. pellita provances and families trials in China systematically, and based on the evaluation, more than 200 excellent families were selected and a advanced generation of breeding population for E. pellita were established [14,15].

1.3. Varieties Improvement is the Way to Solve Disease Problems

Selecting parents with different genetic foundations for artificial hybridization and utilizing heterosis is an important way to improve plant growth, biomass, stress resistance, and adaptability [16]. However, the formation mechanism of heterosis in target traits of E. pellita, especially the disease resistance heterosis, has not been studied yet. Therefore, exploring the molecular mechanisms of genetic recombination and gene expression patterns related to the heterosis of E. pellita resistant to eucalypt leaf blight, and exploring gene information related to eucalypt disease resistance, would improve the development efficiency of new disease resistant genotypes, and that’s the most effective way to solve the problems of poor stress resistance (such as eucalypt leaf blight and typhoon resistance) and limited varieties used in eucalypt plantations.

Conducting a transcriptional analysis of differentially expressed genes across various genotypes will facilitate the exploration of the genetic mechanisms underlying heterosis in eucalypt resistance to eucalypt leaf blight. The analysis of transcription levels plays a crucial role in comprehending gene expression patterns, gene architecture, and function. Transcriptome sequencing (RNA sequencing, RNA-Seq)enables the accurate acquisition of comprehensive transcriptional information and high reproducibility [17]. Furthermore, it allows for the detection of subtle transcriptional level differences, thus facilitating a comprehensive and in-depth analysis of differentially expressed genes and their expression patterns between the hybrids and parents. This approach can contribute to a better understanding of the genetic basis underlying the formation of heterosis [18]. In the current study, two eucalypt hybrids, susceptible and resistant, were screened by leaf blight pathogen testing, which had conducted in our precious study. Transcriptome sequencing for two hybrids and their parents was conducted using RNA-Seq, so as to Explore the genetic basis for the formation of heterosis of resistant to eucalypt leaf blight through gene differential expression analysis, and to explore the relevant genes of resistant to eucalypt leaf blight tentatively, finally, providing a theoretical basis for hybrid breeding of eucalypt resistant to leaf blight.

2. Materials and Methods

2.1. Genetic Material

Eucalypt hybrid leaf blight susceptibility testing was conducted via spray inoculation with spore suspension [5]. After 72 hours of inoculation and infection, the susceptibility index (DI) of leaf blight caused by the fungus of the genus Calonectria on the leaves of hybrid eucalyptus was calculated. The percentage of susceptible area of the tested leaves was scanned using the software Leaf Doctor, and the degree of susceptibility was assigned a value from 0 to 5 based on the percentage. The susceptibility index of the tested eucalypt hybrids was calculated after the statistics were completed [19]. This part has been completed in our previous research.

Based on our previous eucalypt leaf blight susceptibility test results [20], eucalypt leaf blight disease-susceptible hybrid EC333 and the disease-resistant hybrid EC338 were screened, and their corresponding parents were found, for a total of 5 eucalypt genotypes used for transcriptome analysis (Table 1). There were three biological replicates for each genotype. During the peak growing season of eucalypt, healthy young leaves with consistent physiological development at the top of the trunk were collected and quickly placed in liquid nitrogen.

2.2. Transcriptome Sequencing and Analysis

RNA Extraction: RNA was isolated utilizing the TRIzol kit (Tiangen, Beijing) and mRNA purification was carried out through oligo-T-linked magnetic beads [21]. The assessment of RNA integrity was conducted with precision employing the Agilent 2100 bioanalyzer, utilizing RNA integrity and total quantity as key reference parameters [22].

Library building and sequencing: First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase(RNase H-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ ends of DNA fragments, Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 370~420 bp in length, the library fragments were purified with AMPure XP system [23]. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Novaseq platform and 150 bp paired-end reads were generated [24].

Gene Differential Expression Analysis: Raw data (raw reads) of fastq format were firstly processed through fastp software. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30 and GC content the clean data were calculated. All the downstream analyses were based on the clean data with high quality. Reference genome and gene model annotation files were downloaded from genome website directly. Index of the reference genome was built using Hisat2 v2.0.5 and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5 [25]. The mapped reads of each sample were assembled by StringTie (v1.3.3b) [26] in a reference-based approach. FeatureCounts v1.5.0-p3 was used to count the reads numbers mapped to each gene. And then FPKM of each gene was calculated based on the length of the gene and reads count mapped to this gene [27]. (For DESeq2 with biological replicates) Differential expression analysis of two conditions/groups (two biological replicates per condition) was performed using the DESeq2 R package (1.20.0). DESeq2 provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P-value <=0.05 found by DESeq2 were assigned as differentially expressed. (For edgeR without biological replicates) Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edgeR program package through one scaling normalized factor. Differential expression analysis of two conditions was performed using the edgeR R package (3.22.5). The P values were adjusted using the Benjamini & Hochberg method. Corrected P-value of 0.05 and absolute foldchange of 2 were set as the threshold for significantly differential expression [28].

The correlation of gene expression levels among different eucalypt genotypes is a significant indicator in assessing experimental reliability and sample selection. A higher correlation coefficient indicates greater similarity and closeness in expression patterns. The Encode Project suggests a Pearson correlation coefficient R² > 0.92 (under optimal sampling and experimental conditions) (ENCODE Project Consortium, 2004). Besides genotype correlation, the correlation between biological replicates during testing (R² > 0.8) is also vital for the reliability of differential gene analysis. By calculating correlation coefficients based on FPKM values and visualizing them in a heat map, differences between sample groups and duplications within groups can be visually displayed.

Gene Ontology (GO) enrichment analysis of differentially expressed genes was implemented by the clusterProfiler R package, in which gene length bias wascorrected. GO terms with corrected Pvalue less than 0.05 were considered significantly enriched by differential expressed genes. We used clusterProfiler R package to test the statistical enrichment of differential expression genes in KEGG (http://www.genome.jp/kegg/) pathways (Table S2, Figure S1) [29]. The genes were categorized based on their GO annotations according to the subnodes of the three major GO classifications (TERM) (Figure S2).

We use the local version of the GSEA analysis tool http://www.broadinstitute.org/gsea/index.jsp, GO, KEGG data set above were used for GSEA independently. In this investigation, the gene list of E. grandis was retrieved from the database to establish a network, while rMATS (4.1.0) was utilized to scrutinize alternative splicing occurrences.

Alternative Splicing is an important mechanism for regulate the expression of genes and the variable of protein. rMATS(4.1.0) software was used to analysis the AS event. GATK (v4.1.1.0) software was used to perform SNP calling. Raw vcf files were filtered with GATK standard filter method and other parameters (cluster:3; WindowSize:35; QD < 2.0 ; FS > 30.0; DP < 10.

WGCNA (Weighted correlation network analysis) is a systematic biological method used to describe the gene association modes among different samples. It can be used to identify gene sets that are highly synergistic changed, and identify candidate biomarkers or therapeutic targets based on the coherence of gene sets and the correlation between gene sets and phenotypes. The R package WGCNA is a set of functions used to calculate various weighted association analysis, which can be used for network construction, gene screening, gene cluster identification, topological feature calculation, data simulation and visualization. One input file is sample information, that is, a matrix describing the traits of the sample: the traits used for association analysis must be numeric; If it is a regional or categorical variable, it needs to be converted to a 0-1 matrix. The other is gene expression data. For transcriptome sequencing, FPKM can be used as gene expression data [30].

The disease resistant eucalypthybrid is EC338, and the susceptible hybrid is EC333. The paternal parent of EC338 is P9060 (E. pellita), and the maternal parent is W1767 (E. wetarensis). The maternal parent of EC333 is H1522(E. urophylla × E. pellita), the paternal parent of EC333 is unknown, so it was not sequenced together. The target disease resistant genes are most likely to be differentially expressed genes between disease resistant genotype and susceptible genotype, and the disease resistant genes can be traced back to their parents. Therefore, in transcriptome analysis, this study mainly compared the differential gene expression information between disease resistant and susceptible eucalypt hybrids, as well as between hybrids and their parents.

3. Results

3.1. Quality Control of Sequence Data (QC)

3.1.1. Sequencing Data Quality

After performing initial data filtering (Figure S3), examining sequencing error rates, and analyzing the distribution of GC content, high-quality clean reads were obtained for further analysis. From the correspondence between base identification and Phred scores, the evaluation of base calling accuracy in eucalypt samples revealed a Q20 range of 97-98% and a Q30 range of 93-94% (Table S1), indicating reliable data quality with an over 90% correct base identification rate. The sequencing error rate was found to be low, ranging from 0.02 to 0.03, and the GC content distribution was stable around 50% for all samples (Table S1). The results affirm the dependability of the sequencing data and the appropriateness of the refined reads for subsequent analysis.

3.1.2. Mapping Sequencing Information to E. grandis Reference Genome and Its Regional Distribution

Sequencing fragments were randomly interrupted by mRNA. To determine which genes transcribed these fragments, it is necessary to compare the clean reads after quality control to the reference genome of E. grandis [31], a species similar to the material studied in this paper. The transcript sequencing FPKM is paired-end sequencing. A gene is represented by a pair of reads. Perform subsequent quantitative data analysis using the number and percentage of reads aligned to the unique map of the reference genome of E. grandis. The unique map reads of EC338 were 37 971 737. The unique map reads of EC333 were 34 174 054. The clean reads for P9060 and W1767 were 34 581 445 and 37 443 452, respectively. The clean reads for H1522 were 36 358 616 (Table S2).

The proportion of reads in exonic regions, intronic regions, and intergenic regions of the genome were counted separately. Generally model species with better gene annotation (e.g., human and mouse) have a high proportion of mapping to exonic regions. The distribution of sequencing reads across genomic regions for all eucalypt genotypes in this paper showed that not many reads were mapped to intergenic regions (Figure S4), which might be due to contamination by ncRNAs or a few DNA fragments, or the gene annotations might not be well developed. The reads mapping to the reference genome might have originated from precursor mRNAs or introns stranded by variable splicing events.

3.2. Quantitative Analysis of the Genes Sequenced

The number of reads covered by each newly predicted gene from the start to the end region were counted based on the positional information of the gene mapping on the reference genome of E. grandis. The reads with mapping quality value less than 10, reads on non-composite mapping, and reads with mapping to multiple regions of the genome were filtered out. The quantification of gene expression level was performed for each sample separately and then combined to obtain the expression matrix of all samples [32] to obtain the information of detected genes in eucalypt parents and their hybrids (Table S3).

3.2.1. Detection of Expression Distribution of Genes Sequenced in Different Eucalypt Genotypes

Gene expression quantification using FPKM (expected number of fragments per kilobase of transcript sequence per million base pairs sequenced) is the number of fragments per million fragments per kilobase length from a gene. This metric takes into account the effects of sequencing depth and gene length on the count of fragments and converts the read count number into FPKM. It plots the overall expression distribution of the samples to view the expression level of the samples and compare the overall expression distribution of the different samples. Box plots can be utilized to illustrate the dispersion of a set of data (Figure 1A). Violin plots demonstrate the expression distribution and its probability density, with the inflated portion of the plot indicating the most concentrated region of expression in the entire sample (Figure 1C).

The results of three repeated biological tests indicated that the distribution of gene expression levels among different genotypes were similar. With the exception of a slight difference in the parental W1767, the median and the degree of dispersion of the distribution of gene expression levels of the four genotypes of hybrids EC338, EC333, P9060 and H1522 were relatively similar (Figure 1A, Figure 1C). Therefore, the potential for differential gene screening due to the amount of gene expression can be ruled out. The peak of the FPKM density distribution curve is commonly used to analyze the gene expression level in RNA sequencing data. The peaks of the curves represent the most common FPKM values in the dataset and can provide information about gene expression patterns. FPKM values were calculated by dividing the number of fragments mapped to a gene by the length of the gene and the total number of mapped reads, and then normalizing to the library size.

The density distribution curve is a graphical representation of FPKM values. The peaks of the curves represent the highest frequency FPKM values, indicating the most common gene expression levels in the samples. The distribution of the curves for the five genotypes were generally consistent, with most genes having a log10 (FPKM) between 0 and 5. The first and highest peak of the curve is in close proximity to the origin; however, the density plot in this study had two peaks, indicating that two regions with the highest concentration of gene expression occurred in five genotypes (Figure 1B). This also reflected the two expanded portions of Figure 1C, and in fact Figure 1C was a reexamination and summary of the basis of Figure 1A and Figure 1B.

3.2.2. The Correlation of Gene Expression of Eucalypt Genotypes Sequenced

The gene expression correlation heat map analysis conducted in this study (Figure 2) revealed that the eucalypt hybrid EC333 and EC338 exhibit the highest correlation coefficient (R² > 0.9), indicating a strong similarity in their expression patterns and close relationship. The high correlations were observed among the three replicate groups of each genotype enhance the reliability of subsequent differential gene analysis.

3.2.3. Principal Component Analysis of Gene Expression

In this study, PCA was carried out on the gene expression values (FPKM) of all sequenced eucalypt genotypes (Figure 3). The expected outcome of PCA is to have samples from different groups dispersed and samples within groups clustered together in the plot. The PCA plot in Figure 3 demonstrated that the five eucalypt genotypes formed distinct clusters, reflecting their genetic relatedness. The eucalypt hybrids were clustered closely together, with a smaller distance between them compared to their distance from the parental genotypes.

3.3. Gene Co-Expression Venn Diagram and Differential Gene Venn Diagram

There were 18 107 differential genes between EC338 and EC333 (Figure 4A), 15 969 differential genes between EC338 and W1767 (Figure 4B). A total of 17 351 differential genes were identified between EC338 and P9060 (Figure 4C). Therefore, target resistance genes might exist in genes shared by these three gene sets. In addition, genes can interact additively or reciprocally, causing the function of recessive disease susceptibility genes in a parent to become visible in their hybrids. The female parent (H1522) of EC333 is not susceptible to leaf blight, while the male parent was unknown, so its susceptibility might be attributed to gene recombination in cross. Then the 17 411 differential genes between EC333 and H1522 (Figure 4D) were hardly to say associated with susceptibility or not. The four pairs of comparison groups described above were distinguished by a common differential gene of 2143 (Figure 5), indicating that these genes were likely associated with resistance to eucalypt leaf blight.

3.4. Gene Expression Analysis (Differential Gene Screening and Clustering)

The results indicated that 1631 genes were up-regulated, while 2281 genes were down-regulated (Figure S5). In the volcano plot, the horizontal axis represents the difference in gene expression, with up-regulated genes depicted to the left of the dashed line as red dots and down-regulated genes as green dots. For samples with biological replicates, the adjusted p-value (padj) must be less than 0.05, which means that both red and green dots in the volcano plot must be above 1.301. Conversely, the smaller the corrected p-value is, the larger the corresponding log10 (p-value) value is, indicating that the more significant the difference between the genes in the comparison combination, and thus the more deviated from the horizontal axis upward, the more worthy of our attention. Among these genes that deviate up the horizontal axis, we need to find the common differential genes in the four pairs of comparative combinations mentioned above (Figure S6).

In this study, the FPKM values of genes are clustered and analyzed using mainstream hierarchical clustering, and genes or samples with similar expression patterns in the heat map are clustered together. The color in each square reflects not the gene expression value, but the value obtained after homogenization (Z-score) of the rows of expression data (range in -3 to 3). Vertical comparison refers to the expression of different genes in the same sample. Horizontal comparison refers to the expression of the same gene in different samples.

The expression results among different eucalypt genotypes revealed that numerous genes exhibited highly significant differential expression in W1767 (E. wetarensis). When the expression level of the differential genes were low in W1767, it is highly expressed in EC338 (Figure 6). All of these genes with large differences in expression between parents and hybrids were likely to be associated with disease resistance and could be the focus for subsequent screening of resistant genes.

3.5. Enrichment Analysis of Differential Genes in Comparison Groups

Given the significant differences in susceptibility to eucalypt leaf blight, the primary focus should be on the differentially expressed genes of the resistance genotype EC338 and the susceptibility genotype EC333, which should be the focus of the analysis of eucalypt leaf blight resistant genes. Furthermore, in order to investigate the inheritance pattern of the parents, the genes shared in the differential expression pathway between EC338 and P9060, EC338 and W1767, and EC333 and H1522 were also examined. This approach was employed in order to analyze the reasons for the formation of disease resistance heterosis, and attempt to identify the key genes related to eucalypt leaf blight resistance.

The GO and KEGG enrichment analysis was conducted to identify the pathways of differentially expressed genes in the comparison group. It is important to note that only pathways with ‘details’ can be extracted to specific genes in Gene Set Enrichment Analysis (GSEA). Therefore, we have primarily selected the pathway genes with ‘details’ after enrichment as the basis for determining the genes related to eucalypt leaf blight resistance.

The GO enrichment analysis indicates that the differentially expressed gene pathways of EC338 (resistant) and EC333 (susceptible) predominantly involve carbohydrate binding, pattern binding, and polysaccharide binding (Table 2 and Figure 7A). These pathways might be key to endowing hybrids with disease resistance advantages, with carbohydrate binding demonstrating the highest probability, because this pathway existed in four pairs of comparison groups, and only carbohydrate binding had details. Terpene synthase activity and Carbon-oxygen lyase activity, acting on phosphates were present in three pairs of genotypically differentially expressed pathways (Table 2, Figure 7B, Figure 7C and Figure 7D). Carbohydrate binding (Table 2, Figure 7A and Figure 7D) and carbon-oxygen lyase activity (Table 2, Figure 7B and Figure 7D) were present in two of the genotype pairs compared.

Figure 7. Bar (top) and bubble (bottom) plots of the GO analysis for four pairs of differential genes set.

Figure 7. Bar (top) and bubble (bottom) plots of the GO analysis for four pairs of differential genes set.

Figure 8. Bar (top) and bubble (bottom) plots of KEGG analysis for four differential gene pairs.

Figure 8. Bar (top) and bubble (bottom) plots of KEGG analysis for four differential gene pairs.

The Gene Set Enrichment Analysis (GSEA) enrichment analysis utilizing GO annotation primarily focuses on carbohydrate binding (GO: 0030246) (Figure 9A) of EC333 in EC338 VS EC333, as well as terpene synthase activity (GO: 0010333) (Figure 9B), carbon-oxygen lyase activity, acting on phosphates (GO: 0016838) (Figure 9C) and carbon-oxygen lyase activity (GO: 0016835) (Figure 9D) of EC338 in EC338 VS P9060. The specific differentially expressed genes that were the focus of this pathway were those listed in Table S4, Table S5, Table S6, and Table S7, respectively.

Figure 9. GSEA enrichment analysis using GO annotations (The top row is Snapshot of enrichment results. The Enrichment Map displays the green ES line graph at the top, revealing the highest peak at the ES value denoting the gene set ES value. In the hits barcode map located in the middle, the leading edge subset represents the genes that make the most significant contribution to the enrichment score, ranging from 0 to the ES value. The top grayscale enrichment of the bottom sorting level value indicates the level of genotype expression. The Distribution plot in the bottom row illustrates that if the final ES value (marked by the black dotted line) is relatively small, the pathway gene set is randomly distributed, whereas a relatively large ES value suggests non-random distribution. Subsequently, the hypothesis regarding substitution is then examined.).

Figure 9. GSEA enrichment analysis using GO annotations (The top row is Snapshot of enrichment results. The Enrichment Map displays the green ES line graph at the top, revealing the highest peak at the ES value denoting the gene set ES value. In the hits barcode map located in the middle, the leading edge subset represents the genes that make the most significant contribution to the enrichment score, ranging from 0 to the ES value. The top grayscale enrichment of the bottom sorting level value indicates the level of genotype expression. The Distribution plot in the bottom row illustrates that if the final ES value (marked by the black dotted line) is relatively small, the pathway gene set is randomly distributed, whereas a relatively large ES value suggests non-random distribution. Subsequently, the hypothesis regarding substitution is then examined.).

The provided table illustrates that the four differentially expressed gene pathways satisfy the criteria of |NES|>1, NOM p-val<0.05, and FDR<0.25, suggesting promising enrichment outcomes deserving further investigation (Table 3).

The KEGG enrichment analysis revealed that the differentially expressed gene pathways of EC338 (resistant) and EC333 (susceptible) contained phenolylpropanoid biosynthesis, flavonoid biosynthesis, Plant-pathogen interaction, and Sesquiterpenoid and triterpenoid biosynthesis (Table 2 and Figure 8A). These pathways were likely linked to susceptibility to eucalypt leaf blight. Phenylpropanoid biosynthesis (EGR00940) was present in all four pairs of genotype differentially expressed gene pathways, while Sesquiterpenoid and triterpenoid biosynthesis (EGR00909) was found in three pairs of genotype comparisons (Table 2, Figure 8A, Figure 8B and Figure 8C). Flavonoid biosynthesis (EGR00941) and Plant-pathogen interaction (EGR04626) pathways (Table 2 and Figure 8) were present in two pairs of genotype comparisons(Table 2, Figure 8A and Figure 8C).

The GSEA enrichment with KEGG annotation in Table 4 suggested that the main focus is on sesquiterpenoid and triterpenoid biosynthesis (EGR00909) (Figure 10A) of EC338 in EC338 VS EC333, plant-pathogen interaction (EGR04626) (Figure 10C) of W1767 in EC338 VS W1767, as well as phenylpropanoid biosynthesis (EGR00940) (Figure 10B) and Flavonoid biosynthesis (EGR00941) (Figure 10D) of EC338 in EC338 VS W1767 these four pathways. The specific differentially expressed genes that were the focus of this pathway were those listed in Table S8, Table S9, Table S10, and Table S11, respectively.

3.6. Alternative Splicing Event Analysis

The current study identified exon skipping events (SE) with the most significant IncLevelDifference values in the two samples from a pool of the top 20 with the lowest False Discovery Rate (FDR) in four comparison groups: EC338 VS EC333, EC338 VS P9060, EC338 VS W1767, and EC333 VS H1522. It was concluded that the gene LOC104426602 exhibited differential alternative splicing between EC338 and EC333 (Figure 11A), as well as LOC104435785 between EC338 and P9060 (Figure 11B), and LOC104448051 between EC338 and W1767 (Figure 11C). Among various genes with an IncLevelDifference absolute value of 1 in EC333 VS H1522, the exon skipping isoform LOC104418168 was the only one expressed in this differential comparison group (Figure 11D). It is worth noting that LOC104426602, LOC104435785, and LOC104418168 had high count of exon skipping isoforms, whereas LOC104448051 had number of exon inclusion isoforms. The substantial differences in alternative splicing events between the two opposite resistant hybrids could result in variations in transcript and protein structures, as well as functional polymorphisms. Consequently, it was possible that the mentioned four genes are associated with the resistance to eucalypt leaf blight.

3.7. SNP Variation Loci Analysis

P9060 and W1767 exhibited a lower number of silent mutation loci (silent) and missense mutation loci (missense), with counts of 95917 and 98045, respectively. In contrast, other genotypes demonstrated a higher prevalence of functional mutation loci. Missense mutations lead to the substitution of various amino acids, potentially affecting the functionality of the resultant protein. Here we primarily focused on the 2988 distinct missense mutation loci between EC338 (resistant) and EC333 (susceptible). Notably, EC338 showed discrepancies of 12031 and 9903 missense mutation loci compared to its parents (Figure 12A), implying a potential for the emergence of disease resistance-related genes.

INDELs usually cause more significant genetic modifications than SNPs because they involve a larger number of nucleotide sequences. INDELs were used to investigate variant loci between the disease-resistant hybrid EC338 and its parents in present study. It was worth noting that EC338 had 626 more impactful mutation sites than P9060, and 396 more than W1767 (Figure 12B).

The TRANSCRIPT region had the lowest number of nucleotide variations for both INDEL and SNP. Nonetheless, scrutiny of the SNP variation loci demonstrated that the EXON region displayed the highest number of variations across all genotypes. The INDEL variation analysis identified the DOWNSTREAM region had great number of variations. The analysis indicated that EC338 possessed 9631 more mutation loci in the EXON region compared to EC333. Furthermore, the female parent W1767 exhibited 259 additional mutated loci in the DOWNSTREAM region in comparison to EC338, whereas the male parent P9060 displayed 3107 fewer mutated loci in the DOWNSTREAM region than EC338. Gene loci in these regions might be related to eucalypt leaf blight and warrant further exploration.

3.8. Weighted Gene Co-Expression Network Analysis

By setting the correlation coefficient to 0.8 and the soft threshold to 5(Figure 13A), the network was constructed in accordance with the predetermined soft threshold, and the hierarchical clustering of genes within the network was conducted using the dissimilarity between genes to generate a hierarchical clustering tree. Subsequently, the tree was divided into modules using dynamic shearing. The minimum number of genes in the module was 30, merged modules with correlation coefficients greater than 0.75 and anisotropy coefficients less than 0.25 (Figure 13B). Clustered 30 modules based on their eigenvalues. The extraction of the eigenvector genes of each module as the first principal component gene of a particular module also represents the overall level of gene expression within that module. Then performed hierarchical clustering on the characteristic genes of different modules. The positive correlation between MEdarkgrey and MEbrown is the strongest, with a correlation coefficient of 0.55 and a P-value of 0.03 between the eigenvalues of the two modules. The negative correlation between MEdarkgrey and MEgrey is the strongest, with a correlation coefficient of -0.53 and a P-value of 0.04 for the eigenvalues of the two modules. The P-values of the eigenvalues for both pairs of modules were below 0.05, indicating statistical significance.

The two modules with the strongest positive correlation with EC338 were MEpink and METan, with coefficient of 0.485301597 and 0.4389636, respectively. Conversely, the modules MEsaddlebrown and MEred demonstrated the most negative correlations, with coefficient of -0.204643774 and -0.17722817. When considering EC333, the modules MElightgreen and MEsaddlebrown displayed the strongest positive correlations at 0.532894171 and 0.320682081, while modules MEred and MEgreenyellow revealed the most negative correlations at -0.177673509 and -0.170610241 (Figure 14). In theory, the same gene in the module with opposite correlation to the two hybrid species was the gene related to eucalypt leaf blight.

It is crucial that genes shared in modules positively correlated with EC338 and negatively correlated with EC333 Genes shared in modules negatively correlated with EC338 and positively correlated with EC333 were of interest. In this way, of these comparisons, only the MEsaddlebrown module fulfilled the criterion, so only genes 104448928 and others in this module were of interest (Table S12).

4. Discussion

The screening of susceptible and resistant hybrids in this study was conducted through two infection tests on candidate eucalypt hybrids using three Calonectria pathogens. Based on the results of the two infection tests, the resistant genotype EC338 and the susceptible genotype EC333 were screened [20]. Although susceptibility tests had shown differences in the susceptibility index of eucalypt hybrids to different Calonectria pathogens, the relative value of the susceptibility index of the same eucalypt hybrid after being infected twice by three pathogens was basically the same. This basically eliminates the randomness of the susceptibility test results, indicating that the difference in susceptibility of the tested hybrid was due to its inherent genetic differences, that was, the disease resistance was strongly controlled by genetics. Although heterosis is generally controlled by multiple genes and the mechanisms are very complex, the genetic basis for the formation of heterosis varies among different plants, traits, and populations with different genetic backgrounds [33,34]. However, by tracking the parents of disease resistant hybrids and conducting scientific parent combination, it is highly possible to achieve the goal of improvement of disease resistant eucalypt varieties.

Studies had shown that the expression or transgenosis of Carbohydrate binding Modules (CBMs) could change cell wall properties and regulate plant growth and development. Such as Keadtidumrongkul et al. found that CBM2a promoted the growth of Arabidopsis thaliana, and the growth promotion of E. camaldulensis was reflected in the increase of plant height, xylem area, xylem fibers, and conduit cells, there was no significant growth-promoting effect on tobacco [35]. But there are also studies indicating that this is not always the case. Such as the overexpression of certain CBMs did not result in enhanced growth of N.tabacum [36]. And in present study, the enrichment analysis based on GO annotation revealed that the differentially expressed gene pathway between the disease-resistant genotype EC338 and the susceptible genotype EC333 was associated with carbohydrate binding (GO: 0030246), indicating that CBMs might also have a correlation with leaf blight resistance in eucalypt.

In this study, GSEA enrichment based on GO annotation revealed a differential pathway of terpene synthetase (TERPENESYNTHASEACTIVITY) in the disease resistant hybrid EC338 compared to its parent E. pellita P9060. Terpenes are critical components of plant defense against herbivores, insects, pathogenic fungi, and bacteria [37,38]. In many forest tree species, terpenes are vital chemical, and physical defense compounds [39,40]. Conifer terpenes are a primary component of a durable, quantitative defense mechanism against stem boring insects and fungal pathogens [41]. Constitutive and induced terpenes are critical for resistance [42,43,44]. This finding might indicate a certain correlation between terpene synthase and resistance to eucalypt leaf blight.

Furthermore, two differential pathways between EC338 and its parent P9060 were found, namely, the carbon-oxygen cleavage activity on phosphates (CARBON-OXYGEN-LYASE-ACTIVITY_ACTING-ON-HOSPHATES) and carbon-oxygen lyase activity (CARBON-OXYGEN-LYASE-ACTIVITY). A study shown that Cold-resistant DEGs-DEPs had abundant hydrolase activities that could act on glycosyl bonds, carbon-oxygen lyase activity, and iron binding [45]. The MCRA proteins from carbon-oxygen lyase family could be classified as fad-containing double bond hydratases and play a role in the virulence of Streptococcus pyogenes M49 [46]. Based on the results of current and other study, it could be found that the above pathways were indeed related to the stress resistance mechanism of plants.

The current study, GSEA enrichment based on KEGG annotation revealed differential pathways in sesquiterpene and triterpenoid biosynthesis (SESQUITERPENOID-AND-TRITERPENOID-BIOSYNTHESIS) between the resistant hybrid EC338 and the susceptible hybrid EC333. Studies had shown that this pathway was related to the stress resistance mechanism of plants. Such as the study on the molecular response mechanism to drought stress in Atractylodes chinensis had revealed 15 single gene encoding enzymes in the biosynthesis pathways of sesquiterpenes or triterpenes [47]. And a finding done by Singh et al. suggested a relationship between sesquiterpenes, phenolics, flavonoids, and the antioxidant capacity of Ferula plants [48]. Moreover, a study had proposed that the synthesis of sesquiterpenes, triterpenes, flavonoids, and phenylpropanoids might play a crucial role in the response to mechanical damage in Aquilaria sinensis [49].

GSEA enrichment based on KEGG annotations also revealed a differential pathway for phenylpropanoid biosynthesis (PHENYLPROPANOID_BIOSYNTHESIS) in the disease-resistant hybrid EC338 compared to its parent. Gan et al. determined that phenylpropanoid biosynthesis might play an important role in salt tolerance in mulberry at proteomic level [50]. Additionally, a research indicated that Pichia galeiformis enhances resistance to diseases and reduces disease incidence and lesion diameters in chalcone synthase by promoting the activation of genes involved in phenylpropanoid biosynthesis [51].

KEGG enrichment revealed that there is a differential plant pathogen interaction pathway (PLANT_PATHOGEN_INTERACTION) between the resistant hybrid EC338 and the susceptible hybrid EC333. Under biotic stress, plants trigger double immunity against pathogens through plant-pathogen interaction mechanisms, i.e. pathogen-associated molecular pattern-triggered immunity (PTI) and effector-triggered immunity (ETI) [52]. Recent research indicates that PTI and ETI work together to enhance plant immunity against pathogens, thereby enhancing the overall defense system of plants [53,54]. Nonetheless, pathogens have developed sophisticated mechanisms to circumvent plant defenses, leading to successful infection by exploiting host susceptibility [55,56]. Whether it involved plant autoimmunity or pathogen evasion, immunity was intricately linked to the plant-pathogen interaction pathway.

In addition, exploring the causes of plant disease resistance heterosis can also evaluate with phenotypic characteristics such as leaf surface wax, cuticle layer, morphological structure, and chlorophyll content. Recent study conducted by Bonora et al. reveals that resistance mechanism of Corymbia citriodora to different pathogens might be related to leaf phenotypic traits, such as the deposition of polyphenols and tannins on the upper/lower epidermis, and the proportion of monoterpenes, steroids, monounsaturated hydrocarbons, and long-chain hydrocarbons [57]. And, increasing the chlorophyll content in cotton can promote photosynthesis and enhance the plant's ability to resist diseases [58]. So, in the subsequent research of this study, in addition to verifying the identified functional genes, attempts could also be made to study the differences in leaf phenotype characteristics between disease resistant and susceptible genotypes, and analyze the reasons for their disease resistance differences from multiple perspectives.

The susceptibility test in the early stage of this study was conducted under controlled conditions on seedlings aged 6 months, using the pathogen of Calonectria to artificially infect and select susceptible hybrids. However, in the 1.5-year field experiment, there was no significant difference in susceptibility between the two hybrids. This might be related to seedling age and environment. As seedling age increases, the temporal expression of genes causes protein structure modification or isomerism, or different proteins interact to resist adverse environments [59].Yang et al. (2022) illustrated, through comparative transcriptomic and proteomic analysis that P. notoginseng mediated resistance to pathogens by regulating the lignin biosynthesis pathway [60]. Therefore, the protein products most likely associated with the gene pathway associated with eucalypt leaf blight discovered in this study can be isolated, and substrate specificity and kinetic parameter analysis could be performed to further explore the formation mechanism of leaf blight resistance.