Thrips samples were collected from 16 distinct sites on Hainan Island, encompassing various crops, including cucurbits, with sampling conducted between July and December. Total RNA was extracted separately from each thrips population at these sites, with each population containing a minimum of 50 individuals. When multiple thrips species were present at a site, individuals were pooled by species for RNA extraction and analysis. RNA extraction was carried out using the Trizol method, and RNA quantity was assessed with a NanoDrop 2000 spectrophotometer.

To ensure data representativeness and scientific rigor, RNA from the 16 thrips populations was pooled into six composite RNA samples based on thrips species and host plant. This was due to limited RNA, high costs, and our broad research aim. These pooled samples were prepared for metatranscriptomic sequencing. Total RNA samples underwent ribosomal RNA removal using the Ribo-Zero Gold kit (specific for human/mouse/rat rRNA), and paired-end 150 bp sequencing was performed on the Illumina NovaSeq^™ 6000 platform. Library preparation and sequencing were completed by Lianchuan Biotechnology Co., Ltd. (Hangzhou, China), yielding 81 GB of raw sequencing data. This sequencing depth provides robust data for the identification and comprehensive analysis of thrips virome diversity.

After assessing the quality of raw sequencing data from each library, data trimming was performed using Trimmomatic v0.32 (Bolger et al., 2014) to remove low-quality reads and adapter sequences, resulting in a set of clean reads. To efficiently detect viral fragments, these clean reads were initially aligned to the host genome (F. intonsa reference genome) using Bowtie2 v2.4.4 (Langmead and Salzberg, 2012), allowing for the removal of host-derived sequences. The remaining, non-aligned reads were then assembled de novo using Trinity v2.8.5, which generated preliminary viral contigs.

The assembled contigs were annotated using the blastX algorithm in DIAMOND v3.12 with an E-value threshold of 1 × 10⁻⁵, balancing sensitivity and specificity to ensure reliable viral identification while minimizing false positives. The selected viral contigs were subsequently matched against the non-redundant (NR) protein database for initial virus identification. Open reading frames (ORFs) within the viral contigs were identified using NCBI ORFfinder,¹ and homologous proteins were identified using HMMER v3.3 against the Pfam protein family database and the RNA-dependent RNA polymerase (RdRP) database of RNA viruses. This approach allowed for further classification and identification of novel viruses based on the latest ICTV species classification standards.

Finally, to estimate viral abundance across the libraries, the quality-controlled paired-end reads were mapped to the identified viral genomes using Bowtie2. Alignments were processed using Samtools v0.1.1 (Li et al., 2009), and RPKM and FPKM values for each contig were calculated using Salmon, providing quantitative insights into viral prevalence within the libraries.

Primers specific to each identified virus were designed for Sanger sequencing. Total RNA from the six libraries obtained in section 2.2.1 was subjected to RT-PCR to synthesize cDNA, with an expected product size of 400 to 600 bp and an annealing temperature (Tm) of 55–60°C. To address potential challenges in primer design, such as complex secondary structures or low copy numbers, four primer pairs were designed for distinct regions of each viral sequence. Each primer pair was tested sequentially to achieve the desired bands; if none of the four pairs yielded a band, the virus was considered absent in that sample. Primer design was conducted using Primer-BLAST,² and primer sequences are provided in Supplementary Table 2. Detailed reaction conditions are as follows: The Premix Taq PCR reaction system consisted of the following components with their respective volumes. Premix Taq was 12.5 μL, cDNA was 1 μL, F-primer was 1 μL, R-primer was 1 μL, and ddH₂O was added to make the total volume up to 25 μL.

The PCR conditions included an initial denaturation at 95°C for 5 min, followed by 35 cycles of 30 s at 95°C, 40 s at 55°C, and 40 s at 72°C, concluding with a final extension at 72°C for 10 min. PCR products were analyzed using 1.2% agarose gel electrophoresis. If the target band was observed, it was excised from the gel, and the cDNA was purified. Purified products were sent to BGI-Shenzhen for sequencing. Sequencing results were aligned using Geneious v9.0.2, and sequences with similarity above 95% were considered indicative of the virus’s presence.

Phylogenetic analysis was performed by aligning conserved viral sequences identified in this study with known sequences from the same viral families using MAFFT v7.407 with the E-INS-I algorithm. Special focus was placed on RNA-dependent RNA polymerase (RdRp) sequences of RNA viruses, as these are key conserved regions useful for differentiating RNA virus types. Unaligned sequences and poorly aligned regions were removed with TrimAl v1.4 to enhance phylogenetic accuracy. The phylogenetic tree was then constructed using the maximum likelihood (ML) method in IQ-TREE v2.0.3, with bootstrap analysis conducted 1,000 times to validate the reliability of the results. The amino acid substitution model was selected based on the best-fit analysis provided by ModelTest.

