Introduction
============

Retroviruses (family *Retroviridae*) are viruses of major medical significance as some are associated with severe infectious disease or are oncogenic ([@vez001-B10]; [@vez001-B2]; [@vez001-B25]). Retroviruses are also of note because of their ability to integrate into the host germ-line, generating endogenous retroviruses (ERVs) that then exhibit Mendelian inheritance ([@vez001-B21]; [@vez001-B13]). ERVs are widely distributed in vertebrates ([@vez001-B11]; [@vez001-B4]; [@vez001-B10]; [@vez001-B25]) and constitute important molecular fossils for the study of retrovirus evolution. ERVs related to all seven major retroviral genera (alpha-, beta-, delta-, epsilon-, gamma-, lenti-, and spuma-) have been described ([@vez001-B10]), although some of the more complex retroviruses, such as lenti-, delta- and foamy viruses, rarely appear as endogenous copies.

As well as being agents of disease, foamy viruses are of importance because they exhibit long-term virus-host co-divergence ([@vez001-B22]). Endogenous foamy viruses (EFVs), first discovered in sloths (class Mammalia) ([@vez001-B16]) also co-diverge with their hosts, and have also been reported in primates and the Cape golden mole ([@vez001-B8], [@vez001-B9]). The subsequent discovery of a EFV in the coelacanth genome indicated that foamy viruses may have an ancient evolutionary history ([@vez001-B7]), likely spanning hundreds of million years ([@vez001-B2]). Although EFVs or foamy-like elements have been reported in fish, amphibians, and mammals, to date they have not been reported in genomes of two other major classes of vertebrates---reptiles and birds ([@vez001-B23]; [@vez001-B12]; [@vez001-B10]; [@vez001-B25]).

Materials and methods
=====================

Genomic mining and consensus genome construction
------------------------------------------------

