An Evolutionarily Young Polar Bear (Ursus maritimus) Endogenous Retrovirus Identified from Next Generation Sequence Data

Transcriptome analysis of polar bear (Ursus maritimus) tissues identified sequences with similarity to Porcine Endogenous Retroviruses (PERV). Based on these sequences, four proviral copies and 15 solo long terminal repeats (LTRs) of a newly described endogenous retrovirus were characterized from the polar bear draft genome sequence. Closely related sequences were identified by PCR analysis of brown bear (Ursus arctos) and black bear (Ursus americanus) but were absent in non-Ursinae bear species. The virus was therefore designated UrsusERV. Two distinct groups of LTRs were observed including a recombinant ERV that contained one LTR belonging to each group indicating that genomic invasions by at least two UrsusERV variants have recently occurred. Age estimates based on proviral LTR divergence and conservation of integration sites among ursids suggest the viral group is only a few million years old. The youngest provirus was polar bear specific, had intact open reading frames (ORFs) and could potentially encode functional proteins. Phylogenetic analyses of UrsusERV consensus protein sequences suggest that it is part of a pig, gibbon and koala retrovirus clade. The young age estimates and lineage specificity of the virus suggests UrsusERV is a recent cross species transmission from an unknown reservoir and places the viral group among the youngest of ERVs identified in mammals.


Introduction
Retroviruses are a large and diverse family of single stranded positive sense RNA viruses [1,2]. Reverse transcription of retroviral RNA genomes into double stranded DNA (dsDNA) is a key step in the retroviral replication cycle. The resulting dsDNA genome is then integrated into the genome of the infected organism in the form of a provirus [2][3][4]. Proviral integration occasionally occurs in the germline such that viral transmission changes from horizontal (via infection) to vertical as a Mendelian trait [4]. Vertically transmitted retroviruses are designated endogenous retroviruses (ERVs) and represent the successful colonization of the host germline [5]. The result of a multitude of such colonizations is that a large percentage of vertebrate genomes are comprised of ERVs [6].
A blastn search was performed on all 72,214 scaffolds of the draft polar bear genome sequence using, as query sequences, brain derived transcriptome retroviral sequences with the highest percentage similarities to MDEV and PERV [21]. MDEV and PERV sequences attracted the most non-overlapping transcriptome reads, deriving from different parts of the retroviral genome, from the RNA-seq virome data and were therefore further analyzed [23]. The blastn search revealed four loci in four different polar bear scaffolds with significant similarity to the potential retroviral sequences (Table S1). Presumed endogenous retrovirus harboring loci were extracted from the respective polar bear scaffold sequences including~15,000 bp both up and downstream of the blastn identified scaffold regions. The exact boundaries of the identified ERV loci were then determined using Repeatmasker (Table S2), and additional sequence comparisons, as described in [32][33][34]. Repeatmasker results indicated that the identified proviruses belonged to the gammaretrovirus genus of retroviruses and were annotated as koala retrovirus (KoRV) sequences indicating the high similarity to KoRV, specifically KoRV_I-int (Table S2) [35,36].
The Retrotector script was used as an independent verification method for the presence of retroviral motifs in the polar bear scaffolds. Retrotector analysis revealed retroviral motifs similar to known gammaretroviral sequences [33,34,37] identifying proviral genes gag, protease (pro), polymerase (pol), and envelope (env) as shown in Figure 1. A dUTPase domain was not identified. Retrotector generated reconstructed retroviral proteins (called putein) sequences from the defective reading frames of the proviral loci by introducing 10 frameshifts. GAG and PRO proteins were predicted from Retrotector for all four scaffolds with a percentage identity of 94.8% and 80%, respectively. The majority of the differences found in the alignments were identified in scaffold 7 that also contain three frameshifts in the GAG putein ( Figure 1). Excluding scaffold 7 sequences, the percentage sequence identity increased to 99% for both alignments. Putein sequences for the pol genes of scaffold 1 and 162 were generated from provirus regions harboring assembly gaps and were therefore incomplete, while proviral sequences within scaffold 7 and scaffold 200 appeared to be intact with regard to the pol gene region with a percentage identity of 95.6% (Figure 1). Retrotector analysis corrected 5 frameshifts during the scaffold 1 POL putein prediction, further indicating the incompleteness of the identified ORF. The pol gene in scaffold 1 was interrupted by an inserted Gypsy like element and a SINE-MIRb element encompassing nucleotides 4652-5549 of the provirus.
Motifs conserved in retroviral envelope proteins were detected by Retrotector in all proviral sequences, with the exception of the scaffold 1 proviral sequence that was lacking the env gene region. The env gene putein domains in scaffolds 7 and 162 appeared to be affected by mutations. Generation of puteins for both scaffolds identified frameshifts, whereas the provirus in scaffold 200 displayed an intact ENV ORF ( Figure 1). Envelope puteins, as generated by Retrotector, displayed 91.8% among sequence identity for alignable regions.
The National Center for Biotechnology Information (NCBI) Conserved Domain database analysis of proviral scaffolds identified the same motifs and frameshifts as Retrotector for all four proviral sequences and further highlighted the degree of similarity with respect to conserved motifs in other retroviral proteins [38]. Further analyses of all proviral loci indicated tRNA Pro as the primer binding site for viral replication initiation. Sequence analysis further indicated that proviral loci in scaffold 1, 7 and 162 are unlikely to produce functional viruses due to frameshifts, deletions, insertions, and repetitive elements that interrupt the sequences. The provirus in scaffold 200, on the other hand, harbored intact retroviral motifs in the correct reading frame, an indication that functional viral proteins could potentially be produced from this locus ( Figure 1).

