Discovery of novel non-retroviral endogenous viral elements reveals their long-term co-evolution with spiders

doi:10.21203/rs.3.rs-4996788/v1

Download PDF

Article

Discovery of novel non-retroviral endogenous viral elements reveals their long-term co-evolution with spiders

https://doi.org/10.21203/rs.3.rs-4996788/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Endogenous viral elements (EVEs) are widespread in the genomes of various organisms and have played a crucial role in evolution. Historically, research on EVEs primarily focused on those derived from retroviruses; however, the significance of non-retroviral EVEs (nrEVEs) has gradually gained recognition. In this study, we employed a novel approach that combines protein structure prediction with sequence analysis to identify a large group of previously unrecognized nrEVEs across spider genomes. Additionally, we identified nrEVE-related messenger RNAs, microRNAs, and PIWI-interacting RNAs in spiders, suggesting that these nrEVEs may be functionally active. We also experimentally confirmed the presence of spider nrEVEs and their transcripts in individual spiders. Evolutionary analysis suggests that these spider nrEVEs originated from an ancestral nuclear arthropod large DNA virus (NALDV) belonging to the order Lefavirales, class Naldaviricetes, approximately 270 million years ago. The integration of these nrEVEs occurred prior to the last common ancestor of the Araneae, indicating a long-term co-evolutionary relationship between these nrEVEs and spiders. Our findings reveal a novel group of nrEVEs and provide valuable insights into their evolutionary relationship with arthropods.

Biological sciences/Evolution/Coevolution

Biological sciences/Microbiology/Virology/Viral evolution

The genomic integration of viral elements has been reported across all three domains of life^1,2. The genomes of organisms commonly harbor transposons and endogenous viral elements (EVEs), and retroviruses have historically been recognized as the major promoters of EVE integration³. Recent advances in genome sequencing and bioinformatics have revealed numerous non-retroviral endogenous viral elements (nrEVEs) in various organisms^4,5. The mechanisms of nrEVEs integration remain largely unclear; however, nrEVEs evidently undergo selection by the host, who sometimes exploit them for various functions^6,7. These apparent remnants of ancient viral infections, which serve as molecular fossils, provide valuable insights into the coevolution dynamics between hosts and viruses over evolutionary timescales^2,8,9.

Arthropods comprise the most diverse group of animals on Earth and include two major subphyla: Mandibulata (e.g., insects and crustaceans) and Chelicerata (e.g., horseshoe crabs, scorpions, and spiders)¹⁰. Arthropods host a diverse and complex array of viruses owing to their rich species diversity. Among these, a distinctive group of viruses, known as nuclear arthropod large DNA viruses (NALDVs), are exclusively distributed in insects and crustaceans. NALDVs, now classified within the class Naldaviricetes, possess circular double-stranded DNA (dsDNA) genomes ranging from 80 to 500 kb, and replicate within the host cell nucleus. They encompass four identified families: Baculoviridae, Hytrosaviridae, Nimaviridae, and Nudiviridae¹¹, and recently proposed Filamentoviridae¹², as well as the unclassified Apis mellifera filamentous virus (AmFV). Hundreds of NALDV species have been isolated, forming a monophyletic group separate from other arthropod large DNA viruses classified as nucleocytoplasmic large DNA viruses (NCLDVs) within the phylum Nucleocytoviricota¹². The virion envelope protein known as peroral infectivity factors (PIFs), unique core proteins identified solely in Naldaviricetes and bracoviruses of Polydnaviriformidae, are shared across all NALDVs¹¹. Within the class Naldaviricetes, Baculoviridae, Hytrosaviridae, Nudiviridae, and Filamentoviridae are further grouped to collectively form the order Lefavirales as they all share the so-called late expression factors (LEFs)¹³.

nrEVEs originating from NALDVs have been identified in several arthropods and represent a rare integration of large DNA viral genome segments into eukaryotic genomes. Genomic integration of nudiviruses into cultured cells has been reported and fossilized remnants of nudiviruses have been detected in the genomes of numerous arthropods^14,15. Bracoviruses are derived from an endogenous nudivirus within the wasp genome¹⁶ and exhibit remarkable symbiotic relationships with wasp hosts. Bracoviruses are essential for wasps, as they can suppress the host immune response and create a suitable environment for wasp larval development¹⁷. Endogenous nimaviral genomes were discovered in host crustacean genome¹⁸, and endogenization and domestication of filamentoviruses and AmFV were identified in multiple hymenopterans¹⁹. In summary, the nrEVEs of NALDVs appear to be common feature.

Spiders, classified as chelicerates (subphylum Chelicerata), diverged from common ancestors of insects and crustaceans early in evolutionary history²⁰. Currently, more than 50,000 documented spiders (order Araneae) are members of two suborders, Mesothelae and Opisthothelae; the latter contains two infraorders, Mygalomorphae and Araneomorphae²¹. Members of the suborder Mesothelae form a sister group to all other living spiders, as they retained ancestral characters, such as a segmented abdomen with spinnerets in the middle and two pairs of book lungs. Spiders in the suborder Mygalomorphae, including tarantulas and trapdoor spiders, are robust and “hairy,” with large chelicerae and fangs. The Araneomorphae lineage, also called “true spiders,” comprises the vast majority (~ 93%) of living spiders, including orb-web spiders, wolf spiders, and jumping spiders²¹. Research on spiders has predominantly focused on silk and venom²². Only a limited number of RNA viruses and EVEs have been identified in spiders²³, and their impact on spider populations and ecosystems remains unknown.