Non-viral reads were classified and annotated using Kraken v2.0.7-beta 47 (Wood and Salzberg, 2014). For accurate microbial identification, rRNA reads were filtered from the dataset using SortmeRNA (Kopylova et al., 2012). The filtered rRNA reads were then compared to the nucleotide database using the NCBI BLAST suite, applying an E-value threshold of 10⁻³. The microorganism with the lowest E-value match was identified as the closest species.

Between July 2021 and December 2022, we conducted a survey on the occurrence of viral diseases in cucurbit and other vegetable fields in Hainan, China (Supplementary Table 1). During this period, both thrips and plant samples were collected simultaneously, allowing for a direct comparison of virus presence in thrips populations with the incidence of viral symptoms in corresponding plant samples. A total of 543 leaf samples exhibiting symptoms such as dwarfism, mosaic patterns, malformation, yellow spots, or necrosis were collected and stored at −80°C for further analysis. DNA and RNA were extracted using a plasmid DNA extraction kit (Tiangen Biotech, Beijing, China) and the Trizol method, respectively.

Four target plant viruses (WSMoV, MYSV, WGMMV, TeMV), previously identified in the thrips virome, were detected in these plant tissue samples using reverse transcription and PCR. PCR products were analyzed via agarose gel electrophoresis, with target DNA bands excised and purified for sequencing. Detected plant viruses were then identified through sequence alignment with the NCBI nucleotide database. Detection rates for suspected viral samples were grouped by plant family to evaluate the prevalence of target viruses in cucurbit and other vegetable crops. These four viruses were selected as they were the only plant viruses identified in our initial virome analysis of thrips.

We conducted a Pearson correlation analysis to examine the detection of four plant viruses within the T. palmi and F. intonsa viromes. Specifically, we assessed the correlation between the detection rate of MYSV in cucurbit crops and the occurrence of T. palmi. Additionally, we examined the correlation between the detection rate of WSMoV in both Cucurbitaceae and Fabaceae and the abundance of T. palmi. In both analyses, virus detection in plant samples was based on presence/absence rather than incidence rates. Similarly, we analyzed the correlation between the detection rate of WGMMV in cucurbit crops and the occurrence of F. intonsa, as well as the correlation between the detection rate of TeMV in leguminous crops and the occurrence of F. intonsa.

From each of the six thrips sample libraries, a total of 12.88 to 13.80 GB of data were generated, reflecting the data range across different species. After filtering out low-quality reads, 65,532,884 clean reads were retained. Through de novo assembly, reference genome alignment, and BLASTx searches, we identified 36,269 initial viral contigs, with lengths ranging from 200 to 20,390 nt and an average length of 502 nt. Following iterative extension using Bowtie2 and Samtools, and after rigorous screening and quality control, we obtained 30 high-quality viral sequences from the six sample pools, identifying 19 distinct viruses. These viruses were categorized into nine viral families, with three unclassified viral groups comprising one double-stranded RNA virus (Totiviridae), 11 positive-sense single-stranded RNA viruses, six negative-sense single-stranded RNA viruses, and one unclassified RNA virus (Supplementary Table 3). Among the identified viruses, negative-sense single-stranded RNA viruses were the most prevalent. Picornavirales accounted for 35.63% of total viral reads, followed by Amarillovirales (27.80%) and Martellivirales (16.06%) (Table 1).

Based on host categories, the viruses carried by thrips were classified as fungal viruses (one species), plant viruses (four species), arthropod viruses (13 species), and vertebrate viruses (one species) (Supplementary Table 3). The four plant viruses identified were MYSV, WSMoV, TeMV, and WGMMV, with MYSV and WSMoV detected in T. palmi, and TeMV and WGMMV in F. intonsa (Figure 2A).

Analysis of the viral community composition in the three thrips samples revealed significant differences in virus species and abundance among samples (p < 0.05). In the T. palmi sample (HL23), 11 viruses were detected, with Thrips iflavi-like virus 10 being the most abundant, accounting for 25.26% of the total. In the F. intonsa samples (DH12, HH22, QH32, C42), between 2 and 6 viruses were detected, each dominated by a different virus: Thrips flavilike viruses 5 and 6 in DH12, Thrips picorna-like virus 20 in HH22, Thrips virga-like virus 1 in QH32, and SY-Telosma mosaic virus in C42. In the M. usitatus sample (DD11), four viruses were identified, with Thrips iflavi-like virus 4 being the most prevalent, comprising 66.84% of the total (Figure 2A).

Principal component analysis (PCA) demonstrated that viral community compositions differed significantly among the three thrips species (p < 0.05), while viral composition was similar among thrips on the same host plant family. This finding suggests that viral communities are primarily influenced by the host plant (Figures 2B,C).

