Transposon insertion profiling by sequencing (TIPseq) for mapping LINE-1 insertions in the human genome

Transposable elements make up a significant portion of the human genome. Accurately locating these mobile DNAs is vital to understand their role as a source of structural variation and somatic mutation. To this end, laboratories have developed strategies to selectively amplify or otherwise enrich transposable element insertion sites in genomic DNA. Here we describe a technique, Transposon Insertion Profiling by sequencing (TIPseq), to map Long INterspersed Element 1 (LINE-1, L1) retrotransposon insertions in the human genome. This method uses vectorette PCR to amplify species-specific L1 (L1PA1) insertion sites followed by paired-end Illumina sequencing. In addition to providing a step-by-step molecular biology protocol, we offer users a guide to our pipeline for data analysis, TIPseqHunter. Our recent studies in pancreatic and ovarian cancer demonstrate the ability of TIPseq to identify invariant (fixed), polymorphic (inherited variants), as well as somatically-acquired L1 insertions that distinguish cancer genomes from a patient’s constitutional make-up. TIPseq provides an approach for amplifying evolutionarily young, active transposable element insertion sites from genomic DNA. Our rationale and variations on this protocol may be useful to those mapping L1 and other mobile elements in complex genomes.

Characterizations of transposable element sequences remain incomplete in part because their highly repetitive nature poses technical challenges. Using these high copy number repeats as probes or primer sequences can create signals or products in hybridization-based assays and PCR amplifications that do not correspond to discrete genomic loci. Moreover, both the absence of many common insertion variants from the reference genome assembly * Correspondence: jef.boeke@nyumc.org; kburns@jhmi.edu 4 Institute for Systems Genetics, NYU Langone Health, New York, NY 10016, USA as well as the presence of hundreds of thousands of similar sequences together complicate sequencing read mappability. Detecting insertions that occur as low frequency alleles in a mixed sample presents an additional challenge, such as occurs with somatically-acquired insertions. Nevertheless, several recent studies describe strategies for mapping these elements and highlight LINE-1 continued activity in humans today. These methods include hybridization-based enrichment [24][25][26][27][28][29]; selective PCR amplification [6,17,[30][31][32][33][34][35][36][37][38][39]; and tailored analyses of whole genome sequencing reads [10,11,18,19,40,41].
Here we present a detailed protocol to amplify and sequence human LINE-1 retrotransposon insertion loci developed in the Burns and Boeke laboratories, Transposon Insertion Profiling by sequencing (TIPseq) [22,23,[42][43][44]. This method uses ligation-mediated, vectorette PCR [45] to selectively amplify regions of genomic DNA directly 3′ of L1Hs elements. This is followed by library preparation and Illumina deep sequencing (see Fig. 1a). TIPseq locates fixed, polymorphic, and somatic L1Hs insertions with base pair precision and determines orientation of the insertion (i.e., if it is on the plus (+) or minus (−) strand with respect to the reference genome). It detects, though does not distinguish between, both full length and 5′ truncated insertions as short as 150 bp. TIPseq is highly accurate in identifying somatic L1 insertions in tumor versus matched normal tissues, and allows sequencing coverage to be efficiently targeted to LINE-1 insertion sites so it is an economical way to process samples for this purpose. We have used TIPseq to demonstrate LINE-1 retrotransposition in pancreatic [22] and ovarian [23] cancers, and to show that somatically-acquired insertions are not common in glioblastomas [44]. Together with the machine learningbased computational pipeline developed in the Fenyӧ Lab for processing TIPseq data, TIPseqHunter [23], this protocol allows researchers to map LINE-1 insertion sites in human genomic DNA samples and compare insertion sites across samples.

Experimental design
Starting material and optimal reaction size High molecular weight genomic DNA is the starting material for TIPseq. This can be isolated from fresh or frozen tissues or cells. We typically use gDNA from phenol:chloroform extractions and ethanol precipitations, or from silica column preparations. This protocol uses reaction sizes producing consistent results in our hands with starting material of 10 μg genomic DNA (gDNA). We have successfully used a 3.3 μg gDNA input 'scaled-down' protocol with comparable results to the full-scale protocol. However, we caution that smaller reaction volumes will magnify effects of sample evaporation or slight inaccuracies in pipetting. It is important to maintain accurate reaction volumes at each step of the protocol. See Additional file 1: Table S1 for scaled-down reactions that start with as low as 3.3 μg of gDNA.

Restriction enzyme selection
TIPseq uses 6 different restriction enzyme digests run in parallel to maximize the portion of the genome that is cut to a PCR-amplifiable fragment in at least one of the reactions. The combination of enzymes was selected using a greedy algorithm to maximize genomic fragments 1-5 kb long. An L1Hs insertion occurring at any location in the genome is highly likely then to be represented by a fragment 1-3 kb in size in at least one of these parallel digests. This size balances informativeness and amplification efficiency; longer fragments include more sequence, but shorter fragments amplify more efficiently. For vectorette PCR to be successful, restriction enzymes should: 1) have a recognition cut site that occurs at the right genomic frequency (many 5-or 6-base pair cutters work well); 2) cut efficiently and independent of CpG methylation, 3) leave "sticky-end" overhangs for ligation of the vectorette adapters, and 4) be able to be heat inactivated. Most importantly, no restriction enzyme should cut in the retroelement insertion at any position 3′ of the forward primer sequence. This would prevent PCR amplicons from extending into unique gDNA downstream of the element.

