Two diploid species and their allotetraploid hybrid, namely A. stenosperma V10309 (PI666100) ( Figure 1A ), A. valida GK30011 (PI468154) ( Figure 1B ), and [A. valida GK30011 × A. stenosperma V10309 (PI695393)]^4x ( Figure 1C ) were used in this study ( Additional File 1: Figure S1 ) (Chu et al., 2021; Gao et al., 2021; Chen et al., 2023). A. valida is a diploid species with the BB genome; A. stenosperma is a diploid species with the AA genome; and induced allotetraploid ValSten1 has AABB genome. Seeds of these genotypes were treated with two to three ml of a 0.5% solution of Florel^® Growth Regulator (Monterey Lawn and Garden, Fresno, California, USA) and incubated in a petri dish at 28°C for 18-24h to break seed dormancy. Seeds were sown in individual 4” pots with commercial potting mix Promix (Premier Horticulture Inc, Quakertown, PA, USA). The plants were kept in thrips-proof cages (47.5 cm³) (Megaview Science, Taichung, Taiwan) at 25-30°C, 80-90% RH, and a photoperiod of L14: D10 in the greenhouse. Seeds of the allotetraploid cultivar Georgia Green were pre-geminated in moistened paper towel and incubated in a growth chamber kept at 28°C for two to three days and used for thrips maintenance. One-to-two-week-old seedlings with one-to-two nodes and up to 16 leaflets of each genotype were used for TSWV transmission.

The hybrid ValSten1 plants were developed based on the protocol described in Gao et al. (2021). Briefly, in the greenhouse, A. valida plants were emasculated and pollinated with fresh pollen of A. stenosperma. Hybrid plants were identified by a series of pollen traits and tests as described in Gao et al. (2021). Once the hybrid plants were identified, whole genome duplication using small 20-cm lateral branch sections and colchicine was undertaken. Cuttings and resulting plants were then maintained in the greenhouse as stated in Gao et al. (2021). Pods harvested from these plants were assessed by cytological and phenotypic analysis. Three morphological variations viz., flower width, branch angle, and pod wieight variations further confirmed the induced allotetraploid status of ValSten1 plants.

Non-viruliferous thrips and viruliferous Frankliniella fusca thrips were maintained in separate growth chambers. Non-viruliferous thrips were maintained on leaflets of non-infected plants (cv. Georgia Green) within Petri dishes stuffed with a wet cotton round. Colonies were maintained by successive releases of ten adult female thrips, allowed to oviposit for 48h on a peanut leaflet dusted with a trace of pine pollen, and placed in growth chambers at 28-30°C and a photoperiod of L14: D10. Fresh leaflets and water were added to the Petri plates three times a week until emergence of the F₁ generation. TSWV viruliferous thrips colony was maintained similarly on TSWV-infected leaflets collected from the field in a separate growth chamber as described previously (Shrestha et al., 2013). During the off-season, viruliferous thrips were maintained on TSWV-infected leaflets generated by mechanical inoculation in the greenhouse (Marasigan et al., 2015; Shrestha et al., 2015).

TSWV viruliferous and non-viruliferous nature of thrips colonies was periodically tested by RT-qPCR using N-gene-specific primers as previously described with appropriate controls (Rotenberg et al., 2009; Shrestha et al., 2012; Shrestha et al., 2017). At each instance, a subset (~ten each) of viruliferous and non-viruliferous thrips were evaluated for TSWV infection status. All the viruliferous thrips evaluated tested positive and all the non-viruliferous thrips tested negative for TSWV. These indicated that the thrips colonies were true to their infection status or lack thereof.

F. fusca-mediated inoculation was conducted as per the established protocol previously (Shrestha et al., 2015). The experiment included two treatments: mock inoculation via non-viruliferous F. fusca thrips (non-infected) and TSWV inoculation via viruliferous thrips (TSWV-infected). Inoculated plants were maintained in thrips-proof cages (47.5 cm³) in the growth chamber at 27°C and ~80% humidity (Conviron, Pembina, ND, USA). After two weeks, the first fully expanded leaf of inoculated peanuts (ca 0.03 g) was tested by RT-qPCR following methods described previously (Shrestha et al., 2015; Chen et al., 2023) to assess TSWV-infection status.