In three F. intonsa sample libraries (C42, DH12, HH22) and one T. palmi sample library (HL23), we identified a novel double-stranded RNA (dsRNA) virus in the Totiviridae family, tentatively named Thrips toti-like virus 3 (TTV3). The full genome of TTV3 is 4,858 nucleotides long and contains two open reading frames (ORFs). This specific ORF arrangement in the TTV3 genome is significant as it provides the basis for further exploring the virus’s replication mechanism and its interaction with the host thrips. The first ORF encodes a protein of unknown function with a molecular mass of 81.072 kDa, while the second ORF encodes an RNA-dependent RNA polymerase (RdRp) (Figure 3A). Phylogenetic analysis shows that TTV3 forms a distinct evolutionary branch within the Totivirus genus, exhibiting 36.83% amino acid sequence similarity to Camponotus nipponicus virus, a virus found in tree-dwelling ants (Figure 3B).

Twelve positive-sense single-stranded RNA (+ssRNA) viruses were identified, representing various families, including Virgaviridae, Potyviridae, Flaviviridae, Picornavirales, and Hepeviridae. In one F. intonsa sample (HH22), a sequence closely matching Watermelon green mottle mosaic virus (WGMMV) was detected, with 99% amino acid homology (Figure 3C). Additionally, a novel virus related to Virgaviridae, identified in multiple thrips samples (DH12, HH22, QH32, C42, HL23), was named Thrips virga-like virus 1 (TVV1). The TVV1 genome is approximately 11,938 nucleotides in length and contains three open reading frames (ORFs), one of which includes both Hel and RdRp domains. TVV1 shares 42.8% amino acid similarity with F. intonsa-associated virga-like virus 2 (FAVV2) and forms a distinct evolutionary branch with FAVV2 (Figure 3C).

In two F. intonsa samples (C42 and HH22), sequences of Telopea mosaic virus (TeMV) from the Potyvirus genus were identified and designated SY-TeMV. The SY-TeMV genome is 10,057 nucleotides long, with a 9,630-nucleotide ORF encoding 3,209 amino acids. Phylogenetic analysis reveals that SY-TeMV forms a monophyletic branch with other Chinese isolates (Figure 3D).

Furthermore, two novel Flavivirus sequences were isolated from F. intonsa and M. usitatus samples, designated Thrips flavi-like virus 5 (TFV5) and Thrips flavi-like virus 6 (TFV6), respectively. The genomes of these viruses are notably longer than those of other members in the Flaviviridae family (Figure 3C). Phylogenetic analysis suggests that TFV5 and TFV6 are distantly related to other insect-specific flaviviruses (Figure 3E).

In the Picornavirales order, Thrips iflavi-like virus 4 (TIV4) and Thrips picorna-like virus 20 (TPV20) were detected in M. usitatus and F. intonsa samples, respectively. Additionally, four novel viruses were isolated from the T. palmi sample, including Thrips iflavi-like virus 10 (TIV10), Thrips iflavi-like virus 11 (TIV11), and Thrips picorna-like virus 12 (TPV12). These viruses belong to the Iflaviridae family but occupy distinct evolutionary branches (Figure 3F).

Lastly, a novel Hepeviridae virus, named Thrips hepe-like virus 16 (THV16), was identified in the T. palmi sample. The THV16 genome is 5,802 nucleotides long and contains three ORFs. Phylogenetic analysis indicates that THV16 is distantly related to other hepatitis viruses (Figure 3G).

Four novel virus sequences, suspected to belong to the negative-sense single-stranded RNA (−ssRNA) viruses, were identified. Among these, one virus was detected in a T. palmi sample (HL23) and was named Thrips mononega-like virus 17 (TMV17), potentially a new member of the Lispiviridae family. The other three viruses were detected in M. usitatus (DD11) and T. palmi (HL23) samples, showing amino acid sequence homology with viruses from the Rhabdoviridae family, ranging from 42.36 to 52.62%. These viruses were named Thrips rhabdo-like virus 7 (TRV7), Thrips rhabdo-like virus 14 (TRV14), and Thrips mononega-like virus 15 (TMV15), respectively. Due to limited sequencing depth, their full genome structures remain unresolved, and they were identified as vOTUs (virus operational taxonomic units) based solely on the conserved sequence of RNA-dependent RNA polymerase (RdRp).

Phylogenetic analysis revealed that TMV17 is most closely related to Lishi Spider Virus 2 and F. intonsa associated mononegavirales virus 1, with nucleotide and amino acid sequence similarities of 29.71 and 29.73%, respectively. TRV7, TRV14, and TMV15 clustered with other arthropod-hosted Rhabdoviridae members, distinct from viruses infecting vertebrate and plant hosts. Notably, TRV14 is most closely related to Soybean Thrips rhabdo-like virus 3, which infects soybean thrips, with 49.42% amino acid similarity. TMV15 shows 41.12% amino acid similarity to F. intonsa associated mononegavirales virus 3, which infects F. intonsa, suggesting that these thrips-associated viruses may represent a new evolutionary branch (Figure 3H).

Additionally, two full-length sequences with high similarity to the LRNA segment of Orthotospovirus were obtained from the T. palmi sample (HL23). These sequences were most closely related to WSMoV and MYSV, with nucleotide and amino acid sequence similarities of 95.06 and 96.42%, respectively (Figure 3I).