Genomic analysis of P elements in natural populations of Drosophila melanogaster

Estimation of P element content in isofemale strain and pool-seq genomic samples

Whole genome shotgun sequences for D. melanogaster from three geographic regions—N. America (Athens, Georgia, USA), Europe (Montpellier, France), and Africa (Accra, Ghana), described in Verspoor & Haddrill (2011)—were downloaded from the European Nucleotide Archive (ERP009059) (Bergman & Haddrill, 2015). These genomic data were collected from adult females using a uniform library preparation and sequencing strategy (thus mitigating many possible technical artifacts) and include data from both individual isofemale strains and pools of single flies derived from isofemale strains (see (Bergman & Haddrill, 2015) for details). The median insert size of libraries analyzed in this study ranged from 453–490 bp with a mean (±S.D.) of 475  ± 8 bp. The mean depth of coverage of the libraries analyzed in this study ranged from 48×–56× with a mean (±S.D.) of 53 ×  ± 1.8. In total, 50 isofemale strain genomes were analyzed: N. America (n = 15), Europe (n = 20), and Africa (n = 15). Two pool-seq samples were analyzed for each population (N. America (15 and 30 strains), Europe (20 and 39 strains), and Africa (15 and 32 strains)). The smaller pools for each population include one individual from the same isofemale strains that were sequenced as strain-specific samples from that population; the larger pools contain one individual from the same strains as the smaller pools, plus individuals from additional isofemale strains from the same population that were not sequenced as strain-specific samples.

Estimates of the number of P elements per sample were generated using two different non-reference TE detection methods (TEMP (Zhuang et al., 2014) and RetroSeq (Keane, Wong & Adams, 2013)), as well as by computing coverage of the P element normalized to the D. melanogaster reference genome. Non-reference TE detection methods like TEMP and RetroSeq can identify individual insertion sites, but cannot differentiate full-length from truncated copies and have detection biases that are not fully characterized. Therefore, normalized coverage across diagnostic sections of the P element consensus sequence was used to provide a complementary estimate of copy number and to estimate the abundance of full-length and truncated P element copies.

