Comprehensive profiling of retroviral integration sites using target enrichment methods from historical koala samples without an assembled reference genome

Samples and ancient DNA extraction

A total of thirteen museum samples were examined (Table 1). DNA extractions were performed in the ancient DNA (aDNA) laboratory of the Department of Wildlife Diseases of the Leibniz Institute for Zoo and Wildlife Research in Berlin, Germany. The laboratory is dedicated to aDNA work and has never been used for molecular work on modern samples. The room is UV irradiated 4 h every night by ceiling-mounted UV lights. All work performed in the facility follows procedures designed to minimize the possibility of contamination, such as the use of laminar flow hoods and use of protective clothing. The samples used in this study were all derived from museum skin samples and thus no living koalas were sampled at any point during this study.

DNA from approximately 250 mg of koala skin tissue (7 mm × 7 mm) per museum specimen was obtained using a silica-based extraction kit for aDNA (GENECLEAN Ancient DNA Extraction Kit, MP Biomedicals, USA). The protocol followed the manufacturer’s instructions and has been successfully applied to a variety of ancient sample types (Wyatt et al., 2008; Roca et al., 2009). Mock extractions were performed with each set of koala museum specimens as negative controls for extraction. Subsequent to each extraction, the isolated DNA was further purified using a MinElute spin column (Qiagen, Hilden, Germany) as described in (Gilbert et al., 2007) to remove potential inhibitors for the subsequent enzymatic reactions. DNA extracts were not quantified because of the small proportion of endogenous DNA compared to exogenous DNA (contaminants such as bacteria, fungi) in typical aDNA samples.

NGS Library preparation

Illumina sequencing libraries were prepared from the extracts using a previously described protocol (Meyer & Kircher, 2010) with the following modifications: (A) All SPRI purification steps were substituted with spin column purification (MinElute PCR purification kit, Qiagen). (B) Adapter concentration in the ligation reaction was reduced to 0.2 mM per adapter. (C) The purification after adapter fill-in was substituted by heat inactivation at 80 °C for 20 min. The libraries were then used directly as a template for subsequent amplification following a two-step strategy, as previously described (Kircher, Sawyer & Meyer, 2012). A quality control strategy (Meyer et al., 2008) was also applied, which consisted of qPCR to quantify the product after each step of library amplification. Based on qPCR results, three samples for which DNA quality was too poor for analysis were excluded from further processing.

In the first round of amplification, AmpliTaq Gold, a non-proof reading enzyme, and indexing primers (Table S1) were applied, adding a distinct P7 index to each library as described in Meyer & Kircher (2010), 10 indices for the 10 working samples and 3 and 4 negative control indices for PEC and SPEX respectively. Adding distinct indices to each library allows for multiple samples to be sequenced in a single sequencing run. The non-proof reading enzyme allows for amplification to be performed on templates containing deoxyuracils, which are common with aDNA (Der Sarkissian et al., 2015). After removal of 1 µL for qPCR as a library quality control, the libraries were used as template in 100 µL PCR containing 1x Taq buffer II (Applied Biosystems), 5U AmpliTaq Gold (Applied Biosystems), 250 mM each dNTP and 100 nM each indexing primer. Cycling conditions followed manufacturer’s instructions: the pre-denaturation step lasted 12 min at 95 °C, followed by 12 cycles of denaturation at 95 °C for 20 s, annealing at 60 °C for 30 s and elongation at 72 °C for 40 s, with a final extension step of 72 °C for 5 min. PCR products were purified using the QIAquick PCR purification kit (Qiagen, Hilden, Germany).

In the second round of amplification, 5 µL of the purified PCR product from the first round PCR was used as a template for a second PCR. This involved 50 µL reactions containing Herculase II Fusion DNA Polymerase (Agilent Technologies Catalog 600677), which has proof reading activity, and primers IS5 and IS6 (Meyer & Kircher, 2010) at a final concentration of 400 nM each. Cycling conditions included an activation step of 3 min at 95 °C, followed by 15–20 cycles of denaturation at 95 °C for 20 s, annealing at 60 °C for 25 s and elongation at 72 °C for 30 s, with a final extension step at 72 °C for 3 min. The number of cycles used in the PCR for every sample was dependent on the concentration of each of the libraries as determined by the qPCR assay. The PCR amplified libraries were then purified using the QIAquick PCR purification kit. Each library was separately used in subsequent PEC and hybridization capture experiments.

