The fresh leaves of A. spinulosa and A. metteniana were gathered in the afternoon of 20 December 2020, at Chishui alsophila national nature Reserve (109°45'∼ 106° 03′n, 28°23'∼ 28° 27′e), Chishui, Guizhou Province, P. R. China. The samples were placed in microcentrifuge tubes, rapidly frozen in liquid nitrogen, and stored at −80°C until metabolite and RNA isolation. Three biological replicates were included for each plant.

Total RNA was extracted using TRIzol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. Degradation and contamination of RNA were verified using 1% RNase-free agarose gel electrophoresis, and the purity and integrity of RNA was assessed on an Agilent 2,100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). High-quality RNA (RNA Integrity Number [RIN] scores >7.5) was used for subsequent experiments. The messenger RNA (mRNA) was enriched using Oligo(dT) beads, followed by fragmentation into short fragments using fragmentation buffer and reverse transcribed into cDNA with random primers. Second-strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP and buffer. Subsequently, the cDNA fragments were purified with QiaQuick PCR extraction kit (Qiagen, Venlo, Netherlands), end repaired, A base added, and ligated to Illumina sequencing adapters. The ligation products were size selected (∼300 bp) by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina NovaSeq 6,000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

RNAseq raw reads were trimmed by fastp (v 0.18.0) (Chen et al., 2018) to remove adapter contamination and reads with high uncertainty (N > 10%) or low base quality with default parameters. The gene expression levels were quantified and differential expression analysis was conducted based on the reference genome of A. spinulosa (Huang et al., 2022), following methods described in previous studies (Sun et al., 2020). Briefly, index of the reference genome was built and paired-end clean reads were aligned to the reference genome using Hisat2 (v 2.2.1) (Kim et al., 2015) with parameters of “--sensitive --no-discordant --no-mixed -I 1 -X 1000”. The reads numbers matrix was generated using htseq (v 0.12.4) (Anders et al., 2015). We used fragments per kilobase of exon model per million reads mapped (FPKM) algorithm (Roberts et al., 2011) to obtain transcriptional profile. Differentially expressed genes (DEGs) of A. metteniana (H) vs. A. spinulosa (S) was identified using the DEseq2 R package (Love et al., 2014). Genes with FDR values ≤0.05 and FPKM values showing at least a 2-fold difference among samples were considered as DEGs. The hierarchical cluster was performed using R software with hclust function.

The Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were using clusterProfiler (Wu et al., 2021).

The Q1, Q3, retention time, declustering potential (DP) and collision energy (CE) were used for metabolite identification. The SCIEX OSV1.4 software was used to open the downtime mass spectrum file, and the chromatographic peaks were integrated and corrected. The chromatographic peaks were screened according to the minimum peak height of 500, signal-to-noise ratio of 5, smoothing number of 1 and other information. The peak area of each chromatographic peak represents the relative content of the corresponding substance. Finally, the integral data of all chromatographic peaks are derived to obtain the qualitative and quantitative results of metabolites (Bai et al., 2021). These metabolites were annotated using the KEGG database (Kanehisa and Goto, 2000), HMDB database (Wishart et al., 2013), and Lipidmaps database (Zhu et al., 2013), separately. PLS-DA analysis was applied to calculate the corresponding variable importance in projection (VIP) value. The condition of differential metabolites was VIP value of the PLS-DA model ≥1 and independent sample t-test’s p-value ≤0.05. Heatmap clustering analysis was performed in the R software with ComplexHeatmap package. KEGG pathway enrichment analysis of DAMs was examined using KOBAS (v2.0.12) (Xie et al., 2011). DAMs were considered to be significantly enriched in metabolic pathways when their p-values were <0.05.

We adopted a similar method as described in Yang, Sun et al. (Yang et al., 2023) to integrate metabolomics and transcriptomic analyses. The expression correlation between genes and metabolites was evaluated using ‘cor.test’ function under R software (v 3.5.1). We deemed a correlation significant if the absolute value of Pearson correlation coefficient ≥0.9 with a corresponding p-value ≤0.01.