Conventional methods of EVEs discovery rely primarily on homologous sequence searches for known viruses within the host species. Hence, EVEs with low similarity to known viral sequences are unlikely to have been discovered. Recent, advancements in artificial intelligence technology, particularly in protein structure prediction and comparison, have offered new avenues for the discovery and classification of distantly related viral sequences and previously unknown protein components^24–26. In this study, we combined protein structure prediction and sequence alignment to identify baculovirus-related sequences in databases and uncover novel nrEVEs in spider genomes. We further investigated the distribution of these nrEVEs in spider genomes and uncovered the messenger RNAs (mRNAs), PIWI-interacting RNAs (piRNAs), and microRNAs (miRNAs) related to these nrEVEs in the available transcriptional data on spiders. We also experimentally verified the presence of nrEVEs and their transcripts in the spider samples. Finally, we reconstructed the evolutionary relationships between spider nrEVEs, NALDVs, and their host arthropods.

Discovery of abundant baculoviral homologs in spider genomes using a combined strategy of structure and sequence comparisons

We developed a comprehensive approach (Fig. 1a) to identify the uncovered NALDVs nrEVEs. First, by querying the predicted baculoviral protein structures in the AlphaFold Protein Structure Database cluster (AFDB50), we identified a distinct spider species, Araneus ventricosus, harboring multiple baculoviral homologs dispersed across 283 scaffolds of the genome (Table S1). Structural alignment revealed significant similarities between the proteins identified in A. ventricosus and their baculoviral homologs. For example, three proteins showed a homology probability over 0.95, with the conserved baculovirus core proteins (PIF-1, PIF-4, and LEF-4) of Autographa californica multiple nucleopolyhedrovirus (AcMNPV), while their amino acid sequence identities were below 20% (Fig. 1b). Subsequently, the open reading frames (ORFs) of selected scaffolds (Table S2) with large DNA virus-coding strategies (e.g., higher ORF density) were extracted, and their structures were predicted. After structural comparison and sequence alignment of these proteins with the predicted baculoviral protein structures, we identified over 30 homologs of baculoviruses, including 22 baculoviral core proteins (conserved among all baculoviruses) (Fig. 1a). These core proteins included those involved in DNA replication and transcription (e.g., helicase, dnap, p47, lef-4, lef-5, lef-8, and lef-9), as well as virion structure and other functions (e.g., pif-0, pif-1, pif-2, pif-3, pif-4, pf-5, pif-6, pif-8, vp39, 38k, 49k, odv-ec43, odv-ec27, p33, and p48) (Table S2).

To comprehensively investigate the occurrence and distribution of baculoviral homologs in spider genomes, we screened 22 baculoviral core homologs of A. ventricosus against all available arachnid genome assemblies in the National Center for Biotechnology Information (NCBI) Whole Genome Shotgun (WGS) database, encompassing 92 species from diverse orders, including Araneae, Trombidiformes, Sarcoptiformes, Ixodida, Mesostigmata, Scorpiones, Opiliones, and Pseudoscorpiones (Table S3). Our analysis revealed 10641 baculoviral homologs across 89 species (Table S4). Among the 32 assembled spider genomes, 10091 baculoviral homologs were identified in 1742 scaffolds, spanning 13 distinct families (Fig. 1c), with 12 families from the infraorder Araneomorphae and one family (Theraphosidae) from more ancient infraorder Mygalomorphae. Most baculoviral homologs were identified in the order Araneae, whereas a few baculoviral homologs with lower identity levels were identified in other orders (Fig. 1d).

Broad distribution of baculovirus-related nrEVEs in spider genomes

Spider scaffolds harboring baculoviral homologs (Table S5) displayed a wide spectrum, ranging from several kilobase pairs to hundreds of millions of base pairs (Fig. 2a). Several species containing hundreds of baculoviral homologs were further analyzed (Fig. 2b). All chromosomes of T. clavata from SpiderDB²⁷ (Table S6), and selected chromosomes of Hylyphantes graminicola, Latrodectus elegans, Dolomedes plantarius, Meta bourneti, Metellina segmentata, and Parasteatoda lunata from the WGS database contain clustered distribution of baculoviral homologs. Chromosomal regions with baculoviral homologs usually contained direct repeats, inverted repeats, large fragments of repetitive sequences, and host genes, confirming that these baculoviral homologs were endogenous in the spider genome (Fig. 2c). Compared to the adjacent regions, the baculoviral homology rich regions displayed a notable increase in the percentage of ORFs coverage, suggesting sequential integration of large viral DNA fragments (Fig. 2c). Collectively, the distribution and localization of baculoviral homologs confirmed the widespread existence of baculovirus-related nrEVEs in spider genomes.

Spider nrEVEs-related mRNAs, piRNAs and miRNAs exist in spiders

We further conducted an extensive survey of the baculoviral homologs in the NCBI Transcriptome Shotgun Assembly (TSA) database. In total, 1689 TSA datasets covered 1138 different spider species from 76 families (Table S7). Remarkably, mRNAs with baculoviral homologs were widespread in Araneae (Table S8). In addition, the number of baculoviral homologs in different spider families was strongly correlated with the number of datasets available for each family (Fig. 3a). However, the detection rate of core baculoviral homologs varied among the different families, and no obvious correlation was identified between the detection ratio and spider classification (Fig. S1a).

We further examined the transcriptomic data of T. clavata available in SpiderDB, which contains the most detailed data for different spider tissues. Nearly full-length mRNA of certain baculoviral homologs were recovered (Fig. S1b). Although the expression levels of EVE-related mRNAs were relatively limited (FPKM < 3), they displayed a broad expression pattern in the transcriptome data from different spider tissues, with a higher detection ratio in different parts of the major silk gland (Fig. S1c).

