Retrotransposons Are the Major Contributors to the Expansion of the Drosophila ananassae Muller F Element

Abstract

The discordance between genome size and the complexity of eukaryotes can partly be attributed to differences in repeat density. The Muller F element (∼5.2 Mb) is the smallest chromosome in Drosophila melanogaster, but it is substantially larger (>18.7 Mb) in D. ananassae To identify the major contributors to the expansion of the F element and to assess their impact, we improved the genome sequence and annotated the genes in a 1.4-Mb region of the D. ananassae F element, and a 1.7-Mb region from the D element for comparison. We find that transposons (particularly LTR and LINE retrotransposons) are major contributors to this expansion (78.6%), while Wolbachia sequences integrated into the D. ananassae genome are minor contributors (0.02%). Both D. melanogaster and D. ananassae F-element genes exhibit distinct characteristics compared to D-element genes (e.g., larger coding spans, larger introns, more coding exons, and lower codon bias), but these differences are exaggerated in D. ananassae Compared to D. melanogaster, the codon bias observed in D. ananassae F-element genes can primarily be attributed to mutational biases instead of selection. The 5' ends of F-element genes in both species are enriched in dimethylation of lysine 4 on histone 3 (H3K4me2), while the coding spans are enriched in H3K9me2. Despite differences in repeat density and gene characteristics, D. ananassae F-element genes show a similar range of expression levels compared to genes in euchromatic domains. This study improves our understanding of how transposons can affect genome size and how genes can function within highly repetitive domains.

Keywords: Drosophila; Wolbachia; genome size; heterochromatin; retrotransposons.

Figure 1

Results of the manual sequence improvement of the D. ananassae D and F elements. (A) Dot plot comparisons of the scaffolds in the original CAF1 assembly (y-axis) vs. the scaffolds in the improved assembly (x-axis). Dots within each dot plot denote regions of similarity between the CAF1 assembly and the improved assembly. The diagonal lines in the dot plots for the D-element scaffold improved_13337 (left) and the F-element scaffold improved_13010 (middle) show that the overall CAF1 assemblies for these regions are consistent with the corresponding assemblies following manual sequence improvement. However, the high density of dots in the middle of the dot plot for improved_13010 corresponds to a collapsed repeat within the CAF1 assembly (red ←). Manual sequence improvement also identified a major misassembly in the second improved region (improved_13034_2) within the F-element scaffold_13034 (right), where part of the scaffold was inverted compared to the final assembly (red box). The misassembled region is part of the fosmid 1773K10 (bottom inset). (B) The Consed Assembly View for the improved fosmid project 1773K10 shows that the misassembled region contains multiple tandem and inverted repeats. (Top) The gray bar within the Assembly View corresponds to the improved fosmid assembly, and the pink Δ’s denote the ends of the fosmid. The purple and green boxes underneath the gray box correspond to tags (e.g., repeats and comments), and the dark green line corresponds to the read depth. The orange and black boxes above the gray bar correspond to tandem and inverted repeats, respectively. These orange and black boxes indicate that the improved assembly contains multiple tandem and inverted repeats that are located adjacent to each other. (Bottom) The improved assembly for this fosmid was supported by consistent forward-reverse mate pairs (blue Δ’s) and by multiple restriction digests (Figure S2A in File S7).

Figure 2