TEMP (revision d2500b904e2020d6a1075347b398525ede5feae1) and RetroSeq (revision 700d4f76a3b996686652866f2b81fefc6f0241e0) were run using the McClintock pipeline (revision 3ef173049360d99aaf7d13233f9d663691c73935; Nelson, Linheiro & Bergman, 2017). McClintock was run across the major chromosome arms (chr2L, chr2R, chr3L, chr3R, chr4, chrY, and chrX) of the UCSC dm6 version of the Release 6 reference genome (Hoskins et al., 2015) using the following options: -C -m “retroseq temp” -i -p 12 -b. We note that TEMP and RetroSeq both base their predictions on a common BAM file generated by the McClintock pipeline, and thus read mapping is controlled for in our analysis. Reference TE annotations needed for TEMP were generated automatically by McClintock using RepeatMasker (version open-4.0.6). The D. melanogaster TE library used by McClintock to predict reference and non-reference TE insertions is a slightly modified version of the Berkeley Drosophila Genome Project TE data set v9.4.1 (https://github.com/cbergman/transposons/blob/master/misc/D_mel_transposon_sequence_set.fa; described in Sackton et al., 2009).

In addition to providing “raw” output for each component method in the standardized zero-based BED6 format, McClintock generates “filtered” output tailored for each method (Nelson, Linheiro & Bergman, 2017). McClintock filters TEMP output to: (i) eliminate predictions where the start or end coordinates had negative values; (ii) retain predictions where there is sequence evidence supporting both ends of an insertion; and (iii) retain predictions that have a ratio of reads supporting the insertion to non-supporting reads of >1/10. Likewise, McClintock filters RetroSeq output to: (i) eliminate predictions where two different TE families shared the same coordinates; and (ii) retain predictions assigned a call status of greater than or equal to six as defined by Keane, Wong & Adams (2013).

Isofemale strain genomes from these populations have been reported to contain residual heterozygosity (Lack et al., 2016), presumably due to incomplete inbreeding in these strains prior to genome sequencing. To control for potential incomplete inbreeding in strain-specific sequences, we also generated estimates of P element abundance based on the sample-frequency weighted-sum of TEMP predictions. Specifically, we weighted each P element prediction in the raw or filtered TEMP BED files by the frequency of that prediction in the original TEMP insertion.refined.bp.summary output file for that sample. Frequency weighted estimates could not be generated for RetroSeq since this method does not provide an estimate of sample frequency for non-reference TE predictions.

Normalized P element coverage was computed with SAMtools (v.1.4.1) (Li et al., 2009) using the BAM file generated by McClintock that was used for TEMP and RetroSeq predictions above. The RepeatMasker annotation of TEs in the dm6 reference genome generated by McClintock was used to generate a complementary BED file of unique regions in the major chromosome arms (chr2L, chr2R, chr3L, chr3R, chr4, chrY, and chrX). SAMtools depth was then used to calculate the mean coverage of the unique regions of the dm6 genome and across the canonical P element sequence. Mean coverage per site was then averaged across the length of P element exon 0 (positions 153–442) and exon 2 (positions 1222–1947), and the resulting value was divided by mean coverage of the unique regions of the dm6 genome to generate normalized coverage of each exon. Normalized coverage of exon 0 was used to estimate overall P element copy number. Normalized coverage of exon 2 was used to estimate the number of full-length P elements (FP). Normalized coverage of exon 0 minus exon 2 was used to estimate the number of truncated P elements that are predominantly of the KP type (Black et al., 1987). Coverage profiles across the entire P element sequence were plotted as the average [±95% C.I.] across all isofemale lines and as the loess smoothed (span  = 0.01) values for pool-seq samples.

Gonadal dysgenesis phenotypes

We re-analyzed GD assay data from Ignatenko et al. (2015) for 43 of the 50 isofemale strains analyzed for genomic P element content above. Definitions of A and A* crosses in Ignatenko et al. (2015) are inverted relative to those proposed by Engels & Preston (1980), and were standardized prior to re-analysis here. Cross A measures the activity of tester strain males mated to M-strain Canton-S females; cross A* measures the susceptibility of tester females mated to a P-strain Harwich males. P, P′, Q, and M-strains were defined according to Kidwell, Frydryk & Novy (1983) and Anxolabéhère & Quesneville (1998).

Graphical and statistical analyses were performed in the R programming environment (version 3.4.0).

P element content varies among isofemale strains from N. America, Europe and Africa

We initially generated estimates of P element abundance in the genomes of 50 isofemale strains from three worldwide populations using two non-reference TE detection methods: TEMP (Zhuang et al., 2014) and RetroSeq (Keane, Wong & Adams, 2013) (Files S1, S2). TEMP and RetroSeq were chosen for P element detection because they showed the highest performance to accurately detect non-reference TE insertions in a previous simulation study (Nelson, Linheiro & Bergman, 2017). We also generated estimates of P element abundance using a complementary approach based on coverage of reads mapped to the canonical P element sequence normalized to coverage of the D. melanogaster reference genome (File S1). To account for potential false positive predictions, we also investigated the effects of the default filtering of TEMP and RetroSeq output performed by McClintock (Nelson, Linheiro & Bergman, 2017), a meta-pipeline that runs and parses multiple TE insertion detection methods. Additionally, to deal with potential residual heterozygosity in strain-specific genome sequences, we investigated weighting TEMP predictions by their sample frequency. We note that neither TEMP nor RetroSeq can differentiate full-length from truncated insertions in their output, and conversely normalized coverage cannot resolve the location of individual insertions. Moreover, we note that heterochromatic contigs were omitted from TEMP and RetroSeq analysis, however heterochromatic insertions can contribute to estimates of P element content based on coverage.

As shown recently in related work by Srivastav & Kelleher (2017), we find that overall numbers of P elements predicted in strain-specific genome sequences were well correlated across all bioinformatics approaches investigated (r ≥ 0.65; P < 2.5e − 7) (Fig. 1). The lowest correlation between distinct methods was for raw TEMP predictions and normalized coverage (r = 0.71) (Figs. 1G and 1AQ), and the highest correlation was for the filtered TEMP and filtered RetroSeq datasets (r = 0.96) (Figs. 1M and 1AK). McClintock filtering reduced the numbers of predicted P elements per strain substantially for TEMP but less so for RetroSeq (Fig. 1, Table 1) and generally increased the correlation of TEMP and RetroSeq predictions with other methods (Fig. 1). These results suggest that all bioinformatic approaches investigated capture similar aspects of variation in P element abundance across strains, and that filtering by McClintock improves the consistency of TE predictions made by TEMP and RetroSeq on isofemale strains.

Despite providing correlated estimates of abundance across strains, the absolute numbers of estimated P elements in these strain-specific genome data differed substantially across methods (Table 1). Raw output from TEMP and RetroSeq indicate that flies from these isofemale lines contain ∼80–250 P element insertions, whereas McClintock filtered predictions for TEMP and RetroSeq suggest they contain ∼70–120 P element insertions. In contrast, estimates of abundance based on frequency-weighted TEMP predictions and normalized coverage indicates that flies from these samples only contain ∼20–45 P elements per haploid genome, with weighted TEMP predictions giving slightly higher estimates of P element abundance than coverage. Previous results based on Southern blotting suggest 30–50 P elements per haploid genome (Bingham, Kidwell & Rubin, 1982; Ronsseray, Lehmann & Anxolabéhère, 1989; Bonnivard & Higuet, 1999; Itoh & Boussy, 2002). We interpret the discrepancy between unweighted non-reference TE detection methods and either weighted TEMP or normalized coverage, together with the similarity between the latter approaches and Southern hybridization data, as being consistent with overestimation by unweighted non-reference TE detection methods due to residual heterozygosity stemming from incomplete inbreeding in these strains (Lack et al., 2016). If multiple haplotypes persist at many loci in these strains due to partial inbreeding, this could elevate the number of insertion sites per strain predicted by unweighted non-reference TE detection methods (relative to the number expected from a haploid genome), but not increase estimated copy number based on weighted non-reference TE detection methods or normalized coverage (because polymorphic insertion sites would only contribute to abundance for these methods in a manner that is proportional to their frequency). The relatively large decrease in the number of predictions observed in raw versus filtered output for TEMP but not RetroSeq is also consistent with this interpretation since TEMP (but not RetroSeq) is capable of identifying polymorphic insertions in a genomic sample. Since estimates of P element abundance based on weighted TEMP predictions and coverage are concordant with independent estimates based on Southern hybridization, we conclude that these methods provide more accurate estimates of P element abundance per strain and therefore focused further analysis on these two approaches.

Despite differences in how they capture aspects of TE abundance, estimates of P element abundance based on weighted TEMP predictions or coverage show significant differences in the number of P elements per strain across populations (One-way ANOVA; F > 35.27; 2 d.f., P < 1.3e − 9), with strains from the N. American population showing a clear reduction in P element abundance relative to the European and African populations (Table 1; Figs. 2A–2C). Similar results are observed for abundance estimates based on unweighted TEMP and RetroSeq predictions (One-way ANOVA; F > 9.26; 2 d.f., P < 5e−4) (Fig. S1). Moreover, by contrasting normalized coverage in exon 0 versus exon 2 of the P element to estimate numbers of full-length (FP) and truncated (KP) copies, we found clear differences among populations in the number of full-length (One-way ANOVA; F = 32.4; 2 d.f., P = 4.33e−9) and truncated (One-way ANOVA; F = 102.1; 2 d.f., P < 2e−16) copies, with slightly more full-length and many fewer truncated copies in the N. American population relative to the European and African populations (Figs. 2D–2E). Together, this variation in P element structure leads to a strong difference in the KP∕FP ratio across populations (One-way ANOVA; F = 20.0; 2 d.f., P < 9.52e−7) (Fig. 2F).

Differences across populations in overall P element abundance and the relative abundance of full-length and truncated copies can be visualized by plotting normalized coverage averaged across isofemale strains over the length of the P element canonical sequence (Fig. 3). Relative to European and African strains, N. American strains have lower coverage at the termini of the P element including exon 0 (positions 153–442, present in truncated and full-length copies), but slightly higher coverage in exon 2 (positions 1222–1947, present only in full-length copies). European and African strains also have a steep changes in coverage precisely at the breakpoints expected for the KP element (positions 807 and 2561), which are not observed in N. American strains. Together, these results demonstrate clear differences in P element abundance and structure among isofemale strains from N. American versus European and African populations. Additionally, these coverage profiles reveal positional variation in depth across exon 0 in all populations, which may explain the slightly lower abundance estimates obtained from coverage relative to weighted TEMP predictions.

Lack of association between P element content and gonadal dysgenesis phenotypes

To address whether the variation in genomic P element content that is observed across populations affects geographic variation in hybrid dygenesis, we compared P element abundance from 43 isofemale strains analyzed above with previously-published GD phenotypes for the same strains (Ignatenko et al., 2015). As reported in Ignatenko et al. (2015), isofemale strains from these populations were mainly P, M and Q (Table 2, Fig. S2). Based on our genomic analysis (File S1), all strains in these populations that are defined phenotypically as M are actually M’. For N. American and African populations, the degree of P element activity tends to vary more across strains within populations relative to susceptibility (Table 2, Figs. 2G–2H). However, we found no evidence for systematic differences across populations in the degree of activity (One-way ANOVA; F = 0.06, 2 d.f., P = 0.94) or susceptibility (One-way ANOVA; F = 1.66, 2 d.f., P = 0.2). Taken together with the genomic analysis above, these results suggest that population-level differences in the abundance or structure of P elements do not necessarily translate into differences in the frequency of GD phenotypes across populations.

Integrating our genomic predictions with these published GD data at the level of individual strains, we directly tested whether the number of P elements per strain is associated with either activity or susceptibility phenotypes in GD assays. We found that neither the degree of activity nor the degree of susceptibility was significantly correlated with P element abundance estimated from weighted TEMP predictions or coverage (P > 0.11; Fig. 4). Similar results were obtained using the unweighted output of TEMP and RetroSeq as well (P > 0.11; Fig. S3). Generalized linear models with Gaussian errors incorporating population as a factor also failed to identify significant associations between any measure of P element abundance and GD phenotypes when corrected for multiple testing (P > 0.038). These results suggest that there is no simple relationship between the total abundance or relative abundance of full-length versus truncated P elements and GD phenotypes at the level of individual strains, confirming conclusions inferred from the population level above.

Population difference in P element content can be observed in pool-seq samples

Pooled-strain sequencing (pool-seq) is a cost-effective strategy to sample genomic variation across large numbers of strains and populations (Schlotterer et al., 2014). However, little empirical evidence exists reporting how biological conclusions derived from pool-seq data relate to those from strain-specific sequencing for the same samples or populations. The populations analyzed here provide a unique opportunity to address such questions since pool-seq samples are available for the same set of isofemale strains with strain-specific genomic data. Two pool-seq samples are available for each population that differ in the number of individuals (one per isofemale strain) used: N. America (n = 15 and n = 30), Europe (n = 20 and n = 39), and Africa (n = 15 and n = 32). The smaller pools for each population each include one individual from the same isofemale strains that were sequenced as strain-specific samples from that population; the larger pools contain one individual from the same strains as the smaller pools, plus individuals from additional isofemale strains from the same population that were not sequenced as strain-specific samples. Thus, the smaller pool-seq samples are a nested subset of the larger pool-seq samples, and pool-seq samples from the same population are not fully independent from one another.

To address whether the differences among populations we observed in the number of P elements observed in isofemale strain data can also be detected using different bioinformatics approaches in pool-seq data, we estimated the number of P elements in pool-seq samples from the same populations using weighted TEMP predictions and coverage (Table 3). Because depth of coverage per strain in pool-seq samples is inversely proportional to the number of strains in the pool, estimates of P element content based on weighted TEMP predictions or normalized coverage offer the most straightforward approach to compare strain-specific and pool-seq results. As expected if pool-seq samples accurately reflect the average based on strain-specific sequencing, estimates of P element abundance in pool-seq samples from raw weighted TEMP predictions and coverage (Table 3) are similar to those derived from strain-specific sequencing of the same population (Table 1). For all pool-seq samples, however, estimates of the number of P elements per strain based on raw weighted TEMP predictions or coverage are slightly higher than the average weighted TEMP predictions based on isofemale strain data. In contrast, estimates of P element abundance based on filtered weighted TEMP predictions in pool-seq data are much lower than other methods applied to pool-seq data or for filtered weighted TEMP predictions on strain-specific genomic data, especially for the larger pools. Underestimation of P element abundance from filtered weighted TEMP predictions presumably arises because many real P element insertions in these pool-seq samples exist at a frequency lower than the cutoff of 1/10 required by McClintock filtering. This interpretation is consistent with the fact that the vast majority of non-redundant P element insertion sites in strain-specific genomic data are found only in a single strain (>75%), and these singleton insertions would therefore be at a frequency lower than 1/10 in any of the strain pools.

Importantly, estimates of P element abundance based on raw weighted TEMP predictions or coverage in N. America pool-seq samples show a clear reduction in the estimated number of P elements relative to European and African pool-seq samples (Table 3), consistent with the statistically significant differences among populations observed using samples of strain-specific genomic data above. Furthermore pool-seq samples qualitatively and quantitatively (with the exception of the pool-20 sample in Europe) recover population differences in depth of coverage profiles across the canonical P element sequence (Fig. 3). These conclusions are true both for the smaller pools that precisely match the strain-specific samples, as well as the larger pools that include genomes not sampled in the strain-specific data. These results suggest that frequency-weighted raw TEMP predictions and coverage are viable methods for detecting differences in P element abundance and structure across D. melanogaster populations using pool-seq data.