Additionally, we surveyed the small RNA (sRNA) sequencing data of the major silk glands of T. clavata, revealing highly abundant sRNAs with many baculoviral homologs. Notably, the piRNA levels of odv-ec43, p48, and p47 were significantly higher in the major silk glands, and the sRNA patterns of pif-4, lef-9, and 38k showed a typical miRNA pattern (Fig. 3b). In addition, piRNAs appeared to be mainly present in the ducts of major silk glands, whereas miRNAs were mainly present in the tail (Fig. 3b).

Identification of spider nrEVEs and related mRNAs in different spider individuals

To further confirm the prevalence of nrEVEs in individual spiders, we purified genomic DNA from the collected A. ventricosus spiders and subjected them to high-throughput sequencing. Sequencing quality was confirmed by mapping the reads to A. ventricosus reference assembly scaffolds (accession number: BGPR01), with an overall alignment rate of 82%. The coverage of the top nine scaffolds (AvEVE1-9), which contained the highest number of types of baculoviral homologs, were mapped with a coverage ranging from 42.62–98.24%, and mismatch rates ranging from 1–3.2% (Table 1). To explore how the nrEVEs were selected in the spider genome, we analyzed the ratio (Ka/Ks) of the non-synonymous substitution rate (Ka) and synonymous substitution rate (Ks) between our AvEVE sequences and reference assemblies, revealing that the baculoviral homologs in the AvEVEs underwent strong purifying selection pressure (Fig. 3c). Moreover, we performed PCR assays to identify the presence of baculoviral homologs (lef-8, p33, pif-0, and pif-1 of AvEVE2 and lef-9, pif-5, vp39, and 38k of AvEVE3) in the five A. ventricosus individuals. Remarkably, these genes were consistently detected in all the tested spiders, regardless of sex, indicating the existence of AvEVEs in A. ventricosus genome (Fig. S2).

Table 1

Alignment detail of AvEVE1-9
nrEVEs	Accession	No. of type of BHs	Length (bp)	No. of Reads^*	Mismatch^# (%)	Coverage (%)
AvEVE1	BGPR01003736.1	20	128845	3926	3.1	42.62
AvEVE2	BGPR01006520.1	20	79615	18629	1	98.24
AvEVE3	BGPR01005754.1	18	89329	226558	1.1	87.76
AvEVE4	BGPR01001437.1	7	257312	28483	3.2	56.45
AvEVE5	BGPR01014530.1	7	39014	11448	1.3	96.90
AvEVE6	BGPR01017730.1	7	32001	18304	1.5	77.52
AvEVE7	BGPR01014537.1	6	39004	68159	1.6	95.13
AvEVE8	BGPR01022151.1	6	25695	1327	2.6	72.74
AvEVE9	BGPR01031596.1	5	17954	37204	1.7	70.70
^*The read number was gained from the alignment result of sequencing data mapping to the full assemble scaffolds of BGPR01.
^#The mismatch percent represents the mismatch rate of all reads.

In addition, samples of T. clavata were used to evaluate the expression levels of nrEVEs genes in diverse tissues via quantitative reverse transcription PCR (qRT-PCR). The results demonstrated a markedly higher expression level of baculoviral homologs in the silk gland than in other tissues (Fig. 3d).

Spider nrEVEs likely originated from an unidentified ancestor NALDV

To understand the origin of spider nrEVEs, sequences of the most conserved homologs in Naldaviricetes, PIF-0, were subjected to maximum likelihood phylogenetic analysis. The resulting phylogenetic trees of PIF-0 showed that spider nrEVEs formed a monophyletic group comprising three distinct clades (Clades I, II, and III), that were evolutionarily distant from other NALDVs (Fig. 4a). Phylogenetic analyses were conducted using concatenated sequences of PIFs or LEFs from the same scaffolds. The phylogenetic tree based on concatenated sequences of PIFs and LEFs also showed that spider nrEVEs formed a monophyletic group with high confidence and that members of this group were evolutionarily distant from other NALDVs (Fig. 4b and 4c). Notably, phylogenetic trees based on concatenated sequences of PIFs or LEFs confirmed the separation of Clade I and Clade II of spider nrEVEs but lacked Clade III (Fig. 4b and 4c). This is because in contrast with Clade I and Clade II nrEVEs which contain closely distributed PIFs and LEFs, Clade III typically contains isolated nrEVE sequences.

To understand the origin of spider nrEVEs, we selected nine representative scaffolds harboring conserved PIFs and LEFs spanning lengths within the range of 80 and 500 kb for further analysis (Table S9). Among these, three belonged to Clade I, whereas the remaining six were classified as Clade II. All proteins from these nine scaffolds and representative NALDVs were analyzed via structural comparison (Table S10). Based on our structural comparison results and previously studies^11,28,29, we systematically summarized the conserved proteins of spider nrEVEs and NALDVs (Table S11). Notably, these scaffolds collectively harbored all conserved NALDV and lefaviral proteins, except AC81. In addition, they contained the most conserved baculoviral and nudiviral proteins, except for VLF-1 and P6.9, as well as some baculovirus-conserved proteins (49K, ODV-EC43, ODV-EC27, and P48) (Fig. 4d and Table S11). Five proteins conserved in spider nrEVEs (helicase-2, flap endonuclease, DNA-binding protein, nudix hydrolases, and AC106/107) were found in other NALDVs. Additionally, 11 predicted proteins were conserved in spider nrEVEs that did not exist in NALDVs (Fig. 4e), including homologs of protein kinases and innexins (Fig. 4f). Collectively, our data suggests that spider nrEVEs evolved from an unclassified ancient NALDV belonging to the order Lefavirales, of the class Naldaviricetes.