Bait preparation and integration site enrichment

Three methods were compared for retrieving integration sites: single primer extension (SPEX) primer extension capture (PEC), and hybridization capture. All three have been successfully applied to ancient and historical DNA samples and all are applicable to samples that would not be expected to yield results with conventional methods for integration site analysis. The same set of primers was used in PEC and SPEX experiments (Fig. 2, Table S2). Because the two LTRs of a provirus are identical, the primers designed for enriching the 5′ integrations will also target the 3′ LTR retroviral env gene and the primers designed for targeting the 3′ integrations will also extend targeting the retroviral gag leader sequence (Fig. 2A). For both the 5′ and 3′ KoRV LTR, two 20 bp primers were developed which overlap such that the 3′ end of the first primer overlapped 8 bp with the 5′ end of the second primer (Fig. 2B, primers 5.1 and 5.2 and 3.1 and 3.2 respectively). To avoid known LTR polymorphisms among KoRV proviruses, the two primers on each side of the LTR were located 17 bp from the 5′ end and 50 bp from the 3′ end of the LTRs in conserved regions (Fig. 2B). The baits used for hybridization capture were synthesized to generate 32 bp oligonucleotides that spanned the full length of sequence covered by primers 5.1 and 5.2 (32 bp) on the 5′ LTR and primers 3.1 and 3.2 (32 bp) on the 3′ end.

Primer Extension Capture (PEC)

Indexed libraries were pooled in equimolar ratios for primer extension following a published protocol (Briggs et al., 2009). After each step, 1 µL of the product was quantified by qPCR. To minimize the amplification bias, each of the captured products was amplified in triplicate, using 5 µL of the captured product as template for each reaction, using the same kit and cycling conditions as described previously under NGS library preparation for second round amplification of Illumina indexed libraries, except that we ran 20 cycles of amplification for all samples. Amplified captured libraries were purified using the QIAquick PCR purification kit (Qiagen, Hilden, Germany) and eluted in 50 µL of elution buffer (EB) and used as a template for a second round of PEC.

Single Primer Extension (SPEX)

The SPEX experiments generally followed a published protocol (Brotherton et al., 2007) using DNA extracts prior to Illumina library construction with three modifications: (1) Illumina sequencing adaptors were attached to the 5′ end of the primers used in the first round of partially nested PCR; (2) MyTaq HS Mix (Bioline, BIO-25045) was used instead of Platinum Taq DNA Polymerase High Fidelity in the first round of a partially nested PCR; (3) only one round of a partially nested PCR amplification was performed. The nested PCR products were then quantified by qPCR and indexed using Illumina indexing primers (Table S3). The indexed PCR products were purified using a QIAquick PCR Purification Kit (Qiagen). The amplicons were quantified by qPCR and subjected to a second round of amplification using the same conditions as the first round. The products were purified again using the QIAquick PCR Purification Kit (Qiagen), quantified by qPCR and pooled at equimolar ratios. All PEC and SPEX products were pooled and measured using High Sensitivity DNA chips on an Agilent 2100 Bioanalyzer, then sequenced at the National High-throughput DNA Sequencing Centre, Copenhagen, Denmark using Illumina MiSeq Reagent Kit v2 (300 cycle).

Hybridization capture

The amplified libraries were pooled in equimolar ratios, totaling a final amount of 2 µg. An established protocol was followed (Maricic, Whitten & Pääbo, 2010) except that synthesized oligonucleotide baits were used instead of PCR products and the EB volume for final elution using Qiagen MinElute column was 20 µL instead of 15 µL. After 2 days of hybridization and subsequent elution steps, 1 µL of the final eluate was quantified by qPCR and 5 µL (in total 15 µL) was amplified in triplicate using the same kit and cycling conditions as described in the NGS library preparation for second round amplification of Illumina indexed libraries. The pooled PCR products were purified using the QIAquick PCR Purification Kit and was measured using the Tapestation 2200 (Agilent Technologies Catalog G2964AA). Hybridization capture libraries were sequenced at the National High-throughput DNA Sequencing Centre, Copenhagen, Denmark using Illumina MiSeq Reagent Kit v2 (300 cycles).