As genomic data for A. metteniana is not available, we used the gene set of A. spinulosa as a reference to study R genes. Using a method similar to that described in the eggplant genome study (Li et al., 2021), we employed the RGAugury pipeline (Li et al., 2016) to screen the entire gene set for R gene prediction. The default p-value cutoff for initial R gene filtering was set to le-5 for BLASTP. Protein-protein interaction (PPI) network was analyzed using the Search Tool for the Retrieval of Interacting Gene (STRING) database (Szklarczyk et al., 2018), which included direct and indirect associations of proteins. The hub gene was determined using cytoHubba (Chin et al., 2014) with the MCC algorithm. The amino acid sequences of R genes were computed multiple sequence alignments using MAFFT (v 7.505) (Nakamura et al., 2018). Subsequently, the maximum-likelihood phylogenetic tree was obtained by using iqtree (v 2.0.6) (Minh et al., 2020) with parameters of ‘-bb 1,000 -pre iqtree -nt AUTO -m MFP -bnni’. The best model inferred by iqtree was “VT + R9”. The tree was displayed using Interactive Tree Of Life (Letunić and Bork, 2021).

Quantitative real-time PCR (qRT-PCR) has been emerged as an effective method to verify gene expression. We performed qRT-PCR experiment using a similar approach in the yellowhorn transcriptome (Zhao et al., 2021b). Total RNAs were reverse-transcribed using the PrimeScript first-strand cDNA synthesis kit (Takara, Dalian, China). Specific primers were designed using Primer-Blast tools (Ye et al., 2012). The qRT-PCR experiment was carried out using SYBRGreen Fast qPCR Master Mix (Sangon Biotech, China) on an ABI StepOne Plus Real-Time System (ABI, USA), following the manufacturer’s instructions. The quantitative PCR reaction conditions were as follows: 95°C for 90 s, followed by 95°C for 5 s, 60°C for 15 s, and 72°C for 20 s (45 cycles). Three biological replicates were included for each gene. Actin was used as the internal reference gene to normalize the qRT-PCR expression data. The 2^−ΔΔCT method (Livak and Schmittgen, 2001) was employed to calculate the relative mRNA abundance in each sample for every gene. Finally, we estimated the Pearson correlation coefficient of gene expression between the qRT-PCR and RNA-seq profiles in R software with “cor” function.

After collecting fresh leaves from A. spinulosa (S) and A. metteniana (H) (Figure 1A), the metabolic components were detected and analyzed. To improve the accuracy of analysis, each group included three replicates. We identified 373 and 399 metabolites in positive and negative modes, respectively (Supplementary Table S1). In A. spinulosa, the top 10 metabolites with highest metabolite abundance were “2-Caffeoyl-L-tartaric acid ", “Isocitric Acid”, “Kaempferol-3-O-glucoside-7-O-rhamnoside”, “Kaempferol-3-O-neohesperidoside”, “Skimmin ", “Luteolin-6-C-glucoside ", “Demethyl coniferin”, "γ-Linolenic Acid”, “Methylmalonic acid”, and “Aromadendrin-7-O-glucoside”. In A. metteniana, the top 10 metabolites with highest metabolite abundance were “Skimmin”, “Isocitric Acid”, “6-Methylmercaptopurine”, "γ-Linolenic Acid”, "α-Linolenic Acid”, “Spermine”, “Histidinol”, “N-Benzylmethylene isomethylamine”, “L-Glutamine” and “L-Lysine” (Figure 1B). The quantity of DAMs was statistically analyzed (Supplementary Table S2). Clustering heatmaps of DAMs showed significant differences in 111 metabolites. Of these, 82 DAMs were upregulated, and 29 metabolites were downregulated (Figure 1C). To verify the function of the DAMs, KEGG pathway enrichment was conducted. A total of 14 pathways were significantly enriched, including “valine, leucine and isoleucine degradation”, “pyruvate metabolism”, “2-Oxocarboxylic acid metabolism”, “valine, leucine and isoleucine biosynthesis”, “citrate cycle (TCA cycle)”, “biosynthesis of amino acids”, “biosynthesis of secondary metabolites”, “aminoacyl-tRNA biosynthesis”, “glucosinolate biosynthesis”, “glyoxylate and dicarboxylate metabolism”, “ABC transporters”, and so on (Figure 1D). Notably, the pathway with the most DAM enrichment pathways was the biosynthesis of secondary metabolites, which had 21 DAMs.