Figure 1.
Multiple sequence alignments of the four UrsusERV proviral loci as identified in the polar bear genome sequence are shown. Proviral long terminal repeat (LTR), gag, pro, pol and env regions are indicated. Scaffold 1 and scaffold 162 proviral sequences harbor assembly gaps and lack proviral pol and env regions, indicated by horizontal black lines. The scaffold 200 provirus is structurally complete and could potentially produce functional proteins. Sequence differences relative to the consensus sequence are indicated by black vertical lines with the scaffold 7 provirus exhibiting the greatest overall divergence. Only the scaffold 200 provirus harbors a complete env gene that is partially present in scaffolds 162 and scaffold 7 and completely absent from scaffold 1.
A majority rule proviral consensus sequence of 7878 nt was generated from the four identified proviral loci. A consensus sequence that included LTR, gag, pro, pol, and partial env gene characteristics was used for additional blastn searches of the polar bear draft genome sequence scaffolds. Blastn search using the consensus sequence identified no additional proviral loci. However, 15 solitary LTRs were identified. A multiple sequence alignment of identified UrsusERV LTR sequences (proviral and solitary) revealed two subgroups with group-specific nucleotide differences that we designated LTR-A and LTR-B ( Figure 2). Out of a total of 23 LTRs, 12 belonged to the LTR-A subgroup, 10 belonged to the LTR-B subgroup and one appeared to be a recombinant. The scaffold 1 and 200 proviruses were flanked by LTR-A and the scaffold 7 provirus was flanked by LTR-B.
We observed that the scaffold 162 provirus was flanked by a 5 1 LTR characteristic of an LTR-A sequence, yet the 3 1 LTR was characteristic of an LTR-B sequence for the 5 1 most 200 nt and an LTR-A-like sequence for the 3 1 most 300 nt. UrsusERV LTR sequences were therefore analyzed for recombination events using four computational methods, specifically, GARD, DualBrothers, Recombination Analysis Tool (RAT, Norwich, UK) and RECCO (Saarbrücken, Germany) [39][40][41][42][43]. All four methods identified a recombination event for the 3 1 LTR in the scaffold 162 provirus occurring between nt 177 and 207 of the 3 1 LTR sequence ( Figure S1). Phylogenetic analysis of the GARD recombination script clustered the first sub-region of scaffold 162 LTR-A/B ranging from 1 to 177 bp within the LTR-B clade, while the second region ranging from 207 bp onwards within the LTR-A clade, supporting a recombination event ( Figure S1). The RECCO computational method confirmed the presumably recombined LTR sequence ("scaffold162_1-18,384_LTR-A/B" in Figure S1) with a p-value < 0.001. RECCO also indicated another potential recombination event in the Scaffold 143 LTR (Table 1).
Target site duplications flanking proviral LTRs were also examined. Retroviral integration generally generates a target site duplication (TSD) of 4-5 bp immediately flanking the provirus [44]. Four bp TSDs were identified for the four proviral insertions and in 10 out of 15 solo LTRs ( Table 2). The remaining five LTR sequences were identified in small scaffold sequences with TSD information being unavailable in the current draft genome version.

UrsusERV Age Estimation
Age estimation was based on the sequence divergence of the proviral 5 1 and 3 1 LTR sequences. Once the provirus is integrated into the genome of the host, the 5 1 and 3 1 LTRs are identical in sequence [37,38]. Accumulation of differences between the proviral 5 1 and 3 1 LTR will subsequently occur at the mutational rate of the host genome. Thus, sequence divergence between the proviral 5 1 and 3 1 LTRs can be used to estimate the age of a given provirus [45]. Using an evolutionary rate of 0.0015/nt/myr, as previously determined for bears [25], provided integration age estimations of up to two million years old with the provirus in scaffold 200 lacking differences between its two LTR sequences making an age estimate based on a molecular clock difficult.
The oldest provirus located in scaffold 7 was estimated to be roughly 2 myr old ( Table 3). The scaffold 162 provirus could not be dated with confidence, as the 3 1 LTR is a recombinant between LTR-A and LTR-B types. Proviral distribution analysis of UrsusERV in ancient polar bear whole genome sequence data (130,000 years old fossil) identified sequence reads spanning the Scaffold 200 LTR ends and flanking regions. Presence of the youngest provirus in the ancient polar bear genome suggests it integrated in the polar bear lineage at least 130,000 thousand years ago and likely earlier [31]. We stress in this context that age estimations based on LTR divergence should be regarded with caution, especially when regarding evolutionarily younger proviruses, as this dating method can be biased by, for instance, gene conversion thus potentially underestimating proviral ages [46,47].