Expansion of the D. ananassae F element can primarily be attributed to the high density of LTR and LINE retrotransposons. (A) Total repeat density estimates from de novo repeat finders (Red, WindowMasker, and Tallymer) show that the D. ananassae F element has higher repeat density than the D. melanogaster F element (56.4–74.5% vs. 14.4–29.5%). Both F elements show higher repeat density than the euchromatic reference regions from the base of the D elements in D. ananassae (11.6–19.8%) and in D. melanogaster (5.4–15.2%). The repeat densities of the improved D. ananassae F-element scaffolds [D. ana: F (improved)] are similar to the repeat densities for all D. ananassae F-element scaffolds [D. ana: F (all)]. (B) Results from the tantan analysis show that the five analysis regions from D. melanogaster and D. ananassae have similar simple repeat density (6.0–7.5%). (C) TRF analysis shows that D. ananassae has higher tandem repeats density than D. melanogaster on both the F (5.6 vs. 2.8%) and the D elements (2.6 vs. 1.1%). (D) RepeatMasker analysis using the Drosophila RepBase library shows that the F element has higher transposon density than the D element both in D. melanogaster (28.0 vs. 7.7%) and in D. ananassae (78.6 vs. 14.4%). There is a substantial increase in the density of LTR and LINE retrotransposons on the D. ananassae F element compared to the D. melanogaster F element (42.1 vs. 5.5% and 21.8 vs. 3.8%, respectively). The D. ananassae euchromatic reference region also shows higher transposon density than D. melanogaster, but most of the difference can be attributed to the density of DNA transposons (4.3 vs. 0.6%). The region of overlap between two repeat fragments is classified as “Overlapping” if the two repeats belong to different repeat classes.

Figure 3

The high density of “Wolbachia” sequences in the D. ananassae F element can be attributed to Drosophila transposons in the wAna assembly. (A) RepeatMasker analysis shows that 19.8% of the D. ananassae F element matches the genome assembly for wAna. By contrast, 0.02% of the D. ananassae F element matches the genome assemblies for wRi and wMel. Similarly, the D. ananassae D element and the D. melanogaster F and D elements show a substantially higher density of regions that exhibit sequence similarity to the wAna assembly (0.68–6.51%) than the wRi and wMel assemblies (0.00–0.18%). (B) Distribution of regions with matches to the wAna, wRi, and wMel assemblies in the manually improved region of the D. ananassae F-element scaffold improved_13010. The matches to the wAna assembly are distributed throughout the improved scaffold (blue boxes) but there are no matches to either the wRi or the wMel assemblies. (C) The portions of the 3.6-kb wAna scaffold AAGB01000087 (x-axis) with large numbers of alignments to D. ananassae scaffolds show similarity to Drosophila transposons. The portions of the wAna scaffolds that show similarity to D. ananassae scaffolds were extracted from the RepeatMasker output and collated into an alignment coverage track relative to the wAna assembly (brown graph). Whole-genome Chain and Net alignments show that only the last 216 bp of this wAna scaffold has sequence similarity to the wRi and wMel assemblies. RepeatMasker analysis using the Drosophila RepBase library shows that the first 3.4 kb of this wAna scaffold has sequence similarity to the internal and long terminal repeat portions of the BEL-18 LTR retrotransposon (BEL-18_DAn-I and BEL-18_DAn-LTR), as well the internal portion of the Gypsy-1 LTR retrotransposon from D. sechellia (Gypsy-1_DSe-I). Most of the alignments can be attributed to BEL-18_DAn-LTR, with a maximum of 1835 alignments between the wAna scaffold AAGB01000087 and the D. ananassae scaffolds.

Figure 4

F-element genes show distinct gene characteristics compared to D-element genes. Each violin plot is comprised of a box plot and a kernel density plot. The ● in each violin plot denotes the median and the darker region demarcates the IQR, which spans from the first (Q1) to the third (Q3) quartiles. The whiskers extending from the darker region spans from Q1 = −1.5 × IQR to Q3 = +1.5 × IQR; data points beyond the whiskers are classified as outliers. (A) D. ananassae F-element genes have larger coding spans (start codon to stop codon, including introns) than D. melanogaster F-element genes. (B) The D. ananassae F-element genes have larger coding spans because they have larger total intron sizes than D. melanogaster F-element genes. (C) F-element genes have larger total coding exon (CDS) sizes than D-element genes in both D. ananassae and D. melanogaster. (D) F-element genes have more CDS than D-element genes. (E) F-element genes have smaller median CDS size than D-element genes. (F) The median intron size for D. ananassae and D. melanogaster F-element genes shows a bimodal distribution; this distribution pattern indicates that the expansion of the coding spans of D. ananassae F-element genes compared to D. melanogaster F-element genes can be attributed to the substantial expansion of a subset of introns.