Six cDNA libraries were constructed from fresh leaves of A. spinulosa (S) and A. metteniana (H), yielding a total of 231.57 Mb of clean data, with 38.60 Mb obtained for each sample (Q30 ≥ 92.6%). The mRNA abundance of each gene in each sample was profiled using FPKM method. A total of 29,060 genes exhibited an FPKM expression value above 1.0 in at least one sample (Supplementary Figure S1). The hierarchical cluster dendrogram showed the expression pattern of biological repeats be clustered together (Supplementary Figure S2). Furthermore, a total of 14,639 DEGs were identified in the H vs. S comparison, with 7,764 upregulated genes and 6,875 downregulated genes (Figure 2A). Among the five most differentially upregulated genes, three genes had defined functions: indole-3-pyruvate monooxygenase YUCCA1 (YUCCA1), pleiotropic drug resistance protein 1 (PDR1), and cholesterol 22-hydroxylase CYP90B27 (CYP90B27). Correspondingly, all five of the most differentially downregulated genes had defined functions: 14 kDa zinc-binding protein (ZBP14), expansin-A10 (EXPA10), protein MEN-8 (MEN-8), methyl jasmonate esterase 1 (MJE1), and mannose-specific lectin (dfa). The results of DEGs KEGG enrichment analysis indicated that 38 pathways were enriched (Figure 2B; Supplementary Table S3). Furthermore, GO enrichment analysis revealed that the top 20 most significant GO categories as follows: “regulation of transcription, DNA-templated”, “protein phosphorylation”, “carbohydrate metabolic process”, “defense response”, “transmembrane transport”, “membrane”, “transmembrane transporter activity”, “copper ion binding”, “iron ion binding”, “ADP binding”, “sequence-specific DNA binding”, “ATPase activity”, “monooxygenase activity”, “protein kinase activity”, “microtubule binding”, “carbohydrate binding”, “protein serine/threonine kinase activity”, “DNA binding”, “DNA-binding transcription factor activity”, “oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular” (Figure 2C).

A total of 1,290 R genes were identified in the A. spinulosa genome using a genome-wide scanning pipeline (Li et al., 2016) (Supplementary Table S4). Among these R genes, 62.2% (803) belonged to the RLK category, while there were 110 NBS-related R genes, with 25 being of the TIR type. Further analysis revealed that 839 R genes were expressed in leaf tissue (Supplementary Figure S3). Of these, 606 R genes overlapped with DEGs, including 426 RLK-encoding R genes . Ten of the most differentially expressed genes were as follows: disease resistance protein L6 (L6), receptor-like protein kinase At1g49730 (At1g49730), disease resistance protein RUN1 (RUN1), leucine-rich repeat receptor-like serine/threonine-protein kinase BAM1 (BAM1), leucine-rich repeat receptor-like protein kinase At5g63930 (At5g63930), salt tolerance receptor-like cytoplasmic kinase 1 (OsI_16820), G-type lectin S-receptor-like serine/threonine-protein kinase SD2-5 (SD25), wall-associated receptor kinase 2 (WAK2), disease resistance protein RPS2 (RPS2) and L-type lectin-domain containing receptor kinase VIII.1 (LECRK81). Notably, of the differentially expressed 606 R genes, 551 genes could be further supported by known resistance genes reference from the latest PRGdb (Osuna-Cruz et al., 2017). These genes were significantly enriched in “Plant-pathogen interaction” and “MAPK signaling pathway - plant” pathway. Based on the GO enrichment analysis, these genes were assigned to 16 GO terms, such as “protein kinase activity,” “ATP binding,” “protein phosphorylation,” “defense response,” “protein binding,” and more (Supplementary Table S5). Using these 551 genes’ protein sequences as queries, 250 interactions were identified by using the STRING database (Supplementary Table S6). The first four key genes calculated by sytoHubba with the MCC algorithm were BAM1, inactive leucine-rich repeat receptor-like protein kinase At5g48380 (BIR1), receptor-like protein kinase BRI1-like 3 (BRL3), and leucine-rich repeat receptor-like serine/threonine/tyrosine-protein kinase SOBIR1 (SOBIR1) (Figure 3A). To investigate the phylogenetic relationship between R genes, we constructed an evolutionary tree. The maximum-likelihood phylogenetic tree showed that the closest relative genes of BAM1 were LRR receptor-like serine/threonine-protein kinase FLS2 (FLS2) and receptor-like kinase TMK (TMK) (Figure 3B). PPI network also revealed that BAM1 had interaction with FLS2 and TMK (Supplementary Table S6).