UrsusERV Molecular and Genome Screening of Bear Species
We further investigated the presence of UrsusERV sequences in additional polar bears and other bear species. Majority rule consensus sequences of proviral and LTR sequences were used to screen the giant panda draft genome assembly (BGI-Shenzhen AilMel 1.0 December 2009) available through the University of California, Santa Cruz (UCSC) Genome Browser [48]. Blat searches failed to identify significantly similar sequences indicating that UrsusERV is not present in all members of the Ursidae family, further supporting the age estimations based on LTR sequence divergence. We next examined other bear species, specifically whole genome sequencing data from black bear (Ursus americanus), brown bear (Ursus arctos), grizzly bear (Ursus arctos ssp.), an ancient polar bear jawbone [28,31], and six unrelated modern polar bears that were obtained from the National Centre for Biotechnology Information (NCBI) Sequence Read Archive [27]. Short read data were transformed to fastq files and mapped to the four polar bear retroviral scaffolds using the Burrows-Wheeler aligner [49]. Analyses of the polar bear mapped reads indicated the presence of all four proviruses in the six polar bear whole genome sequenced animals. Reads identical to UrsusERV proviral sequence were also identified in the ancient polar bear jawbone sequence data. Genome sequence data from the black and brown bears identified sequences with 93.9% identity to the identified consensus proviral sequence in both bear species indicating that UrsusERV is present in the Ursinae clade of the Ursidae family [50]. This also suggests integration of UrsusERV like viruses in bears occurred earlier than 2 million years ago (mya, the oldest dated provirus) as the black bear and brown/polar bear clades diverged from a common ancestor over 7 mya.
All the UrsusERV sequences flanking the polar bear proviral LTRs in the polar bear draft genome sequence were identified in the other six polar bear whole genome draft sequences. The scaffold 200 3 1 proviral insertion site was the only one identified in the 130,000-year old polar bear jawbone sample providing evidence that the same insertion site for at least one provirus was present in the ancient polar bear population. Sequence coverage was not as comprehensive as for the modern bear sequences and thus we cannot exclude that the other proviruses were present in this sample [31]. Brown and black bears whole genome draft sequences were also screened for the proviral insertion sites identified in polar bears. While the proviral flanking sequences could be identified in the black bear genome, they lacked a proviral integration e.g., the flanking sequences were uninterrupted by UrsusERV. Therefore, while black bears are positive for UrsusERV sequences, the virus has integrated in distinct locations in the genome relative to polar bears. This could in part explain the differences in age estimates between the polar bear derived ERVs and the divergence time of black bear and brown bear/polar bear lineages in addition to the overall lower percentage identity of UrsusERV derived from black bears relative to those found in brown or polar bears. Out of the four female brown bears examined, two had the same UrsusERV insertion sites as polar bears for all the proviruses except for the provirus located in scaffold 200, which was polar bear specific as all brown bears tested lacked an integration at this genomic location although the flanking sequences were identified (Table S3). Two additional brown bear genome sequences, from Kenai and Admiralty Island, where found to have the same provirus insertions found on polar bear scaffolds 7 and 162, while the polar bear scaffold 1 integration was absent indicating insertion site polymorphisms in brown bears for the scaffold 1 proviral integration (Table S3) [51].
Molecular screening for UrsusERV sequences was also performed by sequence-specific PCR using genomic DNA extracted from a giant panda, spectacled bear (Tremarctos ornatus), Syrian brown bear (Ursus arctos syriacus), black bear, and polar bear to further verify the bioinformatic analysis. Genomic DNA was subjected to PCRs using primers designed to amplify a 312 bp region within the gag gene region showing little sequence divergence between three bioinformatically identified UrsusERV proviruses. PCR products of expected sizes were amplified from polar bear, Syrian brown bear and black bear, but not from spectacled bear and giant panda ( Figure 3). Sanger sequencing of the PCR products confirmed that the target gag region of UrsusERV was amplified ( Figure S2). This approach further demonstrates that UrsusERV is absent from the panda and spectacled bear genomes but present in all examined members of the Ursinae subfamily of bears.
Journal 2015, volume, page-page Brown and black bears whole genome draft sequences were also screened for the proviral insertion sites identified in polar bears. While the proviral flanking sequences could be identified in the black bear genome, they lacked a proviral integration e.g. the flanking sequences were uninterrupted by UrsusERV. Therefore, while black bears are positive for UrsusERV sequences, the virus has integrated in distinct locations in the genome relative to polar bears. This could in part explain the differences in age estimates between the polar bear derived ERVs and the divergence time of black bear and brown bear/polar bear lineages in addition to the overall lower percentage identity of UrsusERV derived from black bears relative to those found in brown or polar bears. Out of the four female brown bears examined, two had the same UrsusERV insertion sites as polar bears for all the proviruses except for the provirus located in scaffold 200, which was polar bear specific as all brown bears tested lacked an integration at this genomic location although the flanking sequences were identified (Table S3). Two additional brown bear genome sequences, from Kenai and Admiralty Island, where found to have the same provirus insertions found on polar bear scaffolds 7 and 162, while the polar bear scaffold 1 integration was absent indicating insertion site polymorphisms in brown bears for the scaffold 1 proviral integration (Table S3) [51].
Molecular screening for UrsusERV sequences was also performed by sequence-specific PCR using genomic DNA extracted from a giant panda, spectacled bear (Tremarctos ornatus), Syrian brown bear (Ursus arctos syriacus), black bear, and polar bear to further verify the bioinformatic analysis. Genomic DNA was subjected to PCRs using primers designed to amplify a 312 bp region within the gag gene region showing little sequence divergence between three bioinformatically identified UrsusERV proviruses. PCR products of expected sizes were amplified from polar bear, Syrian brown bear and black bear, but not from spectacled bear and giant panda ( Figure 3). Sanger sequencing of the PCR products confirmed that the target gag region of UrsusERV was amplified ( Figure S2). This approach further demonstrates that UrsusERV is absent from the panda and spectacled bear genomes but present in all examined members of the Ursinae subfamily of bears. Figure 3. Gel electrophoretic separation of PCR products amplified from five bear species. PCR products and a 100 bp DNA ladder were separated through a 1.5% agarose gel. As indicated by a PCR product of an expected size of 312 bp, UrsusERV sequences were present in the examined polar, brown and black bears but absent from spectacled bears and the giant panda. Relevant sized DNA marker bands are indicated on the left.