nrEVEs and spiders have engaged in long-term coevolution

To understand the relationship between spider nrEVEs and their host spiders, we systematically conducted a co-phylogenetic analysis of the PIF-0 proteins with the classification of spiders. The results showed that the evolutionary route of spider nrEVEs was highly consistent with the classification and evolutionary trajectories of spiders (Fig. 5a and Fig. S3). Specifically, the nrEVEs and transcripts from Clades I and II exclusively exist in the Entelegynae lineage, including the retrolateral tibial apophysis (RTA) clade, Orbiculariae clade, and other lineages (e.g., Palpimanoidea and Eresoidea). The nrEVEs from Clade III exist in more ancient spider lineages, such as Mesothelae, Mygalomorphae, and Haplogynae, and in a few species from the RTA and Orbiculariae clades.

Table 2

Age estimates and their 95% HPD intervals of spider nrEVEs and NALDVs
Node	Time, Mya	95% HPD, Mya	Prior	Classification
1	398.12	300.65-511.22		NALDVs
2	357.31	273.63–455.00		Lefavirus
3	353.97	265.34-458.87
4	332.86	254.63-427.89
5	329.86	248.89-428.24
6	270.17	202.96-349.58	387.00 ± 81.15³¹	Mesothelae Spider nrEVEs
7	217.36	159.03-286.51
8	207.03	162.18-259.81		Nudivirus
9	203.80	152.41-263.95
10	198.47	145.99-261.23		Baculovirus
11	179.97	136.62-230.08
12	179.42	143.04-223.59
13	179.32	126.52–242.00		Nimavirus
14	166.92	114.47-226.42		Filamentovirus
15	161.41	129.20-201.25
16	157.98	114.53-209.03
17	154.96	117.88-197.96	201.00 ± 43.15³¹	Entelegynae Spider nrEVEs Clade I&II
18	137.65	103.80-179.87
19	131.09	91.38-176.67
20	119.21	89.67-154.19		Spider nrEVEs Clade I
21	117.44	83.28-159.22
22	107.33	79.75–138.70
23	104.01	95.50-112.38	103.38 ± 4.41³⁰	Bracovirus
24	102.53	68.71-142.16		Hytrosavirus
25	101.78	73.02-136.77		Spider nrEVEs Clade II
26	95.14	63.96-129.92
27	86.33	61.37–116.10
28	81.98	53.56-118.43
29	81.48	60.50-106.94
30	81.44	56.33-111.03
31	72.99	53.05–95.63
32	63.96	42.56–90.28
33	62.14	41.18–86.96
34	53.79	36.96–72.92
35	51.30	34.18–71.95

We performed molecular dating analyses to elucidate the evolutionary events and time scales of nrEVEs in spiders, and their evolutionary history with NALDVs. By using the time to the most recent common ancestor (TMRCAs) of Chelonus inanitus bracovirus (CiBV) and Cotesia congregata bracovirus (CcBV)³⁰, and the diversion times of spiders³¹ as a calibration point, we estimated the age of spider nrEVEs from three subclades and other NALDVs using the concatenated sequences PIF-0 and LEF-8 (Fig. 5b and Table 2). Our analysis indicated that the emergence of the TMRCAs of NALDVs was ~ 398 million years ago (Mya); the TMRCAs of all lefaviruses were ~ 357 Mya; the TMRCAs of spider nrEVEs and nimaviruses were estimated to be ~ 354 Mya; the divergence times of nudivirus and baculovirus was ~ 333 Mya, which is consistent with previous studies³⁰; and the TMRCAs of hytrosavirus, filamentovirus, and AmFV were ~ 330 Mya. Moreover, our molecular dating analysis revealed that the TMRCAs of all spider nrEVEs were ~ 270 Mya, and the divergence time of the spider nrEVEs was highly consistent with those of the species to which they belonged.

Research on EVEs has predominantly focused on retroviruses, and recently identified nrEVEs also primarily originate from RNA viruses^32,33. Studies on nrEVEs from large DNA viruses have mainly concentrated on bracoviriform, which originate from nudivirus. Other NALDVs such as hytrosavirus, filamentovirus, and nimavirus, have also been found to integrate into insect and crustacean hosts^11,34. In the present study, we identified a significant number of novel nrEVEs in spiders using a combination of structural and sequence analyses (Fig. 1a). The discovery of spider nrEVEs has not only significantly broadened the known diversity of nrEVEs, but also expanded the host range of NALDVs to include two subphyla of arthropods, Chelicerata, and Mandibulata, indicating that ancient NALDVs may have widely infected all arthropods.

Based on the analysis of the distribution and evolutionary relationship of nrEVEs between spiders and NALDVs, it was speculated that an ancient lefavirus was integrated with the common ancestor of all known spiders (Fig. 6). The most related order to spiders in our study, Scorpiones, was found to have an absence of related elements (Fig. 1d), and the divergence time of Scorpiones and Araneae was ~ 397 Mya³⁵. The molecular dating time of the TMRCAs of spider nrEVEs was ~ 270 Mya (Fig. 5b), suggesting that the initial integration events might have occurred during the Devonian and Permian periods. The initially integrated viral genome may have been domesticated and coevolved with the ancestor of spiders. During long-term co-evolution, the spider nrEVEs diverge along with the spiders and formed three different clades, including the more ancient clade III and the more recent clades I and II. Horizontal transfer and integration of viral genomes have been reported to play important roles in species evolution; therefore, spider nrEVEs may play a role in the differentiation of spider species^36,37. The KaKs analysis of A. ventricosus indicated that spider nrEVEs were selected during evolution (Fig. 3c). In addition, the integration event of nudivirus in Braconidae has been confirmed to be happened at ~ 100 Mya³⁰, while the divergent time of nimavirus and endogenous nimaviruses were speculated at ~ 179 Mya in our study (Fig. 5b). Because the phylogenetic analysis indicated that nimaviruses share TMRCAs with spider nrEVEs, we assumed that the LEFs in nimaviruses might have been lost during evolution (Fig. 6). Collectively, these results confirmed that the common ancestor of NALDVs appeared ~ 400 Mya and that their integration events occurred continuously during the evolution of NALDVs and arthropods.