Figure 5

Codon bias in D. ananassae F-element genes can primarily be attributed to mutational biases instead of selection. (A) Violin plots of the Nc show that D. ananassae F-element genes exhibit stronger deviations from equal usage of synonymous codons (lower Nc) than D. melanogaster F-element genes. (B) Violin plots of the CAI show that D. ananassae F-element genes exhibit less optimal codon usage (lower CAI) than D. melanogaster F-element genes. F-element genes in both D. ananassae and D. melanogaster show less optimal codon usage (lower CAI) than D-element genes. (C) Scatterplot of Nc vs. CAI suggests that the codon bias in most D. ananassae and D. melanogaster F-element genes can be attributed to mutational bias instead of selection, as indicated by the placement of most of the genes in the portion of the LOESS regression line (red line) with a positive slope. By contrast, the codon bias in most D. ananassae and D. melanogaster D-element genes can be attributed to selection, as denoted by the LOESS regression line with a negative slope. The dotted line in each Nc vs. CAI scatterplot corresponds to the CAI value for a gene with equal codon usage relative to the species-specific reference gene sets constructed by the program scnRCA (0.200 for D. ananassae and 0.213 for D. melanogaster). Hence this species-specific threshold corresponds to the CAI value when the strengths of mutational bias and selection on codon bias are the same. A smaller percentage of F-element genes in D. ananassae (6/64; 9.4%) have CAI values above this species-specific CAI threshold compared to D. melanogaster (18/79; 22.8%).

Figure 6

Metagene analysis shows that the coding spans of F-element genes have lower median 9-bp T_m than the genes at the base of the D element. The T_m profiles were determined using a 9-bp sliding window with a step size of 1 bp. The metagene consists of the 2-kb region upstream and downstream of the coding span, with the length of the coding spans normalized to 3 kb. While the codons in D. ananassae F-element genes have lower GC content, D. ananassae F-element genes exhibit a T_m profile that is similar to the D. melanogaster F-element genes. The green “M” below the coding span denotes the Methionine at the translation start site, and the red star denotes the stop codon.

Figure 7

Histone modification profiles for D. ananassae and D. melanogaster F-element genes at the third instar larval stage of development. (A) Metagene analysis shows that the region surrounding the 5′ end of F-element genes is enriched in H3K4me2 while the body of the coding span is enriched in H3K9me2. The values in the y-axis within each metagene plot correspond to the log-likelihood ratio between each ChIP sample and input control (assuming a dynamic Poisson model) as determined by MACS2. (B) Differences in the H3K27me3 enrichment patterns for the D. melanogaster ey gene and its ortholog in D. ananassae. The entire coding span of the ey gene is enriched in H3K27me3 in D. melanogaster (top) (for the D. melanogaster gene models, the thick boxes denote the coding exons and the thin boxes denote the untranslated regions). By contrast, only the region surrounding the 5′ end of the ey ortholog in D. ananassae shows H3K27me3 enrichment. The 5′ ends of the A and D isoforms of ey shows enrichment of H3K4me2 and H3K27me3 in both D. melanogaster and D. ananassae. These bivalent domains suggest that these two isoforms of ey are poised for activation at the third instar larval stage of development in both species.

Figure 8

D. ananassae F-element genes show similar expression patterns compared to genes on other Muller elements. RNA-Seq reads from seven samples (adult females, adult males, female carcass, male carcass, female ovaries, male testes, and embryos) were mapped against the improved D. ananassae genome assembly and the read counts for the Gnomon gene predictions were tabulated by htseq-count. The read counts for the seven samples were normalized by library size and then transformed using Tikhonov/ridge regularization in the DESeq2 package to stabilize the variances among the samples. The violin plots compare the distributions of the regularized log₂ expression values for the D. ananassae Gnomon gene predictions on all scaffolds (All), on the F-element scaffolds [F (all)], and on the base of the D element [D (base)] for these different developmental stages and tissues.