Vectorette adapter design
Pairs of vectorette oligonucleotides are annealed together to form double stranded vectorette adapters (see Table 1). At one end of the vectorette, the two strands form compatible "sticky-ends" to the restriction enzyme digestion cut sites which allows for efficient adapter ligation (see Additional file 2: Table S2). The vectorette central sequence is partially mismatched such that the vectorette primer sequence is incorporated on the bottom strand, but its reverse complement is missing from the top strand. This forces first stranded synthesis to occur out of the transposable element to create the vectorette primer binding sequence. After this initial extension, exponential amplification can proceed in subsequent PCR cycles (see Fig. 1b).

Specific primer selection
The transposable element primer responsible for first strand synthesis is positioned in the 3' UTR of the LINE-1 sequence (see Fig. 2a). The primer placement takes advantage of 'diagnostics nucleotides' that define currently active LINE-1. The oligo ends with the ' ACA' trinucleotide located in the 3' UTR specific to the L1PA1 [also known as L1(Ta)] subset of Homo sapiens-specific LINE-1 (L1Hs). This strongly favors amplification of polymorphic and newly acquired somatic insertions and minimizes enrichment of older, "fixed present" elements.

Vectorette PCR conditions
Amplicons initiated within L1Hs insertions must traverse the LINE-1 polyA sequence and extend for a significant distance into downstream gDNA. We use a touchdown PCR program to ensure a balance between promoting primer specificity and achieving high-yields. This program progressively lowers the annealing temperature of each cycle from 72°C to 60°C (see Table 2). These cycling conditions, combined with the robust, proofreading DNA polymerase (ExTaq HS, Takara Bio; Shiga Japan), produces the complex mixture of optimally sized amplicons. ) data analysis. The first seven of these steps are shown adjacent to schematic representations in part b., to the right. b Vectorette adapter annealing is shown first. Mismatched sequences within the hybridized vectorette oligonucleotides are illustrated in red and blue, and create a duplex structure with imperfect base pairing. The sticky end overhang on one strand of the vectorette (here, a 5′ overhang on the bottom strand) is drawn in gray. This overhang in the annealed vectorette complements sticky ends left by genomic DNA digest, and the digest and vectorette ligations are shown in the subsequent two steps. The black box within the gDNA fragment illustrate a LINE-1 element of interest (i.e., a species-specific L1Hs). Most gDNA fragments will not have a transposable element of interest, and thus cannot be amplified efficiently by the vectorette PCR. In vectorette PCR, the L1Hs primer begins first strand synthesis (1) and extends this strand through the ligated vectorette sequence. The reverse primer complements this first-strand copy of the vectorette (2) and the two primers participate in exponential amplification (3) of these fragments in subsequent cycles. c Amplicons are sheared, and conventional Illumina sequencing library preparation steps complete the protocol. Paired-end sequencing reads are required to perform data analysis with TIPseqHunter. d A diagram of read pile-ups demonstrate how there is deep coverage of the 3′ end of L1Hs elements. For elements on the plus (+) strand with respect to the reference genome, the amplified sequences are downstream of the insertion site (i.e., covering genomic coordinates ascending from the transposon insertion). For minus (−) stranded insertions, sequences are recovered in the opposite direction