Samples from plants two-to-three weeks post-inoculation were used. Five replications for each genotype were used. Leaflets were collected from the first fully expanded leaf below the terminal of each plant for RNA extraction. Total RNA was extracted by RNeasy plant mini kit following the manufacturer’s protocol (Qiagen, Valencia, CA, UGA). For each replicate, a leaflet sample was obtained from an individual plant. Thus, a total of 30 RNA samples were prepared for sequencing (three genotypes × two infection status × five replicates) and were stored at -80°C before shipping. Prior to library preparation, each sample’s integrity (RNA integrity number, RIN) was measured by using Aglient 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) for RNA quality control (QC). Two samples failed the QC test; therefore 28 samples were used for library preparation and sequencing.

The complementary DNA (cDNA) synthesis, cDNA libraries (messenger RNA library), and sequencing were undertaken by Novogene Corporation Inc. (Sacramento, CA, USA), as described in Catto et al. (2021). Illumina sequencing libraries were constructed using TruSeq RNA sample preparation kits. Briefly, mRNA was selected, fragmented, and first-strand cDNA was synthesized using random primers and reverse transcriptase. Subsequently, Polymerase I and RNase H were used to make the second-strand cDNA. An Illumina TruSeqLT adapter was ligated to the DNA fragments, and PCR amplification was performed for a minimal number of cycles with standard Illumina primers to produce the final cDNA libraries. Twenty-eight libraries were constructed and sequences using two lanes in the Illumina NovaSeq 6000 platform (pair-end 150 cycle sequencing setting, > 6GB raw data per sample).

In advance of the A. valida, A. stenosperma, and ValSten1 transcriptome assemblies ( Additional File 1: Figure S2 ), FastQC v0.11.9 and multiQC v1.11 were used to check the quality of raw reads before and after trimming (Andrews, 2010; Ewels et al., 2016). Trimmomatic v0.39 software was used with the default setting to remove adapters (Bolger et al., 2014). Also, Sortmerna v4.3.3 software was used with the SILVA database to remove rRNA contamination (Kopylova et al., 2012; Yilmaz et al., 2014; Glöckner et al., 2017). The rRNA decontaminated trimmed reads were converted from interleaved to paired files using BBMap v38.93 software for configuring files and for transcriptome assembly (Bushnell, 2014).

The rRNA decontaminated and trimmed reads from A. valida, A. stenosperma, and ValSten1 were used to generate respective de novo assemblies using Trinity v2.10.0 software with the default parameters (Grabherr et al., 2011). The sra2genes v4 software was used to clean up the assemblies using prior evidence from closely related species to address the possibility of over assembly of the transcriptome. Sra2genes is a complete pipeline to reconstruct genes from RNA data sources, and it includes several tools such as Cluster Database as High Identity of Tolerance (CD-HIT) v4.8.1, Exonerate v2.4.0, Blast+ 2.10.1, and A Genomic Mapping and Alignment Program for mRNA and expressed sequence tag (EST) Sequences – Genomic Short-read Nucleotide Alignment Program (GMAP-GSNAP) (Slater and Birney, 2005; Fu et al., 2012; Wu et al., 2016). CD-HIT v4.8.1 was used for the removal of potentially chimeric or misassembled transcripts from the input reads. Exonerate v2.4.0 was involved in the removal of all duplicated sequences. Blast+ 2.10.1 was used to separate the transcripts as various isoforms. GMAP-GSNAP was used to align the reads to the assemblies. Benchmarking Universal Single Copy Orthologs (BUSCO) v4.0.6 was used to determine assembly completeness before and after cleaning of the de novo assemblies against the Fabales odb10 lineage (n=5,366) (Simão et al., 2015; Seppey et al., 2019; Manni et al., 2021).