UrsusERV Phylogenetic Analysis
We next examined phylogenetic relationships of UrsusERV with other retroviruses. Protein consensus sequences were generated for the GAG, POL, and ENV genes using the Retrotector putein predictions. The protein consensus sequences were searched against the NCBI protein database and sequences with identity >35% were extracted. The UrsusERV consensus, individual scaffold protein sequences, and protein sequences extracted from NCBI protein database were aligned using MAFFT [52]. Resulting protein alignments included murine, primate, koala, bat, avian, pig, and feline gammaretroviruses. Galidia ERV protein sequences, a basal class 1 gammaretroviral group, were used as an outgroup to perform Bayesian phylogenetic analysis (Figure 4) [6]. GAG and POL trees 7 Figure 3. Gel electrophoretic separation of PCR products amplified from five bear species. PCR products and a 100 bp DNA ladder were separated through a 1.5% agarose gel. As indicated by a PCR product of an expected size of 312 bp, UrsusERV sequences were present in the examined polar, brown and black bears but absent from spectacled bears and the giant panda. Relevant sized DNA marker bands are indicated on the left.

UrsusERV Phylogenetic Analysis
We next examined phylogenetic relationships of UrsusERV with other retroviruses. Protein consensus sequences were generated for the GAG, POL, and ENV genes using the Retrotector putein predictions. The protein consensus sequences were searched against the NCBI protein database and sequences with identity >35% were extracted. The UrsusERV consensus, individual scaffold protein sequences, and protein sequences extracted from NCBI protein database were aligned using MAFFT [52]. Resulting protein alignments included murine, primate, koala, bat, avian, pig, and feline gammaretroviruses. Galidia ERV protein sequences, a basal class 1 gammaretroviral group, were used as an outgroup to perform Bayesian phylogenetic analysis (Figure 4) [6]. GAG and POL trees displayed almost identical topologies with UrsusERV and the gibbon ape, koala, swine, bat, killer and whale retroviruses forming a clade sister to the murine and feline retroviruses clade. UrsusERV demonstrated greatest affinity with the PERVs for gag and was a sister clade to PERVs, gibbon ape leukemia viruses (GALVs), KoRVs and Mus caroli endogenous retrovirs (McERV)/Mus dunni endogenous retrovirus (MDERV) for pol and env (Figure 4). Among the individual proviruses, the scaffold 1 and 200 sequences representing the two youngest viral integrations occupied derived positions in the GAG and POL puteins trees. The ERV with the oldest integration date (scaffold 7) and the recombinant ERV on scaffold 162 occupied basal positions. This was not the case for the ENV putein tree where most scaffold sequence branching patterns could not be resolved and the youngest ERV (scaffold 200) occupied a basal position in the UrsusERV clade. However, with the exception of scaffold 200, much or all of the ENV protein was absent complicating phylogenetic analysis. Black and brown bear UrsusERV consensus sequence were constructed using the polar bear UrsusERV consensus sequence as reference. Phylogenetic analysis based on nucleotide sequences was also performed, where alignable, the black bear UrsusERV sequence was basal to the UrsusERV clade, while the brown bear UrsusERV sequence formed a sister clade to the polar bear sequence. Comparison of the UrsusERV clade with the other retroviral nucleotide sequences produced identical results to tree topologies from protein sequences. (Figures S3 and S4). The relative ages of bear lineages with respect to estimated bear lineage age, geological epoch and UrsusERV invasion events are shown in Figure 5.  (Figure 4). Among the individual proviruses, the scaffold 1 and 200 sequences representing the two youngest viral integrations occupied derived positions in the GAG and POL puteins trees. The ERV with the oldest integration date (scaffold 7) and the recombinant ERV on scaffold 162 occupied basal positions. This was not the case for the ENV putein tree where most scaffold sequence branching patterns could not be resolved and the youngest ERV (scaffold 200) occupied a basal position in the UrsusERV clade. However, with the exception of scaffold 200, much or all of the ENV protein was absent complicating phylogenetic analysis. Black and brown bear UrsusERV consensus sequence were constructed using the polar bear UrsusERV consensus sequence as reference. Phylogenetic analysis based on nucleotide sequences was also performed, where alignable, the black bear UrsusERV sequence was basal to the UrsusERV clade, while the brown bear UrsusERV sequence formed a sister clade to the polar bear sequence. Comparison of the UrsusERV clade with the other retroviral nucleotide sequences produced identical results to tree topologies from protein sequences. (Figures S3 and S4). The relative ages of bear lineages with respect to estimated bear lineage age, geological epoch and UrsusERV invasion events are shown in Figure 5.  Posterior probabilities >50% are shown. All sequences and analysis description are included in the material and methods section. UrsusERV (highlighted red) consensus and proviral sequences form a distinct clade that in all three analyses is closely related to PERV sequences. Figure 5. Maximum clade probability tree from mtDNA and nuclear DNA alignment display the Beast chronogram estimation. Clades age was estimated using fossil data points, with log normal distribution, and GTR model. Median divergences ages are shown above the blue horizontal bars that indicate the 95% highest posterior intervals of the estimation. Green arrow indicates the UrsusERV first insertion into the Ursinae clade that affected most likely all members of the clade. Red arrow indicates the second UrsusERV insertion into the brown bear clade, while the blue arrow represents scaffold 200 provirus that appears to be polar bear specific.
Comparison of host and retroviral phylogenies can reveal discordances that provide evidence for cross species transmissions of retroviruses. If an ERV and its hosts have co-evolved, the species tree and retroviral tree should be largely concordant. Discordant trees then indicate retroviral introgression among lineages or independent infection of different lineages by the same or related retroviruses [53]. Comparison of the host and retroviral phylogenetic relationships indicated substantial cross species transmissions in all lineages. UrsusERVs identified in the bears do not exhibit a consistent host pathogen co evolutionary pattern consistent with multiple invasions of the bear lineage from an unknown reservoir ( Figure 6). Comparison of host and retroviral phylogenies can reveal discordances that provide evidence for cross species transmissions of retroviruses. If an ERV and its hosts have co-evolved, the species tree and retroviral tree should be largely concordant. Discordant trees then indicate retroviral introgression among lineages or independent infection of different lineages by the same or related retroviruses [53]. Comparison of the host and retroviral phylogenetic relationships indicated substantial cross species transmissions in all lineages. UrsusERVs identified in the bears do not exhibit a consistent host pathogen co evolutionary pattern consistent with multiple invasions of the bear lineage from an unknown reservoir ( Figure 6). The host tree on the left was based on mtDNA cytochrome b gene sequences, and the retroviral tree on the right was based on the polymerase gene nucleotide sequence. Ursinae species that were positive for UrsuERV and the corresponding ERV sequences are shown in red. Phylogenetic trees were generated using maximum likelihood analysis as implemented in RAxML [54]. The bootstrap consensus trees illustrated were inferred from 500 replicates with the percentage bootstrap given next to each branch. Evolutionary phylogeny comparison of host versus retrovirus illustrates lack of co-evolution and cross species transmission events.