DNA shearing
We use a Covaris focused ultrasonicator (Covaris; Woburn, MA) with the manufacturer's recommended settings to shear the vectorette PCR amplicons to 300 bp prior to library preparation (see Additional file 3: Figure S2B). Shearing PCR amplicons may produce a broader size range than when shearing genomic DNA. If necessary, the treatment time may be modified on a per sample basis to adjust the final size distribution.

Library preparation and size-selection
Library construction may be performed using any kit that is compatible with Illumina next generation sequencing, including Illumina's TruSeq LT or PCR-free DNA sample prep kits (Illumina; San Diego, CA). We recommend using Kapa Library Preparation Kit for Illumina (Kapa Biosystems; Wilmington, MA) and to follow the manufacturer's instructions. If necessary, amplification may be performed during library construction, however, we advise using a PCR-free library preparation. Library adapters add approximately 120 bp of length to the sheared DNA. It may be necessary to perform a size selection during library preparation so that final library size is greater than 400 bp. This will prevent the generation of overlapping read pairs and reads containing adapter sequence. If necessary, we recommend performing dual-SPRI bead selection during library preparation or adding Pippin prep selection (Sage Science; Beverly, MA) after library pooling to remove all fragments smaller than 400 bp.

Illumina sequencing
Our data analysis pipeline, TIPseqHunter, requires 150 bp or shorter paired-end reads for optimal results. Longer reads may be trimmed to meet this requirement. We recommend a minimum of 15-25 million read pairs per sample. For example, for the Illumina HiSeq4000 this corresponds to pooling 12 samples per lane in highoutput mode. These guidelines should result in sufficient coverage and read depth for identifying L1 insertion loci.

Data analysis
TIPseq produces reads that contain LINE-1 sequence, adjacent genomic sequence, or both (junction reads) (see Fig. 2b). TIPseq data analysis reveals precise, base-pair resolution of L1Hs insertions and their orientation). We recommend using our custom bioinformatics program: TIPseqHunter [23]. We developed this program with a machine learning algorithm that uses known insertions as a training set for identifying new insertions. TIPseq-Hunter is available for download at: https://github.com/ fenyolab/TIPseqHunter (see Table 6). It is also available as a Docker image at: https://github.com/galantelab/tip-seq_hunter. This encapsulates all java dependencies, read aligners, genome indexes and biological annotation files needed by both steps of the pipeline. The genome indexes and annotation files in both TIPseqHunter and the Docker image use the human reference genome assembly GRCh37 (hg19). Instructions for use and download can be found in the README file at: https://github. com/galantelab/tipseq_hunter/blob/master/README.md. For sequencing runs of less than 20 million read pairs, 10-20 GB of RAM is suggested, and running time using 8 core processors on a Linux system is approximately 25 h. For runs in excess of 60 million reads, TIPseq-Hunter requires 40-50 GB of RAM, and running time is 1-1.5 h per 1 million reads. TranspoScope, a bioinformatics tool for browsing the evidence for transposable element

De novo insertion validation
TIPseqHunter accurately detects fixed, polymorphic, and de novo L1Hs insertions. Our previous studies have produced validation rates has high as 96% [23]. While users can therefore be confident in TIPseqHunter calls, we recommend validating at least subsets of predicted  Table 7). This will confirm the presence of the insertion and report the length and structure of the element. It is important to use the same high quality gDNA used in the TIPseq procedure to validation insertion candidates.
Normal control DNA should be tested in parallel when validating somatic insertions from tumor-normal studies (see Fig. 3a). L1-specific 3' PCR may be used to validate large insertions that are difficult to span in PCR and to identify possible 3′ transduction events (see Table 8).