To identify foamy viruses in reptiles, the TBLASTN program ([@vez001-B3]) was used to screen relevant taxa from 28 reptile genomes ([Supplementary Table S1](#sup1){ref-type="supplementary-material"}) and 130 bird genomes ([Supplementary Table S2](#sup1){ref-type="supplementary-material"}) (as of October 2018) downloaded from GenBank ([www.ncbi.nlm.nih.gov/genbank](http://www.ncbi.nlm.nih.gov/genbank)). In each case, the amino acid sequences of the Pol genes of representative EFVs, endogenous foamy-like viruses, and exogenous foamy viruses were chosen as queries ([Supplementary Table S3](#sup1){ref-type="supplementary-material"}). As filters to identify significant and meaningful hits, we chose sequences with more than 30 per cent amino acid identity over a 30 per cent region, with an *e*-value set to 0.00001. Genomes that contained only single hits for EFVs were excluded as likely false positives. We extended viral flanking sequences of the hits to identify the 5′- and 3′-long terminal repeats (LTRs) using LTR finder ([@vez001-B26]) and LTR harvest ([@vez001-B5]). Sequences highly similar to foamy virus proteins found in tuatara were termed 'ERV-Spuma.*n*-Spu' (in which *n* represents the number of the sequence extracted from this tuatara genome) according a recently proposed nomenclature for ERVs (Gifford et al. 2018), and aligned to generate an ERV-Spuma-Spu consensus genome (ERV-Spuma.0-Spu) ([Supplementary Table S4](#sup1){ref-type="supplementary-material"}). Conserved domains were identified using CD-Search service in NCBI ([@vez001-B19]).

Molecular dating of integration times
-------------------------------------

The ERV integration time can be estimated using the following simple relation: *T* = (*D*/*R*)/2, in which *T* is the integration time (million years, MY), *D* is the number of nucleotide differences per site between the two LTRs, and *R* is the genomic substitution rate (i.e. number of nucleotide substitutions per site, per year). We used the previously estimated neutral substitution rate for squamate reptiles (7.6 × 10^−10^ nucleotide substitutions per site, per year) ([@vez001-B20]). LTRs less than 300 bp in length were not included in this analysis. Five pairwise LTRs were used for date estimation ([Supplementary Table S5](#sup1){ref-type="supplementary-material"}).

Phylogenetic analysis
---------------------

To determine the evolutionary relationship of EFVs and retroviruses, sequences of the Pol proteins were aligned using MAFFT 7.222 ([@vez001-B15]) and confirmed manually in MEGA7 ([@vez001-B18]). The phylogenetic relationships among these sequences were then determined using the maximum-likelihood (ML) method in PhyML 3.1 ([@vez001-B6]), incorporating 100 bootstrap replicates to determine node robustness. The best-fit models of amino acid substitution were determined by ProtTest 3.4.2 ([@vez001-B1]): RtREV + Γ+I for Pol, and LG + Γ+I + F for concatenated Gag, Pol, and Env. All alignments used in the phylogenetic analyses can be found in [Supplementary data](#sup1){ref-type="supplementary-material"} sets S1 and S2.

Results and discussion
======================

Discovery of foamy viral elements in reptile genomes
----------------------------------------------------

We screened all available reptilian and bird genomes by using the TBLASTN algorithm with various foamy viruses, including EFVs, as screening probes. We only considered viral hits within long genomic scaffold (\>20 kilobases in length) to be *bona fide* ERVs. This genomic mining identified 118 ERV hits in tuatara (*Sphenodon punctatus*) and none in bird genomes. Hence, a total of 118 ERV hits in the tuatara genome were extracted and subjected to evolutionary analysis ([Supplementary Table S6](#sup1){ref-type="supplementary-material"}) and these ERVs were named as ERV-Spuma.*n*-Spu (where *n* = 1∼118).

Genomic organization
--------------------

We extracted all significant foamy viral elements and constructed a consensus genomic sequence of ERV-Spuma-Spu ([Supplementary Fig. S1](#sup1){ref-type="supplementary-material"}, [Table S4](#sup1){ref-type="supplementary-material"}), termed ERV-Spuma.0-Spu. The consensus genome harboured a pairwise LTRs and exhibits a typical spuma virus structure, encoding three main open reading frames (ORF)---*gag*, *pol*, and *env*---and one putative additional accessory gene, ORF 1 ([Fig. 1](#vez001-F1){ref-type="fig"}). Interestingly, this accessory ORF 1 exhibits no sequence similarity to known foamy accessory genes. Notably, by searching the Conserved Domains Database ([www.ncbi.nlm.nih.gov/Structure/cdd](http://www.ncbi.nlm.nih.gov/Structure/cdd)), we identified three typical foamy conserved domains for both the consensus and one of two full-length original ERV-Spuma.23-Spu: (1) Spuma virus Gag domain (pfam03276) ([@vez001-B24]), (2) Spuma aspartic protease (A9) domain (pfam03539) which exists in all mammalian foamy virus Pol protein ([@vez001-B2]), and (3) foamy virus envelope protein domain (pfam03408) ([@vez001-B7]) ([Supplementary Figs S2 and S3](#sup1){ref-type="supplementary-material"}), confirming that ERV-Spuma-Spu is indeed of foamy virus origin.

![Genomic organizations of ERV-Spuma-Spu. LTR, long-terminal repeat; PBS, primer-binding site; Pro, aspartic protease; RT, reverse transcriptase; RH, ribonuclease H; IN, integrase.](vez001f1){#vez001-F1}

Estimated integration times
---------------------------

To broadly estimate the integration time of ERV-Spuma-Spu, we utilized the LTR-divergence method which analyzes the degree of divergence between 5′ and 3′LTRs assuming a known rate of nucleotide substitution ([@vez001-B14]). In total, five pairwise LTRs flanking ERV-Spuma-Spu elements were used for date estimation ([Supplementary Table S5](#sup1){ref-type="supplementary-material"}), from which we estimated an integration time of ERV-Spuma-Spu ranging from 1.3 to 35.47 MYA (million years ago). Although these dates are young relative to the age of reptiles, LTR dating may severely underestimate ERV ages ([@vez001-B17]; [@vez001-B2]), such that all estimates of integration time should be treated with caution.

Evolutionary relationships of ERVs-Spuma-Spu
--------------------------------------------

Sequences of the Pol protein (490 amino acids in length) of ERV-Spuma-Spu were used for phylogenetic analysis. Our ML phylogenetic trees revealed that the EFVs present in the tuatara genome formed a close monophyletic group within the foamy virus clade, indicative of a single origin, and with high bootstrap support ([Fig. 2](#vez001-F2){ref-type="fig"}). The divergent phylogenetic position of ERV-Spuma-Spu is compatible with virus-host co-divergence for the entire history of the vertebrates. However, it is possible that this pattern will change with a larger sampling of taxa such that the EFV phylogeny expands.

![Phylogenetic tree of retroviruses, including ERV-Spuma-Spu, inferred using amino acid sequences of the Pol gene (490aa). The tree is midpoint rooted for clarity only. The newly identified ERV-Spuma-Spu sequences are labelled using a grey-shaded box with their accession numbers (different pol sequences in same contig are numbered in the suffix). The scale bar indicates the number of amino acid changes per site. Bootstrap values \<70 per cent are not shown. The alignment of pol amino acid sequences is provided in [Supplementary data set S1.](#sup1){ref-type="supplementary-material"}](vez001f2){#vez001-F2}

In addition, we inferred a phylogenetic tree of FVs, EFVs, and foamy-like ERVs. This was consistent to those in previous studies ([@vez001-B2]) and revealed that the ERVs reported previously in the tuatara genome are distantly related to foamy virus and belong to class III retroviruses ([Supplementary Fig. S4](#sup1){ref-type="supplementary-material"}) ([@vez001-B23]; [@vez001-B12]). Failure to detect any ERV-Spuma-Spu-related elements in the remaining reptilian genomes suggest that the virus may be not vertically transmitted among reptiles, although this will clearly need to be reassessed with a larger sample size. The absence of EFVs in avian genomes is puzzling and clearly merits additional study.

Previous studies provided strong evidence for the co-divergence of foamy viruses and their vertebrate hosts over extended periods of evolutionary time ([@vez001-B16]). An analysis of concatenated gag-pol-env protein sequences suggests that the reptilian ERV-Spuma-Spu newly described here follow the same pattern of long-term virus-host co-divergence ([Fig. 3](#vez001-F3){ref-type="fig"}). As such, these data imply that ERV-Spuma-Spu may have diverged from the other mammalian foamy viruses along with its tuatara host more than 320 MYA (<http://www.timetree.org/>). The discovery of ERVs-Spuma-Spu therefore fills a major gap in our understanding of the taxonomic distribution of the foamy viruses and their evolutionary history.

![A simplified evolutionary relationship between foamy viruses (left) and their vertebrate hosts (right). The scale bar indicates number of amino acid changes per site in the viruses or host speciation time (million years ago, MYA). Bootstrap support values are provided at each node. The alignment of the FV Gag-Pol-Env amino acid sequences is provided in [Supplementary data set S2.](#sup1){ref-type="supplementary-material"}](vez001f3){#vez001-F3}

Supplementary Material
======================

###### 

Click here for additional data file.

J.C. is supported by National Natural Science Foundation of China (31671324) and CAS Pioneer Hundred Talents Program. E.C.H. is supported by an ARC Australian Laureate Fellowship (FL170100022).

Data availability
=================

All the data needed to generate the conclusions made in the article are present in the article itself and/or the [Supplementary data](#sup1){ref-type="supplementary-material"}. Additional data related to this article may be requested from the authors.

**Conflict of interest:** None declared.

[^1]: Xiaoman Wei and Yicong Chen authors contributed equally to this work.