Discussion
Next generation virome sequence analysis from the brains of two polar bears yielded sequences with similarity to PERV and MDEV. The identified short sequences were used to screen bear genome sequence data, eventually identifying full-length proviruses and solo LTRs belonging to a previously undescribed gammaretrovirus group of ERVs that is specific to the Ursinae subfamily of bears, thus designated UrsusERV. Prior to the advent of genomics, the effort required to fully characterize novel multi-copy retroviruses was substantial. The increasing availability of partially or fully annotated genomes from wildlife is facilitating the characterization of viral sequences by providing a reference for obtaining full length viral genomes based on fragmentary sequences identified by either standard or NGS based approaches. The age estimates for UrsusERV suggests the group, like some PERVs, enJSRVs and KoRV, is overall relatively young, which is also reflected by the relative intact nature of the provirus on scaffold 200 and potentially by the unusual low copy number of elements identified [55]. Genomic invasion occurred subsequent to the divergence of the bear lineages such that basal bear lineages that include giant pandas and spectacled bears are devoid of UrsusERV homologues. Dating of scaffold 162 was not possible due to the recombinant nature of the 3′ LTR, yet the non-recombined regions of the proviral LTRs, that are identical in sequence, indicate a very recent integration of that provirus as well. The LTRs of the provirus in scaffold 200 were identical in sequence, suggesting it is an evolutionarily very young ERV locus. In fact, the integration in scaffold 200 was polar bear specific. Although young, unlike KoRV and some enJSRVs, all UrsusERV loci are likely fully endogenized as the integration sites were fixed, among independent modern polar bear whole genome data examined (n = 6). Fewer were identified in draft genome sequence from an ancient polar bear. However, this more likely reflects the lower coverage of the ancient genome sequence data than actual absence. 10 Figure 6. Tanglegram illustration of the phylogenies of gammaretroviruses and their hosts. The host tree on the left was based on mtDNA cytochrome b gene sequences, and the retroviral tree on the right was based on the polymerase gene nucleotide sequence. Ursinae species that were positive for UrsuERV and the corresponding ERV sequences are shown in red. Phylogenetic trees were generated using maximum likelihood analysis as implemented in RAxML [54]. The bootstrap consensus trees illustrated were inferred from 500 replicates with the percentage bootstrap given next to each branch. Evolutionary phylogeny comparison of host versus retrovirus illustrates lack of co-evolution and cross species transmission events.