Level of expertise required
The first part of the TIPseq protocol and final validations (steps 1-21,31) require basic molecular biology equipment and techniques (digestion, ligation, and PCR). The second part of the protocol (steps [22][23][24][25][26][27][28][29] involves the use of more advanced equipment and methods (DNA shearing, library preparation, and deep sequencing). It is possible to contract 'advanced' steps to sequencing core facilities depending on each user's level of expertise and access to the required equipment, and this is our recommendation for users

Applications of the method
TIPseq was initially adapted from a microarray based approach called Transposon insertion profiling by microarray or TIPchip [9,42], which was first developed for mapping Ty1 elements in Saccharomyces cerevisae [42] .
Although TIPseq is applicable to other transposable elements or species, this protocol is optimized to detect LINE-1 insertions in the human genome, and currently our TIPseqHunter program can only process human LINE-1 TIPseq data. TIPseq may be used for a variety of applications, including: population studies to identify common structural variants, tumor vs. normal comparisons to identify somatically-acquired insertions and trace cellular phylogenies, and in patients with specific phenotypes to evaluate for de novo retrotransposition events. Whole genome sequencing (WGS) can also be used for these purposes, and the main advantage of TIPseq is that insertion sites can be relatively deeply sequenced inexpensively. Targeting sequencing to retrotransposon insertion sites can result in a 400x cost saving for L1Hs mapping, and a 60x cost savings for Alu mapping.

Limitations of the method
Although TIPseq is a highly useful tool for detecting LINE-1 insertions, there are some limitations to the method that should be considered. First, TIPseq relies on restriction enzyme digestion of a large amount of high quality (high molecular weight) genomic DNA. For samples with limited amounts or reduced quality DNA, such as single-cell or fixed tissue, this protocol may need adjusted to work with similar efficiency. Secondly, while this method provides insertion location and orientation information, it does not differentiate between insertion 'types'. This includes classifying full length versus truncated insertions and elements with 5′ inversions or 3′ transductions (see Fig. 2a). While TIPseq will detect these insertions, further analysis, such as gel electrophoresis or Sanger sequencing, is required to confirm insert size and sequence variations. Finally, TIPseq does not distinguish between heterozygous and homozygous insertion alleles. An additional qualitative validation, such as PCR, is needed to confirm zygosity.

Anticipated results
The TIPseq procedure should yield more than 10 μg of purified PCR amplicons depending on vectorette PCR efficiency. The size distribution of these amplicons usually averages 1-3 kb (see Additional file 4: Figure S1A). This size distribution may vary depending on the quality of starting material. Sheared DNA should average around 300 bp (see Additional file 3: Figure S2B). Shearing of PCR amplicons produces a broader size range than when shearing gDNA. If necessary, shearing conditions may be adjusted to alter the final size distribution. The HiSeq4000 generates approximately 300 million read pairs per lane. Pooling up to 12 samples per lane will produce the recommended minimum of 15-25 million read pairs per sample. The final sequencing output consists of reads that align to the 3'UTR of LINE-1 and/or the adjacent genomic DNA. Read pairs will be either L1-genome, genome-genome, L1-junction, or junction-genome, or 'unpaired' genome (see Fig. 2b). On average, approximately 30 to 40% of TIPseq reads will align to LINE-1 sequence. Our validation rates for detecting new L1 insertions are as high as 96% [23]. TIPseq will identify full length and 5′ truncated L1's 150 bp and larger, including elements with 5′ inversions and 3′ transductions. However, additional PCR and Sanger sequencing must be performed to confirm these events (see Table 8).

Conclusions
This protocol describes in detail our approach to transposon insertion profiling by next-generation sequencing (TIPseq). The assay as described targets signature sequences in the 3'UTR of evolutionarily young L1PA1 elements for insertion site amplification. A subset of these elements is active in the modern human genome. Their ongoing activity makes them valuable to map for characterizing heritable genetic polymorphisms, de novo insertions, and somatic retrotransposition activity. While LINE-1 insertion sites can be detected in whole genome sequencing data, selectively amplifying these sites can allow investigators to target their sequencing to insertion locations. This enables LINE-1-directed studies to more efficiently and affordably use sequencing and computational resources. We have demonstrated that variations of this protocol are effective at selectively amplifying other transposable element in humans [i.e., Alu insertions (See Additional file 5: Table  S3), and endogenous retroviruses (ERV-K)], and we expect that similar approaches can be taken to map active mobile genetic elements, other high-copy recurrent sequences, or transgene insertions.

Reagent setup Genomic DNA
TIPseq requires starting with high molecular weight genomic DNA. We recommend isolating fresh gDNA when possible. Poor quality genomic DNA will reduce TIPseq's efficiency. Always avoid vortexing, rough pipetting, and excessive freeze-thaw cycles to ensure gDNA integrity is maintained throughout the protocol.

Oligonucleotide stocks
Vectorette adapter oligonucleotides should be resuspended with TE buffer to stock concentrations of 100 μM. PCR primers should be resuspended with molecular grade water to stock concentrations of 100 μM. Stocks should be stored at − 20°C, thawed and mixed well before use.

Master mix preparations
All master mixes should be prepared on ice immediately before use. We recommend including a 2-3 sample excess when preparing each master mix. See Tables 3, 4, 5 for master mix formulas.

Equipment setup Thermal cycler
We recommend performing the restriction enzyme digestions, inactivation steps, and PCR in a pre-heated thermal cycler with heated lid.
(CAUTION Ethidium bromide is toxic and is a potential mutagen and carcinogen. Use proper protective wear.) The gel should be run at a constant 100 V for 45 min or until separation of the ladder is clearly visible.

Covaris shearing system
The Covaris LE220 shearing system is setup according to the manufacturer's instructions.

Procedure
Steps 1-5: Vectorette adapter annealing (Timing: 2 h) Steps 6-9: Genomic DNA digestion (Timing: 1 h setup and overnight incubation) 6. Dilute 10 μg genomic DNA in 123.5 μL of molecular grade water and aliquot diluted gDNA to each of six 0.2 mL PCR tubes 7. Prepare digestion master mix on ice for the appropriate number of samples plus excess (See Table 3). Mix by gently pipetting the entire volume 5 times and quickly spin to collect. 8. Add 6 μL of digestion master mixes in parallel to each gDNA aliquot. Mix by gently flicking and spinning. 9. Incubate overnight at the appropriate activation temperature in a thermal cycler with heated lid.
Steps 10-14: Vectorette adapter ligation (Timing: 3 h setup and overnight incubation) 10. Inactivate the restriction enzyme digests for 20 min at 80°C in thermal cycler with heated lid. Cool to room temperature. 11. Add 2 μL of the appropriate 1 μM annealed vectorettes adapters to each digest and mix by gently flicking and spinning. Critical: Be sure to add each annealed vectorette to its corresponding enzyme digest. 12. Use a thermal cycler with heated lid to incubate at 65°C for 5 min and then slowly cool to room temperature (0.5°C/min). Move samples to 4°C for at least 1 h. 13. Prepare ligation master mix on ice for the appropriate number of samples plus excess (See Table 4). Mix by gently pipetting the entire volume 5 times and quickly spin to collect.    Table 2). The program can be left to run overnight. Critical: BAM file has to be generated by bowtie2 alignment with "XM" tag Running TIPseqHunter: (1) for quality control, alignment, feature selection, modeling, prediction: ./TIPseqHunterPipelineJar.sh fastq_path output_path fastq_r1 key_r1 key_r2 num_rp Critical: Detailed information is provided in the TIPseqHunterPipelineJar.sh file. Some parameters need to be pre-set. Parameters: fastq_path: path of the fastq files (Note: this is the only path and file name is not included) output_folder: path of the output files (Note: this is the only path and file name is not included) fastq_r1: read 1 file name of paired fastq files key_r1: key word to recognize read-1 fastq file (such as "_1" is the key word for CAGATC_1.fastq fastq file) Critical: key has to be unique in the file name key_r2: key word to recognize read-2 fastq file and replaceable with the read-1 key word to match to read-1 file (such as "_2" is the key word for CAGATC_2.fastq fastq file) Critical: key has to be unique in the file name num_rp: the total number of the read pairs in the paired fastq files (Note: it is the total number of read-pairs, i.e. either the total number of read1 or read2 but not together.) (This number is for normalization purpose) (2) for somatic insertions: TIPseqHunterPipelineJarSomatic.sh repred_path control_path repred_file control_file Critical: Detailed information is provided in the TIPseqHunterPipelineJarSomatic.sh file. Some parameters need to be pre-set. Parameters: repred_path: path of "model" folder under output folder control_path: path "TRLocator" folder under output folder repred_file: file with suffix ".repred" and generated from P11 in repred_path (Note: file name should be ending with ".repred".) ( Troubleshooting: If PCR yield is too low, restart procedure with freshly annealed vectorette adapters, isolate fresh gDNA, or increase the initial amount of gDNA. 21. Run 2 μg of purified DNA on 1.5% agarose gel. Critical: Vectorette PCR amplicons should appear as a smear on the gel averaging around 1-3 kb. (see Additional file 4: Figure S1A). Troubleshooting: The presence of a very high molecular weight smear could indicate primervectorette concatemer amplification. Digest 2 μg of purified vectorette PCR amplicons with BstYI and run on a 1.5% agarose gel. BstYI cuts within the vectorette primer. An intense band around 50 bp indicates the presence of vectorette-primer concatemers in the PCR product (see Additional file 4: Figure S1B).  Figure S2B).
Steps 26-28: Library preparation and quality control (Timing: 1 d) 26. Use 200 ng of sheared DNA to prepare libraries using KAPA Library Preparation Kit for Illumina Critical: Use the same high quality gDNA that served as starting material for TIPseq 4. Run the PCR in a thermal cycler with heated lid using a 10-minute extension for 30 cycles. Critical: It is necessary to use an extension time long enough to amplify a full length, 6kb L1 insertion. 5. Run the PCR product on a 1% agarose gel and excise the band containing the filled allele (see Fig. 3a). Troubleshooting: If no filled band occurs, we recommend trying a 3' L1 specific PCR. (see Table 8). 6. Purify the excised DNA using Zymoclean Gel DNA Recovery Kit following the manufacturer's instructions. 7. Sanger sequence the purified DNA using both the forward and reverse PCR primer Critical: Use the same high quality gDNA that served as starting material for TIPseq 4. Run the PCR in a thermal cycler with heated lid using a 60°C annealing temperature and at least a 30-second extension for 30 cycles. Critical: It is important to use a slightly higher annealing temperature and shorter extension time to reduce the amount of off-target L1 binding and amplification. 5. Run the PCR product on a 1% agarose gel and excise the band from the successful reaction. Critical: Only one of the two PCR reactions should produce a band. A plus stranded L1 insertion will produce a band in the reverse primer reaction, and a minus stranded L1 will produce a band in the forward primer reaction. The size of the band should equal the distance from the genomic primer to the L1 insertion site plus 150bp of L1 and polyA sequence. A band larger than expected could indicate a 3' transduction event has occurred (See Fig. 3b) 6. Purify the excised DNA using Zymoclean Gel DNA Recovery Kit following the manufacturer's instructions. 7. Sanger sequence the purified DNA using the L1 primer and either the forward or the reverse genomic primer, depending on which reaction was successful. Critical: It may be necessary to use internal primers to sequence through the product completely.
according to the manufacturer's instructions without performing dual-SPRI size selection. Critical: Avoid performing library amplification. We recommend avoiding size selection, but dual-SPRI bead selection may be performed. 30. Analyze data using TIPseqHunter (see Table 6).
Troubleshooting: If the data contain a large amount of overlapping read pairs, use Pippin prep selection after pooling (step 28) to remove fragments under 400 bp. 31. Perform PCR validation and Sanger sequencing (see Tables 7 and 8)

Timing
Steps

Troubleshooting
See Table 9 for troubleshooting information.  Figure S1B). Repeat procedure with fresh reagents in an amplification-free area.   Table 8)