Preprocessing of sequence data

Adaptor sequences were removed from sequence reads using cutadapt-1.2.1 (Martin, 2011), and quality trimming was performed using Trimmomatic-0.22 with default settings (Bolger, Lohse & Usadel, 2014). The paired forward and reverse sequence reads were merged using Flash-1.2.5 where possible (Magoč & Salzberg, 2011), and both the merged and unmerged reads were used for further analysis. PCR duplicates (clonality in the sequencing data) with 100% sequence identity were removed using cd-hit-v4.6.1 (Li, Jaroszewski & Godzik, 2001).

Identification of KoRV integration sites

Figure 3 and Table S6 summarize the computational pipeline used for the identification of KoRV integration sites. For its implementation, both existing software and customized perl scripts were used that made use of BioPerl (Stajich et al., 2002). Because the nested primers or bait were designed near the ends of LTR, the primer extension products would include either the first 49 bp of the 5′ LTR or the last 82 bp of the 3′ LTR, which are designated “LTR ends” in Fig. 2A. All sequences with a KoRV flank should contain an LTR end, as a result of the primer extension (Fig. 2B). Therefore, KoRV integration sites could be identified as the sequence beyond the KoRV LTR end since all integration sites would be attached to an LTR sequence. However, due to DNA degradation in museum samples, some primer extension products may not have a complete LTR end. Furthermore, minor deletions at the end of the integrated LTRs may be present (Fields, Knipe & Howley, 1996); for example, a 19 bp deletion was found in a KoRV provirus (Ishida et al., 2015). To get around these potential issues, identification of the LTR ends relied on sequentially selecting sample sequences that contain defined LTR segments; this was done in separate steps for the 5′ and 3′ flank-containing sequences. The LTR end was divided into two segments, designated A and B (Fig. 2B): the B segment corresponds to the last 19 bp of the LTR and is referred to as 5B or 3B in the 5′ and 3′ LTR ends, respectively. The A segment is the remaining section of the LTR end, which has a length of 30 bp in the 5′ end (5A) and 63 bp in the 3′ end (3A).

Initially, sequences containing either of the two A regions in the KoRV LTR end (5A or 3A in Fig. 2B) were identified. For this step, optimal local pairwise sequence alignments (Smith-Waterman, EMBOSS Rice, Longden & Bleasby, 2000) were computed between each sample sequence and the A region in either the 5′ or 3′ LTR end. Sequences were used for further analysis if they could be aligned to at least 20 bp of the 30 bp 5A segment with at least 90% identity, or if they could be aligned to at least 43 of the 63 bp 3A segment with at least 90% identity (Table 2). Despite the differences in the lengths of the 5′ and 3′ A segments, this alignment criteria was selected as it resulted in approximately the same sequence identity threshold (∼60%) for both ends. Sequences not passing these criteria were discarded as artifacts. The LTR ends of all sequences meeting these criteria were trimmed to the distal 19 bp and then used for further analysis. A higher sequence identity threshold was not chosen due to potential DNAse degradation of the molecules or ancient DNA based damage lowering homology. A 20 bp sequence length was the minimum on the 5A segment that allowed for LTR identification, whereas the 3A region was longer allowing for a longer minimum segment.

From these sequences, B segments of either 3′ or 5′ LTR ends were identified (3B or 5B in Fig. 2B). For this step, optimal local sequence alignments were computed between each of the trimmed sequence and the B segment in either the 3′ or the 5′ LTR end. Only sequences that could be aligned to at least 12 bp of the 19 bp long B segment (3B or 5B) with at least 80% (Table 2) identity were selected. This criteria was chosen by considering the known polymorphisms in this region of KoRV that originate from the mutagenic properties associated with LTR/host flanking region junctions. The last 19 bp of LTR ends were trimmed from all sequences meeting the selection criteria, leaving LTR free KoRV flanks or KoRV genomic DNA adjacent to the LTR.