Discussion
Next generation virome sequence analysis from the brains of two polar bears yielded sequences with similarity to PERV and MDEV. The identified short sequences were used to screen bear genome sequence data, eventually identifying full-length proviruses and solo LTRs belonging to a previously undescribed gammaretrovirus group of ERVs that is specific to the Ursinae subfamily of bears, thus designated UrsusERV. Prior to the advent of genomics, the effort required to fully characterize novel multi-copy retroviruses was substantial. The increasing availability of partially or fully annotated genomes from wildlife is facilitating the characterization of viral sequences by providing a reference for obtaining full length viral genomes based on fragmentary sequences identified by either standard or NGS based approaches. The age estimates for UrsusERV suggests the group, like some PERVs, enJSRVs and KoRV, is overall relatively young, which is also reflected by the relative intact nature of the provirus on scaffold 200 and potentially by the unusual low copy number of elements identified [55]. Genomic invasion occurred subsequent to the divergence of the bear lineages such that basal bear lineages that include giant pandas and spectacled bears are devoid of UrsusERV homologues. Dating of scaffold 162 was not possible due to the recombinant nature of the 3 1 LTR, yet the non-recombined regions of the proviral LTRs, that are identical in sequence, indicate a very recent integration of that provirus as well. The LTRs of the provirus in scaffold 200 were identical in sequence, suggesting it is an evolutionarily very young ERV locus. In fact, the integration in scaffold 200 was polar bear specific. Although young, unlike KoRV and some enJSRVs, all UrsusERV loci are likely fully endogenized as the integration sites were fixed, among independent modern polar bear whole genome data examined (n = 6). Fewer were identified in draft genome sequence from an ancient polar bear. However, this more likely reflects the lower coverage of the ancient genome sequence data than actual absence.
Black bears carry sequences homologous to UrsusERV but none of the integration sites are orthologous. An UrsusERV nucleotide consensus sequence that was generated from the black Viruses 2015, 7, 6089-6107 bear whole genome data shared 93.5% identity to the polar bear UrsusERV consensus sequence. Phylogenetic analysis of the nucleotide sequences indicated that black bear consensus is basal to the UrsusERV clade ( Figure S4). The basal position of the sequence and the lack of shared integration sites with brown and polar bears suggests that invasion of the black bear genome occurred independently from other ursids tested. However, we cannot rule out that ancient lineage sorting may explain the lack of orthologous integration among different bear species. All brown bears shared the proviral insertion in scaffold 162 consistent with the predicted date of divergence of polar bears and brown bears ca. one million years ago [50]. The provirus in scaffold 7, the oldest UrsusERV locus dated at approximately 2.1 million years, was present, based on draft sequence evidence, in the genomes of three out of four brown bears. The proviral integration on scaffold 1 was identified in two out of four brown bears. Sex chromosomes were excluded as the source of the polymorphism observed, as all brown bears examined were female. The polymorphic nature of these integrations in brown bears could indicate that the ancestor of brown bears was polymorphic for these insertions and different brown bear populations retained or lost the ERVs respectively. In contrast, according to our results, the provirus on polar bear scaffold 200 was absent from all brown bears suggesting invasion happened very recently and is exclusive to the polar bear lineage. The similarity to PERVs, the relative young age and the lack of congruence of viral homology and host specificity suggests that the exogenous UrsusERV equivalent was repeatedly in circulation over at least the last two million years from an unknown reservoir. The general lack of host and viral phylogenetic congruence observed suggests that UrsusERV related viruses have been independently infecting multiple mammalian lineages over the last 12 million years ( Figure 6). Although generally basal in the phylogenetic trees, the killer whale endogenous retrovirus and UrsusERV show relatively strong affinity suggesting that a common reservoir for UrsusERV like viruses are present in the Arctic region. The very young nature of the provirus on scaffold 200 further suggests that exogenous forms could potentially still be in circulation.
The identification of unique LTR groups associated with different proviruses suggests there were invasions of at least two distinct UrsusERV variants. The phylogenetic positions of the scaffold 7 and 162 proviruses, which are basal in both the GAG and POL trees, suggest LTR-B containing viruses invaded the bear genome first. A subsequent invasion of LTR-A type variants represented by proviruses on scaffolds 1 and 200 followed, which have younger estimated integration ages and are more derived in the phylogenetic analysis. This may reflect a shift in viral strains circulating in the source species. The independence of the invasions is also suggested by the absence of orthologous integrations in the black bears for the polar bear scaffolds. Were the black bear loci orthologous with the polar bear loci, the age estimates would be vast underestimates given that the black bear lineage diverged from the polar bear and brown bear lineages at least 7 mya ( Figure 5). If the invasions were independent as it seems, then the polar bear and brown bear UrsusERV younger age estimates are more consistent with divergence date estimates for the brown and polar bear clades. Confounding this analysis is the basal position of the provirus on scaffold 200 in the ENV phylogeny. The provirus on scaffold 1 could not be analyzed for ENV as the gene is deleted. This incongruity could represent recombination post integration or that the ERV is under selective pressure with the env gene having accumulated mutations [56]. An alternative explanation may be an assembly artifact as the scaffold 1 provirus harbors an assembly gap within the pol gene region. Even though independent invasions by LTR-A and LTR-B viruses is the most parsimonious explanation for the observed proviral distribution, the possibility of a single endogenization event, diversification and subsequent lineage sorting cannot be excluded as an explanation based on the current data.
UrsusERV occupied in all phylogenetic analyses a sister taxa position relative to a clade containing the PERVs, MDEV, McERV, GALVs and KoRVs. Only the ENV tree conflicted with this result placing McERV and MDEV in a basal position, which may reflect recombination events in those viruses. A consistent result across all analyses was the basal position of the killer whale endogenous retrovirus, which has an estimated endogenization date in the delphinoid genomes sometime in the past 12 million years [1].
With a young age estimate for PERVs, KoRVs and UrsusERV, and with GALV being only identified as an exogenous virus, the entire clade is relatively young in comparison to many other retroviral groups [47,[57][58][59]. This suggests retroviruses with broad host range were disseminated across a wide variety of environments and species in the last 12 million years and, in the case of UrsusERV, it remains possible that exogenous forms related to the scaffold 200 provirus are still in circulation. As more next generation whole genome sequence data from wildlife becomes available and more endogenous retroviral genomes are characterized, the origins and distribution of gammaretroviruses such as UrsusERV will become clearer.

Samples
Polar bear tissue samples from Knut, a male polar bear, were provided by the Zoological Garden Berlin. Tissue samples from Jerka, a female polar bear, a spectacled bear (Tremarctos ornatus) and a black bear (Ursus americanus) were provided by the Zoological Garden Wuppertal, Tierpark Berlin, and Allwetterzoo Münster, respectively. All tissue samples were frozen at´80˝C until further processing.

Nucleic Acid Preparation and Next Generation Sequencing
DNA and RNA extraction from brain and liver bear tissue samples was performed with approximately 25 mg of tissue using a QIAamp DNA mini extraction kit and RNeasy Lipid kits following manufacturer instructions respectively. RNA sample preparation and next generation sequencing was performed as described in [23]. DNA and transcriptome data generated in [23] were used for data mining. Retroviral sequence data mined from the transcriptomes has been deposited in NCBI Sequence Read Archive (SRA) with accession number SRP065872.

Mining of UrsusERV Sequences in the Polar Bear Genome
Viral database blast searches of the next generation sequence data of the two polar bear transcriptomes were performed as previously described [4]. Viral database screening results revealed reads with similarity to the Mus dunni endogenous retrovirus, Feline sarcoma virus, Murine leukemia virus, and to the Porcine endogenous retroviruses.

Generation of UrsusERV Provirus and LTR Consensus Sequences
Multiple alignments of the proviral amino acid sequence and LTR's nucleotide sequence were generated using MAFFT v7.017 FFT-NS-i alignment strategy [52] as a plugin within the Geneious (Biomatters, Auckland, New Zealand) software. Resulting alignments were manually curated Viruses 2015, 7, 6089-6107 and a majority rule consensus was generated for each of the proviral genes as well as the LTRs. Retrotector [37] was employed to generate putative, reconstructed retroviral proteins, called puteins, for Gag, Pro, Pol, and Env genes based on the proviral consensus sequence. Conserved motifs in resulting consensus putein sequences were further verified using the NCBI Conserved Domain Database [38]. Proviral and LTR aligned sequences can be found in Supplementary fasta files 1 and S2. The UrsusERV consensus sequences can be found in RepBase under UrsusERV-A and B for the two LTR variants identified respectively and UrsusERV-I for the non-LTR sequences.

Recombination Inference Analysis
Recombination analysis was performed for multiple sequence alignments using three computational methods: Recombination Analysis Tool (RAT, Norwich, UK), Genetic Algorithm Recombination Detection (GARD), DualBrothers and RECCO [39][40][41][42][43]. GARD was used as implemented in the datamonkey online server [60] while DualBrothers (dual multiple change model) was used as a plugin in the Geneious software. RAT and GARD were used with default parameters, while DualBrothers was used with default model prior parameters and the following scanning window parameters: Window length = 300, Step size = 10. RECCO was implemented using the default parameters.

Sequence Read Archive and Bioinformatic Analysis
NCBI Sequence Read Archive data were used to determine whether the proviral sequences identified in the polar bear genome and the retroviral reads identified in the transcriptome were present in other bear species and whether the same proviral flanking regions were present in independent polar bear samples. Whole genome sequence data of Ursus americanus (accession number: SRR830685), Ursus arctos (SRR830213, SRR1693654, SRR1693624, SRR518712), Ursus maritimus (SRR827584-5, SRR827600, SRR947747, SRR942297), and ancient Poolepynten Ursus maritimus (SRR518649, SRR518651-7, SRR518704-5) were obtained. Short reads obtained were first converted to fastq format using sratoolkit version 2.3.5 [27]. Resulting fastq files were trimmed for adaptors and read quality using cutadapt version 1.5 and then were aligned to the draft polar bear genome using BWA 0.7.9a-r786 mem algorithm [49,61]. Resulting alignments were further curated using samtools 0.1.19 [62]. The identified proviral sequences from the draft polar bear genome were also used to search the Panda genome (Ailuropoda melanoleuca) by BLAT search at the UCSC Genome browser [48].

Phylogenetic Analysis
Multiple amino acid sequence alignments were generated using MAFFT version 7 and the iterative refinement method [52].
Further curation of the resulting alignments was performed manually.
Phylogenetic analysis was performed using consensus amino acid sequences of the three major proteins (Gag, Pol and Env) of UrsusERV and related protein sequences obtained from GenBank: BaEV (accession numbers: BAP91048.1, BAA89659. Galidia ERV sequences were used as outgroup in the analysis. Protest version 3.4 was used to identify the best model for the phylogenetic analysis [63]. A Jones model with gamma distribution and invariable sites was determined as the optimal model for GAG and POL protein alignments, while for the ENV amino acid alignment the Rtrev model with invariable sites and gamma distributions was chosen. Bayesian inference analysis was performed as described in [64].

Age Estimation of UrsusERV Proviral Sequences
To estimate UrsusERV retroviral integration ages into the polar bear genome we used a molecular clock method. Proviral LTRs, upon integration into the host genome, are identical in sequence and they then acquire mutations independently from each other over time at the mutation rate of the host DNA [20]. Sequence divergence between the 5 1 and 3 1 LTRs of a provirus can thus serve as a molecular clock. Hypermutable CpGs sites were excluded from the age estimations. Sequence divergences were subsequently corrected using the Kimura-2-parameter (K2P) model [65]. The provirus formation dates were estimated using the formula T = D/2* 0.0015 per nucleotide and million years, where D is the corrected sequence divergence of the proviral 5 1 and 3 1 LTRs and 0.0015/nt/myr is the reported mutation rate in polar bears [25].