Trimmed reads were mapped to the respective de novo assemblies (see Data Availability for NCBI assessions) using Bowtie2 v2.4.1 with default mapping parameters (Langmead et al., 2009; Langmead and Salzberg, 2012; Langmead et al., 2019). Gene count estimates were derived from the mapped reads using RNA-Seq by Expectation Maximization (RSEM) v1.3.3 for A. stenosperma ( Additional File 2: Table S1 ), A. valida ( Additional File 2: Table S2 ), and ValSten1 ( Additional File 2: Table S3 ) (Li and Dewey, 2011). Custom R script was used to determine the fragments per kilobase million (FPKM) across all samples on R v4.1.0 using the following R libraries: dplyr, tidyverse, and stringr ( Additional File 1: Figures S3–S5 ) (R Core Team, 2021). DESeq2 was used to measure differentially expressed genes by comparing the gene counts from non-infected samples with virus-infected samples, where genes that had a |log₂ fold change (LFC)| ≥ 4 and a false discovery rate (FDR) < 0.05 were classified as being significantly differentially expressed (Love et al., 2014).

The de novo assemblies for A. stenosperma ( Additional File 3: Table S4 ), A. valida ( Additional File 3: Table S5 ), and ValSten1 ( Additional File 3: Table S6 ) were compared against an Arachis filtered subset of the NCBI database for non-redundant proteins (NR) and RefSeq genes using OmicsBox (Götz et al., 2008; Camacho et al., 2009). The OmicsBox tool also performed Blast2GO and Gene Ontology (GO) mapping to assign functional annotations to genes within each assembly (Conesa et al., 2005; Götz et al., 2008; Mi et al., 2019). Additional annotations were performed using InterProScan and the Kyoto Encyclopaedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000; Jones et al., 2014; Kanehisa et al., 2016). The GO terms were processed with topGO (https://www.bioconductor.org/packages/release/bioc/html/topGO.html) and visualized using rrvgo (https://bioconductor.org/packages/release/bioc/html/rrvgo.html) and the reduced + visualize Gene Ontology (REVIGO) web tool (Supek et al., 2011). GO terms down to level 3 were analysed.

DEGs from two wild peanut species: A. stenosperma and A. valida, their respective hybrid ValSten1, and previously published DEGs from two domestic peanut cultivars: A. hypogaea (SunOleic 97R) and A. hypogaea (Tifguard) (Catto et al., 2021) were used to determine DEG clusters using the online tool OrthoVenn2 (Xu et al., 2019). The parameters for DEG ortholog clustering in OrthoVenn2 were run with the cut-off value of 1e^-5. Overlapping regions were tested for significance using GeneOverlap (https://bioconductor.org/packages/release/bioc/vignettes/GeneOverlap/inst/doc/GeneOverlap.pdf).

Quantitative reverse transcription-polymerase chain reaction (RT-qPCR) was utilized to validate Arachis species transcripts following TSWV infection. Three sequences from each genotype with a |LFC| ≥4 and a false discovery rate (FDR) < 0.05 were randomly selected. The sequences were extracted with the tool seqtk. RT-qPCR was performed on plant samples obtained from four biological repeats from the remaining samples. Primers for targeted DEGs were designed by NCBI primer design (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). Primer sequences are listed in Additional File 1: Table S7 .

The cDNA was synthesized by a Go-Script reverse transcription system (Promega Corporation, Madison, WI) following the manufacturer’s protocol and then diluted 20-fold for quantitative polymerase chain reaction (qPCR). The reaction mix for qPCR included 2x GoTaq qPCR Master Mix, 1 μl of sequence-specific primers (final concentration of 250 mM), 2 μl cDNA of sample, and nuclease-free water for a final reaction volume of 20 μl. The reaction was run at 95°C for 2 min, followed by 40 cycles at 95°C for 15s, 58°C for 20s, and 72°C for 30s. The reaction was extended with a melting curve in a QuantStudio 3 System (applied biosystems by Thermo Fisher Scientific, Waltham, MA) to rule out non-specific binding. Two technical replicates for targeted transcripts and the reference gene (alcohol dehydrogenase class III) (Lai et al., 2021b), and water control were included in each RT-qPCR run. The log₂fold change of each target transcript in infected plants against mock-inoculated plants was calculated after normalization to the reference gene. The log₂ transformed (ratio of infected samples/ratio of non-infected samples) expression of target genes (transcripts) were correlated with Pearson’s correlation using the function “cor” in software R.

Total raw reads obtained from infected plants and non-infected plants of A. valida GK30011 (PI468154), A. stenosperma V10309 (PI666100), and ValSten1 were assembled de novo using Trinity platform. Total raw reads generated from the three genotypes were 222, 234, and 253 million pair reads, respectively, which after trimming amounted to 218, 231, and 250 million pair reads, respectively. The percentage of reads mapped to the de novo assembled transcriptome for A. stenosperma, A. valida, and ValSten1 genotypes were 86%, 87%, and 80%, respectively ( Additional File 1: Table S8 ). These reads were assembled into 141,144 (A. valida), 106,374 (A. stenosperma), and 137,039 (ValSten1) contigs. The assembly of A. stenosperma contained 4,571 (85%) complete BUSCOs, which included 2,571 (48%) single-copy and 2,000 (37%) duplicated orthologs. Similarly, A. valida contained 4,724 (88%) complete BUSCOs, which included 2,545 (47%) single-copy and 2,179 (41%) duplicated orthologs. For ValSten1, there were 4,670 (87%) complete BUSCOs, which included 2,209 (41%) single-copy and 2,461 (46%) duplicated orthologs. One infected sample of A. stenosperma showed low RIN (RNA integrity number) and one non-infected sample of A. valida that showed uneven baseline at QC were not processed from the initial 30 libraries.

The reads obtained from infected and non-infected samples from the three genotypes were normalized and clustered using FPKM and principal component analysis (PCA) for comparison. The PCA clustered TSWV infected samples of the three genotypes separately from the non-infected ones ( Figure 2 ). However, one sample (asten_paired_V3B) in A. stenosperma was removed due to the unexpected clustering in PCA, although it did not have a reduced FPKM value ( Additional File 1: Figure S6 ). Additional checks on infection status were performed by mapping reads, using RSEM and Bowtie2, from A. stenosperma and A. valida to the ValSten1 de novo assembly and clustering the samples via PCA ( Additional File 1: Figure S7 ). Differentially expressed genes (DEGs) observed for A. stenosperma, A. valida, and ValSten1 in response to TSWV were 3,196 (596 overexpressed and 2,627 underexpressed; Figure 3A ), 8,380 (6,332 overexpressed and 2,048 underexpressed; Figure 3B ), and 1,312 (633 overexpressed and 679 underexpressed; Figure 3C ), respectively. TSWV-infected samples of A. valida had more DEGs (8,380) compared with A. stenosperma (3,196) and ValSten1 (1,312). A higher percentage of DEGs for A. stenosperma were underexpressed, whereas more overexpressed genes were identified in A. valida. Similar numbers of underexpressed and overexpressed genes were found within ValSten1.

DEGs observed in the wild species in response to TSWV infection were functionally annotated. The de novo assemblies included 107,043 transcripts, 149,877 transcripts, and 138,389 transcripts (non-significant and significant genes) of A. valida, A. stenosperma, and ValSten1, respectively.

Gene ontology (GO) provided context for the functionality of genes and comprised three level 1 categories: biological process (BP), cellular component (CC), and molecular function (MF). GO terms within the BP category provided biological relevance by attributing biological objectives to gene products. Significantly enriched GO terms were determined by the Revigo tool (Supek et al., 2011) by comparing the GO terms distribution from DEGs to that of the entire transcriptome, also referred to as the background genes. DEG specific GO terms that were overrepresented were considered significantly enriched (p < 0.05) with respect to the background.

In A. stenosperma, 127 BP GO terms were significantly enriched among DEGs across all GO term levels ( Additional File 1: Figures S8A, B ; Additional File 4: Tables S9, S10 ), with 14 terms being classified as levels 2 & 3 ( Figures 4A, B ). In A. valida, 256 BP GO terms were significantly enriched among DEGs across all GO term levels ( Additional File 1: Figures S8C, D ; Additional File 4: Tables S11, S12 ), with 19 terms being classified as levels 2 & 3 ( Figures 4C, D ). In ValSten1, 135 BP GO terms were significantly enriched among DEGs across all GO term levels ( Additional File 1: Figures S8E, F ; Additional File 4: Tables S13 , S14 ), with 9 terms being classified as levels 2 & 3 ( Figures 4E, F ).

To determine the transcriptional changes in each genotype related to TSWV infection, the number of orthologous clusters between A. stenosperma, A. valida, ValSten1, A. hypogaea (SunOleci 97R), and A. hypogaea (Tifguard) (Catto et al., 2021) were compared using the OrthoVenn2 web platform ( Figure 5 ). Orthologous clustering analysis resulted in 3,965 clusters of DEGs that were commonly shared by at least two genotypes ( Additional File 5: Table S15 ) and 15 single-copy DEG clusters from all five genotypes ( Additional File 5: Table S16 ). In total, 71 DEG clusters were found to contain DEGs shared between all five genotypes, with cluster53, cluster179, and cluster185 relating to the putative disease resistance protein RGA3 (Song et al., 2003; Van Der Vossen et al., 2003) (UniProt ID: Q7XA40) and defense response (GO:0006952; Additional File 5: Table S15 ). There were 17 DEG clusters that comprised four of the genotypes, but not in the susceptible A. hypogaea (SunOleic 97R), with cluster674 relating to the TMV resistance protein N (UniProt ID: Q40392) (Whitham et al., 1994; Dinesh-Kumar and Baker, 2000; Dinesh-Kumar et al., 2000; Caplan et al., 2008) and signal transduction (GO:0007165; Additional File 5: Table S15 ). The highest overexpressed gene, with a LFC of 29.6, was found in A. stenosperma and was annotated as linoleate 9S-lipoxygenase (ArasteEVm001500t4). Manual assessment determined that such a large LFC was caused by lack of mapped reads (no detectable expression) in the mock inoculated/non-infected samples. This gene was found to be in cluster4, containing genes from all genotypes, and was functionally annotated as linoleate 9S-lipoxygenase (P38414) (Hilbers et al., 1994) and oxylipin biosynthetic process (GO:0031408; Additional File 5: Table S15 ).

With respect to A. stenosperma, A. valida, and ValSten1, orthologous DEG clustering analysis resulted in 1,507 DEG clusters that were commonly shared by at least two genotypes: A. stenosperma ⋂ A. valida (779), A. stenosperma ⋂ ValSten1 (79), A. valida ⋂ ValSten1 (412), or A. stenosperma ⋂ A. valida ⋂ ValSten1 (237) ( Additional File 1: Figure S9 ). Additionally, 1,574 DEG clusters were found to be specific to A. stenosperma (269), A. valida (1,230), and ValSten1 (75) ( Additional File 1: Figure S9 ). Sixty nine of the 237 orthologous clusters shared by the three wild peanut genotypes were reported as containing single-copy DEGs ( Additional File 5: Table S17 ). All pairwise comparisons of DEG clusters from A. stenosperma ⋃ A. valida (2,404), A. stenosperma ⋃ ValSten1 (1,439), A. valida ⋃ ValSten1 (2,357) showed more overlap than expected by chance (Fisher’s exact test) p=1.1e^-50, p=3.5e^-12, and p=7.4e^-47, respectively.

The phytovirus response DEGs from A. stenosperma (3,196), A. valida (8,380), and ValSten1 (1,312) were grouped into three major categories: defense, phytohormone, and photosynthesis related genes ( Table 1 ). The categories were chosen based on the study with resistant and susceptible cultivated peanuts (Catto et al., 2021). Within the defense related DEGs, the percentage (No. of overexpressed DEGs out of total DEGs within category) in A. stenosperma, A. valida, and ValSten1 were 34% (25/73), 64% (490/763), and 55% (69/126), respectively ( Table 1 ). A similar pattern was observed in the case of phytohormone related DEGs. Upregulation of phytohormone related DEGs of the eight examined categories was higher in A. valida. The percentages (No. of overexpressed DEGs out of total DEGs within category) in A. stenosperma, A. valida, and ValSten1 were 9% (1/11), 51% (100/198), and 36% (16/44), respectively ( Table 1 ). Regarding photosynthesis related DEGs, the percentages (No. of overexpressed DEGs out of total DEGs within category) in A. stenosperma, A. valida, and ValSten1 were 10% (3/29), 46% (249/536), and 38% (28/73), respectively.

Three DEGs from each genotype were randomly selected and their expression values were validated using RT-qPCR ( Additional File 1: Table S7 ). A positive correlation was found between the expression from both RNASeq and RT-qPCR across all three genotypes (cor=0.87, t=4.6, df=7, p=0.002; Additional File 1: Figure S10 ).