All sequences that contained the A region, but for which the B region was not detected using the pairwise alignment strategy, were then subjected to another test. Specifically, these sequences were used as queries for two separate local database searches using BLAST (Altschul et al., 1990). Such sequences represent LTRs that have suffered deletions at the end, a common occurrence in proviruses. One search was against HiSeq sequencing data of a koala from Queensland, Australia with 100X coverage. The data represent raw Illumina sequences and are not annotated or assembled. After adaptor and quality trimming, 6.469 billion reads from this koala, with a mean length of 78 bp, were used for this step. Sequences were considered KoRV integration sites when their non-LTR portion could be aligned with greater than 90% identity to the koala reads over 60% length of the sample sequence. A second search was against the Tammar wallaby (Macropus eugenii) genome (GenBank: ABQO000000000.2), which represents the closest related species to koala for which a genome has been assembled (Renfree et al., 2011). Although the wallaby and koala lineages diverged more than 50 Mya (Meredith, Westerman & Springer, 2009), we expected that some of the koala genomic DNA (flanking KoRV) could be aligned to the homologous wallaby regions. Sequences with at least 70% identity over 50% length of the sample sequence to the wallaby genome were therefore considered to be KoRV integration sites. An assumption made is that KoRV does not frequently or preferentially insert into repetitive sites which could cause us to underestimate the total number of integrations. This will only be resolvable once an annotated koala genome becomes available. For the sequences with a match to the wallaby scaffolds or the koala data, the LTR sequences were trimmed and were then concatenated with the KoRV flanks (obtained in previous steps) for further analysis.

Sorting of sequences representing different integration sites

All sequences with matches to the different segments of the 3′ and 5′ LTR ends and/or to wallaby scaffolds or koala HiSeq data from each of the enrichment techniques were collected. The sequences matching 3′ and 5′ LTR ends were kept separate, resulting in a total of six different data sets for further analysis (two data sets each for the PEC, SPEX and hybridization capture). LTR ends were removed from all sequences in these data sets. Before using these sequences to identify shared and unique integration sites, all KoRV flanks were sorted into three categories by length (Table 3): (1) including KoRV flank sequences shorter than 4 bp, the typical length of a KoRV target site duplication. These sequences were valid and in the right extension direction but too short for any biological interpretation, and thus were excluded from further analysis. (2) KoRV flanks with length of 4–14 bp. (3) KoRV flank sequences with length of 15 bp or longer. Both KoRV flanks 4–14 bp or with length of 15 bp or greater were used for identifying shared and unique integration sites, but only KoRV flanks of minimum length of 15 bp were used for subsequent pairing of 5′ and 3′ integration sites to one KoRV provirus. Additionally, as per the experimental design (Fig. 2A), each of the 5′ and 3′ primer extension products has two directions of extension, a and b. Extension a is towards the KoRV flanks yielding integration sites as expected, while extension b is towards the KoRV proviral genome yielding unwanted products for this integration site study. These sequences were designated as ‘internal KoRV reads’. However, despite not representing integration sites, extension b products still represent correctly enriched products from the specific enrichment technique.

Combining the KoRV flank types 4–14 bp long and 15 bp or longer, the PEC data had 392 5′ flank sequences and 2,347 3′ flank sequences; the SPEX data 6,521 5′ flank sequences and 9,200 3′ flank sequences; and hybridization capture 1,158 5′ flank sequences and 28 3′ flank sequences. A clustering approach was used to sort all sequences in each of the six data sets into groups of similar sequences; each cluster representing a unique integration site. Sequences that did not share significant similarity with any other sequences in the input file were called singletons. For each of the six data sets, all-against-all BLAST comparisons were run, and the BLAST output was used as input for clustering using TRIBE-MCL (Enright, Van Dongen & Ouzounis, 2002), separately for each data set. Different combinations of E-values (all against all BLAST) and inflation values (TRIBE-MCL) were used for this step and the optimal parameter combination for each data set was evaluated. For all combinations of E-values and inflation values, multiple sequence alignments were computed for all clusters using MAFFT v7.127b (Katoh et al., 2002). To assess the quality of the clustering, alignments of the 30 largest clusters of each clustering result were visualized in jalview (Waterhouse et al., 2009) and were verified by eye. An alignment was considered high quality if the total number of mismatches and gaps in every sequence of the alignment was no more than 10% of the sequence length. If all 30 clusters were evaluated to be of high quality, the sequence was further analyzed. The parameter combinations for optimal clustering and related all against all BLAST are listed in Table 4.