Whether exogenous viruses homologous to these nrEVEs exist in spiders is unclear. However, many nrEVEs in spiders maintain their genomic integrity (data not shown) during long-term co-evolution and under the pressure of purifying selection. Analysis of the genome composition of these nrEVEs revealed that, in addition to key viral genes for replication and transcription, they contain almost all the virion structure proteins conserved among various NALDVs, such as PIFs, VP39, 38K, 49K, ODV-EC27, and ODV-EC43³⁸. It is tempting to speculate that nrEVEs can potentially form complete viral particles, and spiders may use these viral particles to perform special functions, just as braconid wasps use braconid viruses to complete their parasitic lives.

However, whether the function of nrEVEs in spiders remains unknown. mRNAs related to spider nrEVEs had higher expression patterns in the main silk glands (Fig. 3d), and the sRNA patterns in different parts of the main silk glands showed obvious preferences (Fig. 3b), indicating that spider nrEVEs may function in the silk gland. Spider nrEVEs were exclusively found in Araneae rather than other orders of arachnids in the present study, suggesting that spider nrEVEs may be one of the factors shaping the highly specialized and diverse silk glands of spiders. More in-depth functional and experimental studies, such as electron microscopy and protein level verification, may provide a clearer understanding of the role of spider nrEVEs.

In previous studies, the discovery of EVEs was primarily relied on sequence searches for known viral proteins, which were limited by sequence similarity^15,18,29. Structural prediction and comparison play an important role in the discovery and classification of spider nrEVEs, suggesting that applying structural comparison could provide new insights for identifying and classifying distant viral sequences. Most recently, study based on predicted viral structural alignments have provided functional insights that were not achievable with previous approaches³⁹. Additionally, our analysis showed that spider nrEVEs exhibited significant differences in ORF coverage compared to their adjacent sequences (Fig. 2c), aligning with previous findings that endogenous sequences of viruses often exhibit heterogeneity⁴⁰. Genomic binning based on viral genome features has also recently been developed^41,42, and differences in the sequence characteristics between viral and eukaryotic genomes may serve as an effective method for identifying endogenous viral sequences.

Summarily, in this study, the widespread presence of nrEVEs in spider genomes is identified through structural comparisons and sequence analyses, highlighting the long-term co-evolution between nrEVEs and spiders. The discovery of spider nrEVEs has significantly expanded our understanding of nrEVEs in arthropods, as well as the host distribution and evolutionary history of NALDVs within these organisms. The integration of nrEVEs into the spider genome may have influenced spider evolution, more specifically in spider silk glands. These findings offer crucial insights into the specific functions, evolutionary mechanisms, and development of spider nrEVEs and silk glands. Furthermore, our study demonstrates the potential of protein structure prediction and comparison, alongside genomic features, for virus classification and the discovery of endogenous viruses.

Protein structure prediction and homologs search