The expression levels of 14,097 DEGs were significantly correlated with the abundance of 111 DAMs. Further analysis showed that 585 significantly differentially expressed R genes were significantly positively correlated with 111 DAMs, while 574 significantly differentially expressed R genes were significantly negatively correlated with 111 DAMs (Supplementary Tables S7, S8). Subsequently, we constructed interaction networks between R genes and DAMs exhibiting significantly positive correlations, as well as for R genes and DAMs showing significantly negative correlations, respectively. According to the results of cytoHubba, in the significantly positively correlated network, the top 10 core genes are probably inactive leucine-rich repeat receptor-like protein kinase At5g48380 (BIR1), receptor-like protein kinase BRI1-like 3 (BRL3), BAM1, phytosulfokine receptor 1 (PSKR1), leucine-rich repeat protein kinase family protein (AT1G27190), leucine-rich receptor-like protein kinase family protein (At5g46330), SOBIR1, lysm-containing receptor-like kinase 1 (LYK1, also named as CERK1), leucine-rich repeat transmembrane protein kinase family protein (AT1G68400), and Leucine-rich repeat protein kinase family protein (AT1G67510). The first 10 core metabolites are N-Benzylmethylene isomethylamine, 1-O-p-Coumaroylquinic acid, Neochlorogenic acid (5-O-Caffeoylquinic acid), Ferulic acid-4-O-glucoside, Scopoletin-7-O-glucuronide, 1-O-Feruloylquinic acid, 3-O-Feruloylquinic acid, Luteolin-4′-O-glucoside, Kaempferol-7-O-glucoside, and Kaempferol-4′-O-glucoside. In the significantly negatively correlated network, the top 10 core genes are BIR1, BRL3, BAM1, At5g46330, AT1G27190, SOBIR1, Leucine-rich receptor-like protein kinase family protein (PSY1R), STRUBBELIG-receptor family 8 (SRF8), LRR receptor-like serine/threonine-protein kinase HSL2 (HSL2), and PSKR1. The first 10 core metabolites are Demethyl coniferin, Luteolin-8-C-glucoside (Orientin), Aromadendrin-7-O-glucoside, 4-O-(6′-O-Glucosylcaffeoyl)-3,4-dihydroxybenzoic acid, LysoPC 20:5, Procyanidin B2, Kaempferol-3-O-glucoside-7-O-rhamnoside, Kaempferol-3-O-neohesperidoside, Quercetin-7-O-rutinoside, and Cinnamtannin B1 (Figures 4A, B).

We utilized qRT-PCR to validate the gene expression profiles identified through Illumina sequencing analysis. Ten DEGs were randomly selected for validation, including wall-associated receptor kinase-like 1, leucine-rich repeat receptor-like serine/threonine-protein kinase BAM2, FLS2, and others (Supplementary Tables S9, S10). The gene expression of all 10 DEGs in the qRT-PCR experiment was consistent with the expression patterns in RNA-Seq results, with a Pearson correlation coefficient of the fold change between the qRT-PCR experiment and RNA-Seq being 0.78 (Supplementary Figure S4). These results confirmed the reliability of Illumina sequencing in this experiment.