Singletons and non-singleton clusters containing sequences derived from a single individual koala were considered to represent unique integration sites. Clusters containing sequences shared by more than one koala were considered to represent shared integration sites (Table S4 and S5). A consensus sequence was computed from the alignment of each sequence cluster. Singletons and consensus sequences were then further evaluated first by computing pairwise alignments between these sequences and the gag or env part of KoRV genome (Fig. 2A) (GenBank: AF151794.2). The sequences that could be aligned to the KoRV genes with at least 90% identity and of any length were categorized as primer extension or flank capture within the KoRV genome. The LTR sequences at the 5′ and 3′ ends of the KoRV genome are identical or nearly so and therefore 50% of the PCR products should extend into the KoRV genome (Fig. 2A). Sequences that could not be mapped to KoRV genome were potential KoRV integration sites and were evaluated further. For such sequences, a length filtering was performed with a threshold of 15 bp, since this is the minimum length that can be effectively identified by BLAST. The sequences longer than 15 bp were first used as query in BLAST to search against the koala shotgun Hiseq data; they were also mapped to wallaby genome (GenBank: ABQO000000000.2) in Geneious version 6.18 (http://www.geneious.com, Kearse et al., 2012). Identified sequences for either one of the two computations were considered to be KoRV integration sites. Sequences shorter than 15 bp are too short for efficient mapping or BLAST; however, because they contained an LTR end, were included in the KoRV specific enrichment statistics (Table 2), although they were not further analyzed.

Pairing of 5′ and 3′ integration sites to one KoRV provirus

Ishida et al. (2015) identified the length of the retroviral target site duplication (a stretch of host DNA directly adjacent to the retrovirus which is duplicated during retroviral integration) for KoRV to be 4 bp. Based on this target site duplication length (Fig. 4), all 5′ and 3′ integration sites were examined for shared four bp target site identity. Only KoRV flanks 15 bp or longer were used for pairing 5′ and 3′ integration sites. The minimum 26 bp (30 bp minus the four bp target site duplication) combined length discriminated true wallaby matches from non-significant blastn results.

The paired 5′-3′ flanking sequences were (1) mapped against the wallaby genome using the mapping tool in Geneious using default settings, where only the paired 5′-3′ integration sequence that could be mapped to the wallaby genome with over 70% of their total length were scored as positively identified; (2) used as query to search in the HiSeq data of Queensland wild koala using BLAST. Here, only the paired 5′-3′ integration sites that could be aligned with over 90% identity with the koala HiSeq reads were considered positive.

Statistical analysis of shared integration sites

Statistical tests were performed to check if the occurrences of KoRV at sampled integration sites increased as the samples became younger among the 10 museum koala samples. Two logistic regression models were employed: one for 5′ integration sites and one for 3′ integration sites. Both models had the same structure. The occurrence was considered (binary: 1 = presence, 0 = absence) as the response variable and time as a continuous fixed effect. Because results were qualitatively similar irrespective of expressing “time” as rank or directly as years, for the sake of simplicity, only the latter was reported. The identity of koalas and of insertion sites were considered as two Gaussian random effects, making this logistic regression a Generalised Mixed effect Model (GLMM). The GLMM was fitted using the function HLfit from the R package spaMM 1.4.1 (Rousset & Ferdy, 2014), considering a Binomial error structure. The effect of time was tested by performing an asymptotic Likelihood Ratio Test (LRT) using the function anova.HLfit from the same package.

Single primer extension (SPEX)

Using SPEX to target the 5′ LTR flanks, 66 integration sites unique to a single koala, and 15 integration sites shared by at least two koalas were identified across the 10 koala samples, for sample descriptions see Table 1. These integration sites derived from consensus sequences generated from sequence clusters with at least 4 bp of sequence (representing the length of the target site duplication of KoRV) (Ishida et al., 2015) flanking the KoRV LTR (categorized as either KoRV flank sequences 4–14 bp long or 15 bp or longer in Table 3). An additional 15,822 sequences were less than 4 bp; these could not be further analyzed since their length was shorter than the target site duplication, these are listed as KoRV flanks shorter than 4 bp in Table 3. Additionally, 212 reads were identified as part of the envelope gene of the KoRV genome. This results from the presence of identical primer target sites in the 5′ and 3′ LTRs (Fig. 2A), since KoRV 5′ and 3′ LTRs are identical or nearly so (Ishida et al., 2015). Thus, approximately 50% of the sequences are expected to be internal KoRV proviral reads that extend from the LTR into the proviral genome rather than into the host flanking region. For clarity, we term these sequences ‘internal KoRV reads’. These sequences that extended into the KoRV genome were categorized separately but included in the total enrichment efficiency evaluation because they still represent correct enrichment of target sequences. SPEX for integration sites next the 3′ LTR also identified 182 unique and 28 shared 3′ LTR flanks; with 1,527 sequences being too short for further analysis (less than 4 bp of flank sequence) and 223 internal KoRV reads that matched the KoRV genome (Table 3).

Primer extension capture (PEC)

PEC was designed to identify flanking regions 5′ of integration sites and detected 126 unique and 17 shared integration sites. An additional 496 sequences included less than 4 bp of flank that was too short for further analysis, while 141 internal KoRV reads extended into the KoRV genome. PEC targeting regions downstream of the 3′ LTR integration sites identified 538 unique and 134 shared integration sites. An additional 1,806 reads were less than 4 bp, while 1 internal KoRV read was identified that matched KoRV (Table 3).

Hybridization capture

Using the 5′ LTR region as bait, 862 unique and 25 shared 5′ flanking regions were identified by hybridization capture. An additional 191 sequences included less than 4 bp of flank and 151 internal KoRV reads were characterized. Using the 3′ LTR region as bait, only 24 unique and no shared integration sites were identified by hybridization capture. The strong bias of this method towards identifying 5′ integration sites has been previously observed (Tsangaras et al., 2014b). Additionally, 41 sequences included less than 4 bp of flank, while 14 sequences were classified as internal KoRV reads (Table 3).

Summary of computational data processing

At each step of our bioinformatics pipeline, we recorded for each experiment the number of sequences that met our screening criteria (Fig. 3). The mean length, minimum length and maximum length of sequences were also calculated at each step (Table S6). Before any screening criteria were applied, SPEX produced 7,628 million, PEC produced 6,956 million reads, and hybridization capture produced 31,096 million. After screening and PCR duplicate removal of this sequencing data, 9.37% of the initial sequencing reads were kept for SPEX, 17.08% for PEC, and 37.55% for hybridization capture. Clonal sequences i.e., duplicate sequences resulting from PCR bias in amplification were more prevalent for products of SPEX than for products of either PEC or hybridization capture.

Two rounds of bioinformatics LTR end identification was performed (Fig. 3). After the first round of LTR end identification, 142,577 (19.94% of the reads after pre-processing) LTR positive sequences were identified for SPEX, 31,787 (2.67% of the reads after pre-processing) for PEC, and 5,648 (0.05% of the reads after pre-processing) for hybridization capture. Sequences passing the second round of 5′ LTR end selection were 22,542 for SPEX, 1,029 for PEC, and 1,495 for hybridization capture, while the sequences passing the second round of 3′ LTR end selection were 10,950 for SPEX, 5,559 for PEC, and 83 for hybridization capture. No KoRV LTR ends were detected in negative controls, extraction or PCR controls lacking template, for any experiment.

Cross-technique comparisons

The efficiency of target enrichment for each technique was calculated as the total number of identified flank sequences divided by the total number of sequences after removal of clonal sequences. The total number of target enrichment products included KoRV flanking sequences of any length and internal KoRV reads.

As shown in Table 3, PEC enriched the highest total number of 3′ integration sites, 672, whereas hybridization capture enriched the most 5′ integration sites, 887. As a percentage of the total sequences retrieved, SPEX achieved the highest target enrichment efficiency (4.68%). Both PEC and hybridization capture exhibited lower enrichment percentages (0.55% and 0.01% respectively).

Due to a phenomenon known as CapFlank (Tsangaras et al., 2014a), koala genome sequences near the integration sites may be enriched together with KoRV flanks by concatenation of library molecules on the baits. To estimate the numbers of such target flanks, after PCR clonal sequence removal, all sequences were screened using BLAST against the wallaby, which represents the phylogenetically closest species to the koala with an assembled genome. Hybridization capture exhibited the lowest efficiency of on-target enrichment (0.01%, Table 3) and highest ratio of CapFlank enrichment (16.41%), while SPEX achieved the highest efficiency of on-target enrichment (4.68%) and lowest ratio of CapFlank enrichment (0.23%).

As illustrated in Fig. 5, for the 5′ LTR integration sites, hybridization capture yielded the highest total number of integration sites, 887, and contained 91.36% of the integration sites identified in the SPEX data set and 86.71% of the integration sites identified in PEC data set. The 3′ LTR integration data followed a different profile with PEC generating the highest total number of integration sites, 692, containing 85.07% of the integration sites in the SPEX data set and 91.67% of the integration sites in the hybridization capture data set. The expected enrichment profile and approximate location of the recovered reads based on each of the three methods performance is shown in Fig. 6.

Shared and unique integration sites

After identical integration sites across the data sets generated by the 3 techniques were combined, 52 shared and 865 unique 5′ KoRV host flanks could be identified. Shared integration sites accounted for 5.7% of the total number identified using 5′ flanking host sequences, a similar percentage as estimated in a previous study (Tsangaras et al., 2014b). Among the 3′ flanking regions, 146 shared and 570 unique integration sites were identified, with shared sites accounting for 20.4% of total integration sites identified using 3′ host genomic sequences.

Pairing of 5′ and 3′ flanking regions to identify individual proviral integration sites

KoRV typically produces a 4 bp target site duplication upstream and downstream of its integration site (Ishida et al., 2015). All 4 bp putative 5′ target site duplications were compared to all 4 bp putative 3′ duplications. In cases where there was an exact match at the 4 bp, the two flanking regions were concatenated to simulate the sequence that a virus-free host would have at that locus, assuming that the target site duplications actually were from the same locus. A total of 1,690 concatenated 5′ and 3′ host flanking sequences were used to query the koala HiSeq data to identify proviral integration sites (Fig. 4). There were 63 matches, indicating that the 5′ and 3′ flanks actually corresponded to integration sites at the same KoRV proviral locus. Of these 63 loci, 40 corresponded to proviral integration sites present in a single koala (Data S1), whereas 23 corresponded to a proviral integration site detected in at least two koalas (Data S2).

The comparison of integration sites across different studies

The KoRV integration sites identified by this study were compared to those reported by Tsangaras et al. (2014b) and Ishida et al. (2015) (Table S3). Each study used a different set of koalas, and there was no overlap in koala individuals examined by the three studies. All but one of the koala specimens used by Tsangaras et al. (2014a) and Tsangaras et al. (2014b) and compared to our results were museum samples. By contrast, all the koalas examined in Ishida et al. (2015) were from modern samples.

For the 3′ integration sites, no sharing of integration sites between the museum samples in this study and museum samples in Tsangaras et al. (2014b) was detected. Two integration sites were found to be shared between the two youngest museum samples of the current study and Pci-SN265 (the only modern koala in Tsangaras et al. (2014b)). Moreover, one integration site was found shared between a modern koala (Pci-SN248) in Ishida et al. (2015) and Pci-SN265 of Tsangaras et al. (2014b). One integration site was also found shared between two museum koalas in the current study and modern koalas in Ishida et al. (2015) (Table S4).

Among 5′ integration sites, three were shared between the museum samples in this study and those used in Tsangaras et al. (2014b). Two integration sites were found to be shared between the museum samples of this study and Pci-SN265, and two integration sites were found shared between modern koalas (including Pci-SN248) in Ishida et al. (2015) and Pci-SN265. Additionally, four integration sites were found shared between relatively young museum koalas in this study and modern koalas in Ishida et al. (2015). A 5′ integration site (KoRV-5-shared_7) was shared by 9 koalas, including 4 museum koalas in this study, 4 museum koalas in Tsangaras et al. (2014b), and one modern koala (Pci-SN404) in Ishida et al. (2015) (Table S5). Statistical modeling of shared KoRV integration sites among 10 koalas showed an increased sharing of integration sites over time. The details are described in Article S1.