Coding sequences of representative baculovirus species (Table S9) were extracted from the NCBI viral genomic database (https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/; accessed April 2023). Protein structure models were predicted using LocalColabFold⁴³ (https://github.com/YoshitakaMo/localcolabfold) with default setting, and the top-ranked structure for each sequence was selected for subsequent homologous structure searches. The AFDB50 database was obtained from https://foldseek.steinegerlab.workers.dev on June 23, 2023. Homologous structure searches in AFDB50 were conducted using Foldseek⁴⁴ (v6.29e2557) with default mode (https://github.com/steineggerlab/foldseek), and the homologous cutoff value was set at 'prob' = 0.9. The baculoviral homologous proteins identified in A. ventricosus were queried in the UniProt database, and the corresponding target scaffolds were downloaded from the NCBI WGS database.

All ORFs of the selected scaffolds in A. ventricosus, representative NALDVs and spider nrEVEs (Table S9) were extracted using orfipy⁴⁵ (v0.0.3) with parameters “--min 150 --start ATG --ignore-case.” Subsequently, the structures of these ORFs and RefSeqs of the NALDV proteins were predicted using LocalColabFold with the default setup. The predicted structures were used for “all-against-all” searches using Foldseek. The homologs were conducted through structure alignment results (‘prob’ > 0.9) using cytoscape⁴⁶, and confirmed via sequence alignments utilizing MSAProbs⁴⁷ (https://toolkit.tuebingen.mpg.de/tools/msaprobs).

Survey of baculoviral homologs in NCBI WGS and TSA database

The WGS data for Arachnida and TSA data for Araneae were obtained from NCBI in June 2023 (Table S3 and S7). Chromosomes of T. clavate were downloaded from SpiderDB 1.0 (https://spider.bioinfotoolkits.net). All ORFs in the assemble sequences were extracted using orfipy (v0.0.3) with the parameters “--min 150 and --start ATG --ignore-case.” Subsequently, the database of all peptide was built for homologous search using MMseqs2⁴⁸ (v14.7e284). The sequence of baculoviral homolog proteins identified by structure comparison of selected scaffolds in A. ventricosus were used as queries with default settings plus “--format-output query, target, evalue, bits, qseq, tseq."

Maximum likelihood phylogenetic analysis

For the phylogenetic analysis of single or concatenated PIFs and LEFs, sequences were first extracted from the MMseqs2 (v13.45111) search results using custom R script and aligned using MAFFT⁴⁹ (v7.505). Aligned sequences were trimmed using TrimAl⁵⁰ (v1.4.1) with -gt 0.5, and used for phylogenetic analyses using IQ-TREE⁵¹ (2.2.5) with ModelFinder Plus to determine the best-fit model. Support was computed from 1,000 replicates for Shimodaira–Hasegawa (SH)-like aLRT and UFBoot. As per the IQ-TREE manual, supports were deemed good when SH-like aLRT ≥ 80% and UFBoot ≥ 95%.

Molecular dating

The sequences of concatenated PIF-0 and LEF-8 were used for molecular dating analysis. The divergence times were estimated using the Bayesian inference approach implemented in BEAST (v2.7.6)⁵². The divergence times of CiBV and CcBV³⁰, and different spider families³¹ were used as calibration points, as previously described (Table. 2). We assumed uncorrelated lognormal relaxed clock and free mean substitution rates. The phylogeny was calibrated under a normal distribution and Yule process, as previously described³⁰. Monte Carlo Markov chains (MCMCs) were run for 100 million generations, sampling every 1,000 generation, using a burn-in of 1,000. Age estimates were analyzed using the FigTree (v1.4.4).

DNA and RNA extraction of spiders

The spiders (A. ventricosus and T. clavata) captured in Nanyang, Henan Province, China, were purchased via Taobao (https://shop104568564.taobao.com). For the DNA extraction of A. ventricosus, the spiders were first immersed in liquid nitrogen and macerated prior to DNA extraction. Total DNA extraction was carried out using the E.Z.N.A.® Insect DNA Kit (Omega Bio-tek, Norcross, USA) following the manufacturer’s protocol. For RNA extraction from T. clavata, four fresh tissues (legs, caphalothorax, silk glands, and abdomen without silk glands) were dissected in PBS. The tissues were macerated using a tissue shredding tube with TRIzol reagent (Invitrogen), and total RNA was extracted following the manufacturer’s protocol.

DNA sequence and data analysis

DNA from an adult A. ventricosus was used for library construction and shot-read sequencing at The Beijing Genomics Institute (BGI). The sequence data were aligned to the reference assembly of A. ventricosus (NCBI WGS accession number: BGPR01) by Bowtie2 (v2.5.3)⁵³ with default settings, and the consensus sequences of the scaffolds were extracted using samtools⁵⁴ (v1.20) and iVar⁵⁵ (v1.4.2). KaKs analyses were conducted between the consensus sequences generated by our sequencing data and the reference assemblies of AvEVE1-9 via TB-Tools (v2.085)⁵⁶. The sRNA data of major silk glands of T. clavata were downloaded from the National Genomics Data Center (accession number: CRR644805-CRR644813), after mapping to the nuclear sequence of baculoviral homologs in T. clavata by ‘Bowtie’ with a perfect match, the analysis to followe were conducted by custom R script.

PCR and qRT-PCR

For PCR assay, 2 µL spider genome DNA was used as template and PCR was conducted by Phanta Max Master Mix (Vazyme). The amplicons were electrophoresed on 0.8% agarose gel. For reverse transcription experiments, 0.5 µg total RNA was used to synthesize cDNA by HiScript III RT SuperMix (Vazyme). cDNA was used for quantitative PCR assay using the StepOne Plus Real-time PCR system (Applied Biosystems) with TB Green Premix Ex Taq II (Takara). The primers used for PCR and qPCR are listed in Table S12.

Data and code available

Databases used in our study include: NCBI RefSeq (https://ftp.ncbi.nlm.nih.gov/refseq/); NCBI WGS (Whole Genome Shotgun) genomes (https://ftp.ncbi.nlm.nih.gov/genbank/wgs/); NCBI TSA (Transcriptome Shotgun Assemblies) (https://ftp.ncbi.nlm.nih.gov/genbank/tsa/); UniRef90 (https://ftp.ebi.ac.uk/pub/databases/uniprot/uniref/uniref90/); National Genomics Data Center (https://ngdc.cncb.ac.cn/); AlphaFold Protein Structure Database (https://alphafold.ebi.ac.uk/); and SpiderDB 1.0 (https://spider.bioinfotoolkits.net). The high-throughput sequencing raw data of A. ventricosus genomic DNA were deposited into the Genome Sequence Archive (GSA) of the National Genomics Data Center and are available through accession number: CRA017441. The original data generated in our study, including phylogenetic trees, predicted protein structures, custom codes, etc., have been made publicly available at Figshare.

Author contributions

H.H., Z.H., and M.W. conceptualized and designed the study. H.H. performed all computational analyses and experiments. H.H., Z.H., M.W., and J.M.V. interpreted the results and wrote the manuscript. All authors reviewed and edited the manuscript.

Acknowledgement

This study was supported by the grants from the National Natural Science Foundation of China (No. 37277104 to H.H., 32370182 to M.W., and 32170156 to Z.H.). We thank all the databases and projects that provide the public sequencing data for this study.

Feschotte, C. & Gilbert, C. Endogenous viruses: insights into viral evolution and impact on host biology. Nature Reviews Genetics 13, 283-296 (2012). https://doi.org/10.1038/nrg3199
Holmes, E. C. The evolution of endogenous viral elements. Cell Host Microbe 10, 368-377 (2011). https://doi.org/10.1016/j.chom.2011.09.002
Weiss, R. A. The discovery of endogenous retroviruses. Retrovirology 3, 67 (2006). https://doi.org/10.1186/1742-4690-3-67
Horie, M. et al. Endogenous non-retroviral RNA virus elements in mammalian genomes. Nature 463, 84-87 (2010). https://doi.org/10.1038/nature08695
Houe, V., Bonizzoni, M. & Failloux, A. B. Endogenous non-retroviral elements in genomes of Aedes mosquitoes and vector competence. Emerg Microbes Infect 8, 542-555 (2019). https://doi.org/10.1080/22221751.2019.1599302
Huang, H.-J. et al. Co-option of a non-retroviral endogenous viral element in planthoppers. Nature Communications 14 (2023). https://doi.org/10.1038/s41467-023-43186-2
Frank, J. A. & Feschotte, C. Co-option of endogenous viral sequences for host cell function. Curr Opin Virol 25, 81-89 (2017). https://doi.org/10.1016/j.coviro.2017.07.021
Aiewsakun, P. & Katzourakis, A. Endogenous viruses: Connecting recent and ancient viral evolution. Virology 479-480, 26-37 (2015). https://doi.org/10.1016/j.virol.2015.02.011
Johnson, W. E. Origins and evolutionary consequences of ancient endogenous retroviruses. Nature Reviews Microbiology 17, 355-370 (2019). https://doi.org/10.1038/s41579-019-0189-2
Regier, J. C. et al. Arthropod relationships revealed by phylogenomic analysis of nuclear protein-coding sequences. Nature 463, 1079-1083 (2010). https://doi.org/10.1038/nature08742
van Oers, M. M. et al. Developments in the classification and nomenclature of arthropod-infecting large DNA viruses that contain pif genes. Arch Virol 168, 182 (2023). https://doi.org/10.1007/s00705-023-05793-8
Koonin, E. V. et al. Global Organization and Proposed Megataxonomy of the Virus World. Microbiol Mol Biol R 84 (2020). https://doi.org/ARTN e00061-19 10.1128/MMBR.00061-19
Passarelli, A. L. & Guarino, L. A. Baculovirus late and very late gene regulation. Curr Drug Targets 8, 1103-1115 (2007). https://doi.org/Doi 10.2174/138945007782151324
Lin, C. L. et al. Persistent Hz-1 virus infection in insect cells: evidence for insertion of viral DNA into host chromosomes and viral infection in a latent status. J Virol 73, 128-139 (1999). https://doi.org/10.1128/JVI.73.1.128-139.1999
Cheng, R. L., Li, X. F. & Zhang, C. X. Nudivirus Remnants in the Genomes of Arthropods. Genome Biol Evol 12, 578-588 (2020). https://doi.org/10.1093/gbe/evaa074
Bezier, A. et al. Polydnaviruses of braconid wasps derive from an ancestral nudivirus. Science 323, 926-930 (2009). https://doi.org/10.1126/science.1166788
Strand, M. R. & Burke, G. R. Polydnaviruses: Evolution and Function. Current Issues in Molecular Biology 34, 163-181 (2019). https://doi.org/10.21775/cimb.034.163
Kawato, S. et al. Crustacean Genome Exploration Reveals the Evolutionary Origin of White Spot Syndrome Virus. J Virol 93 (2019). https://doi.org/10.1128/JVI.01144-18
Guinet, B. et al. Endoparasitoid lifestyle promotes endogenization and domestication of dsDNA viruses. Elife 12 (2023). https://doi.org/10.7554/eLife.85993
Legg, D. A., Sutton, M. D. & Edgecombe, G. D. Arthropod fossil data increase congruence of morphological and molecular phylogenies. Nat Commun 4, 2485 (2013). https://doi.org/10.1038/ncomms3485
Kulkarni, S., Wood, H. M. & Hormiga, G. Advances in the reconstruction of the spider tree of life: A roadmap for spider systematics and comparative studies. Cladistics 39, 479-532 (2023). https://doi.org/10.1111/cla.12557
Sanggaard, K. W. et al. Spider genomes provide insight into composition and evolution of venom and silk. Nat Commun 5, 3765 (2014). https://doi.org/10.1038/ncomms4765
Mihara, T. et al. Linking Virus Genomes with Host Taxonomy. Viruses 8, 66 (2016). https://doi.org/10.3390/v8030066
Jumper, J. et al. Highly accurate protein structure prediction with AlphaFold. Nature 596, 583-589 (2021). https://doi.org/10.1038/s41586-021-03819-2
van Kempen, M. et al. Fast and accurate protein structure search with Foldseek. Nat Biotechnol (2023). https://doi.org/10.1038/s41587-023-01773-0
Lin, Z. et al. Evolutionary-scale prediction of atomic-level protein structure with a language model. Science 379, 1123-1130 (2023). https://doi.org/10.1126/science.ade2574
Hu, W. et al. A molecular atlas reveals the tri-sectional spinning mechanism of spider dragline silk. Nat Commun 14, 837 (2023). https://doi.org/10.1038/s41467-023-36545-6
Guinet, B. et al. A novel and diverse family of filamentous DNA viruses associated with parasitic wasps. Virus Evol 10, veae022 (2024). https://doi.org/10.1093/ve/veae022
Ter Horst, A. M., Nigg, J. C., Dekker, F. M. & Falk, B. W. Endogenous Viral Elements Are Widespread in Arthropod Genomes and Commonly Give Rise to PIWI-Interacting RNAs. J Virol 93 (2019). https://doi.org/10.1128/JVI.02124-18
Thézé, J., Bézier, A., Periquet, G., Drezen, J.-M. & Herniou, E. A. Paleozoic origin of insect large dsDNA viruses. Proceedings of the National Academy of Sciences 108, 15931-15935 (2011). https://doi.org/10.1073/pnas.1105580108
Bond, Jason E. et al. Phylogenomics Resolves a Spider Backbone Phylogeny and Rejects a Prevailing Paradigm for Orb Web Evolution. Current Biology 24, 1765-1771 (2014). https://doi.org/10.1016/j.cub.2014.06.034
Wallau, G. L. RNA virus EVEs in insect genomes. Curr Opin Insect Sci 49, 42-47 (2022). https://doi.org/10.1016/j.cois.2021.11.005
Frank, J. A. et al. Evolution and antiviral activity of a human protein of retroviral origin. Science 378, 422-428 (2022). https://doi.org/10.1126/science.abq7871
Gilbert, C. & Belliardo, C. The diversity of endogenous viral elements in insects. Curr Opin Insect Sci 49, 48-55 (2022). https://doi.org/10.1016/j.cois.2021.11.007
Kumar, S. et al. TimeTree 5: An Expanded Resource for Species Divergence Times. Mol Biol Evol 39 (2022). https://doi.org/ARTN msac174 10.1093/molbev/msac174
Arnold, B. J., Huang, I. T. & Hanage, W. P. Horizontal gene transfer and adaptive evolution in bacteria. Nat Rev Microbiol 20, 206-218 (2022). https://doi.org/10.1038/s41579-021-00650-4
Gilbert, C. & Cordaux, R. Viruses as vectors of horizontal transfer of genetic material in eukaryotes. Curr Opin Virol 25, 16-22 (2017). https://doi.org/10.1016/j.coviro.2017.06.005
Jia, X. et al. Architecture of the baculovirus nucleocapsid revealed by cryo-EM. Nat Commun 14, 7481 (2023). https://doi.org/10.1038/s41467-023-43284-1
Nomburg, J. et al. Birth of protein folds and functions in the virome. Nature (2024). https://doi.org/10.1038/s41586-024-07809-y
Hackl, T., Duponchel, S., Barenhoff, K., Weinmann, A. & Fischer, M. G. Virophages and retrotransposons colonize the genomes of a heterotrophic flagellate. Elife 10 (2021). https://doi.org/10.7554/eLife.72674
Kieft, K., Adams, A., Salamzade, R., Kalan, L. & Anantharaman, K. vRhyme enables binning of viral genomes from metagenomes. Nucleic Acids Res 50, e83 (2022). https://doi.org/10.1093/nar/gkac341
Nayfach, S. et al. CheckV assesses the quality and completeness of metagenome-assembled viral genomes. Nat Biotechnol 39, 578-585 (2021). https://doi.org/10.1038/s41587-020-00774-7
Mirdita, M. et al. ColabFold: making protein folding accessible to all. Nat Methods 19, 679-682 (2022). https://doi.org/10.1038/s41592-022-01488-1
van Kempen, M. et al. Fast and accurate protein structure search with Foldseek. Nat Biotechnol 42, 243-246 (2024). https://doi.org/10.1038/s41587-023-01773-0
Singh, U. & Wurtele, E. S. orfipy: a fast and flexible tool for extracting ORFs. Bioinformatics 37, 3019-3020 (2021). https://doi.org/10.1093/bioinformatics/btab090
Shannon, P. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 13, 2498-2504 (2003). https://doi.org/10.1101/gr.1239303
Gabler, F. et al. Protein Sequence Analysis Using the MPI Bioinformatics Toolkit. Curr Protoc Bioinformatics 72, e108 (2020). https://doi.org/10.1002/cpbi.108
Steinegger, M. & Soding, J. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat Biotechnol 35, 1026-1028 (2017). https://doi.org/10.1038/nbt.3988
Katoh, K., Misawa, K., Kuma, K. & Miyata, T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res 30, 3059-3066 (2002). https://doi.org/10.1093/nar/gkf436
Capella-Gutierrez, S., Silla-Martinez, J. M. & Gabaldon, T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25, 1972-1973 (2009). https://doi.org/10.1093/bioinformatics/btp348
Minh, B. Q. et al. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era. Mol Biol Evol 37, 1530-1534 (2020). https://doi.org/10.1093/molbev/msaa015
Suchard, M. A. et al. Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol 4, vey016 (2018). https://doi.org/10.1093/ve/vey016
Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat Methods 9, 357-359 (2012). https://doi.org/10.1038/nmeth.1923
Danecek, P. et al. Twelve years of SAMtools and BCFtools. Gigascience 10 (2021). https://doi.org/10.1093/gigascience/giab008
Grubaugh, N. D. et al. An amplicon-based sequencing framework for accurately measuring intrahost virus diversity using PrimalSeq and iVar. Genome Biol 20, 8 (2019). https://doi.org/10.1186/s13059-018-1618-7
Chen, C. et al. TBtools-II: A "one for all, all for one" bioinformatics platform for biological big-data mining. Mol Plant 16, 1733-1742 (2023). https://doi.org/10.1016/j.molp.2023.09.010

There is NO Competing Interest.

Supplementarytables.xlsx
Supplementary tables
rs2.pdf
Reporting Summary
Supplementaryinformation.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Discovery of novel non-retroviral endogenous viral elements reveals their long-term co-evolution with spiders

Status:

Version 1

Abstract

Figures

Introduction

Results

Discussion

Methods

DNA and RNA extraction of spiders

PCR and qRT-PCR

Declarations

Data and code available

Author contributions

Acknowledgement

References

Additional Declarations

Supplementary Files

Status:

Version 1