The Role of <i>piRNA</i>-Mediated Epigenetic Silencing in the Population Dynamics of Transposable Elements in <i>Drosophila melanogaster</i>

<div><p>The <i>piwi</i>-interacting RNAs (<i>piRNA</i>) are small RNAs that target selfish transposable elements (TEs) in many animal genomes. Until now, <i>piRNAs’</i> role in TE population dynamics has only been discussed in the context of their suppression of TE transposition, which alone is not sufficient to account for the skewed frequency spectrum and stable containment of TEs. On the other hand, euchromatic TEs can be epigenetically silenced via <i>piRNA</i>-dependent heterochromatin formation and, similar to the widely known “Position-effect variegation”, heterochromatin induced by TEs can “spread” into nearby genes. We hypothesized that the <i>piRNA</i>-mediated spread of heterochromatin from TEs into adjacent genes has deleterious functional effects and leads to selection against individual TEs. Unlike previously identified deleterious effects of TEs due to the physical disruption of DNA, the functional effect we investigated here is mediated through the epigenetic influences of TEs. We found that the repressive chromatin mark, H3K9me, is elevated in sequences adjacent to euchromatic TEs at multiple developmental stages in <i>Drosophila melanogaster</i>. Furthermore, the heterochromatic states of genes depend not only on the number of and distance from adjacent TEs, but also on the likelihood that their nearest TEs are targeted by piRNAs. These variations in chromatin status probably have functional consequences, causing genes near TEs to have lower expression. Importantly, we found stronger selection against TEs that lead to higher H3K9me enrichment of adjacent genes, demonstrating the pervasive evolutionary consequences of TE-induced epigenetic silencing. Because of the intrinsic biological mechanism of <i>piRNA</i> amplification, spread of TE heterochromatin could result in the theoretically required synergistic deleterious effects of TE insertions for stable containment of TE copy number. The indirect deleterious impact of <i>piRNA</i>-mediated epigenetic silencing of TEs is a previously unexplored, yet important, element for the evolutionary dynamics of TEs.</p></div

Differences in expression levels between alleles with and without TEs.

<p>1. Median for the differences in expression rank between alleles with and without TEs. Rank differences were calculated as the mean expression rank of “with TE” alleles minus that of “without TE” alleles. A positive rank difference means “with TE” alleles have larger expression rank (and thus lower expression) than “without TE” alleles.</p><p>2. The number of genes whose “with TE” alleles have significantly (permutation <i>p-value</i> < 0.05) larger expression rank (i.e. lower expression) than “without TE” alleles.</p><p>3. The expected number of genes whose “with TE” alleles have significantly larger expression rank (i.e. lower expression) than “without TE” alleles under the null hypothesis with 5% false positive rate.</p><p>Differences in expression levels between alleles with and without TEs.</p

Comparisons of the H3K9me3 density for genes whose nearest TEs have different population frequencies.

<p>TEs that are observed in the North American <i>D</i>. <i>melanogaster</i> population (Observed) have higher population frequencies than those that are not observed in the population (NOT). Genes whose nearest TEs are observed in the population have significantly lower H3K9me3 density. Dashed lines show the median of H3K9me3 density for genes whose nearest TE are not observed in the population. Numbers of genes included in the analysis are 297 (for Observed) and 491 (for NOT). Notations for <i>p-values</i> are * (<i>p</i> < 0.05), ** (<i>p</i> < 0.01), and *** (<i>p</i> < 0.001).</p

The associations between H3K9me3 density of genes and their neighborhood TE content.

<p>(A) Boxplots for the H3K9me3 density of genes that are of different distance from TEs. The H3K9me3 density of genes that have at least one TE in the introns of the gene (in gene), within 1kb, 1-2kb, 2-5kb, and 5-10kb upstream/downstream from the gene are significantly lower than those of genes that have no TE 10kb upstream/downstream (whose median of H3K9me3 density is shown as a dashed line). Numbers of genes in each category are 588 (in gene), 579 (within 1kb), 323 (1-2kb), 859 (2-5kb), 1162 (5-10kb), and 8072 (no TE in 10kb). All comparisons are against genes without TEs in 10kb upsteam/downstream. Result of 0–4 hr embryo is shown. (B) The <i>spearman rank</i> correlation coefficients between H3K9me3 density of a gene and the number of adjacent TEs within a specific window decrease as the distance between genes and TEs increases. The median, mean, and maximum number of TEs in specific windows are (1) in gene: 0 (median), 0.0523 (mean) and 8 (maximum), (2) 0-1kb: 0 (median), 0.0555 (mean), and 4 (maximum), (3) 1-2kb: 0 (median), 0.0341 (mean), 2 (maximum), (4) 2-5kb: 0 (median), 0.101 (mean), and 5 (maximum), and (5) 5-10kb: 0 (median), 0.158 (mean), and 7 (maximum). Notations for <i>p-values</i> are * (<i>p</i> < 0.05), ** (<i>p</i> < 0.01), and *** (<i>p</i> < 0.001).</p

The decay of H3K9me3 density of intergenic sequences adjacent to TEs.

<p>The median (A) and mean (B) of H3K9me3 density in nonoverlapping 1kb windows of intergenic sequences decreases as the windows are farther from TEs. Different colors are for different developmental stages: embryo 0–4 hr (light blue), embryo 4–8 hr (blue), embryo 8–12 hr (light green), embryo 12–16 hr (green), embryo 16–20 hr (pink), embryo 20–24 hr (red), L1 larvae (yellow), L2 larvae (orange), and pupae (purple). Dashed lines are the median/mean of the H3K9me3 density in the window closest to TEs. The observed median (C) and mean (D) H3K9me3 densities of windows adjacent to TEs are higher than those adjacent to 1,000 sets of randomly chosen TE-size sequences (gray lines), particularly for windows that are closer to TEs. Results of 0–4 hr embryos are shown here. Note that the scales of y-axis are different between (A) and (B) and between (C) and (D).</p

Results from the molecular evolution analysis of genes that occur in <i>D. melanogaster</i>.

<p>The dN/dS refers to a <i>D. melanogaster-D. simulans</i> pairwise calculation.</p

Delineating HP1E loss in the <i>obscura</i> group.

<p>We amplified the syntenic region of HP1E in the obscura group and successfully identified intact <i>HP1E</i> genes from <i>D. guanche</i> and <i>D. bifasciata</i>. We found highly pseudogenized versions of <i>HP1E</i> in <i>D. azteca</i>, <i>D. affinis</i> (<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002729#pgen.1002729.s003" target="_blank">Figure S3</a>). These latter four species also share dramatic karyotypic changes specific to this lineage including an <i>X:3L</i> fusion, a <i>Y:4</i> fusion and a neo-<i>Y</i> (indicated as <i>Y</i>′ in figure, note that <i>3L</i> and <i>4</i> = elements “D” and “F”, respectively). Thus, to the level of resolution possible from the available species, <i>HP1E</i> loss coincided with the karyotypic changes in the <i>obscura group</i>. The <i>HP1E</i> cytolocation on chromosome <i>3R</i> (element “E”), post-karyotype evolution, is apparently undisrupted.</p

HP1 diversity in Drosophila genomes.

<p>A. Phylogeny of 12 Drosophila species, which were each queried for HP1-like genes in this study <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002729#pgen.1002729-Clark1" target="_blank">[2]</a>, <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002729#pgen.1002729-Adams1" target="_blank">[42]</a>. Scale bar refers to the approximate divergence time between these species <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002729#pgen.1002729-Clark1" target="_blank">[2]</a>, <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002729#pgen.1002729-Adams1" target="_blank">[42]</a>. B. Schematics of proteins encoded by the various HP1 genes in Drosophila genomes are presented alongside the HP1 gene name. Highlighted in boxes are the canonical chromo (green) and shadow (blue) domains that typify HP1 genes. Note that in some instances, we were unable to confirm the exact gene model and therefore the lengths of the N-terminal tails (these are indicated with dashed lines). We also report the <i>D. melanogaster</i> cytolocation of the gene or if the gene is absent in <i>D. melanogaster</i>, the sytenic location in the <i>D. melanogaster</i> genome based on neighboring genes. The final column reports the species in which the gene is found. Genes shaded gray represent founding HP1 gene family members that were reported in the original <i>D. melanogaster</i> genome sequencing study <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002729#pgen.1002729-Adams1" target="_blank">[42]</a>. “*” refers to an allele that harbors a premature stop codon but conserved C-terminal sequence (<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002729#pgen.1002729.s007" target="_blank">Table S1</a>, <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002729#pgen.1002729.s004" target="_blank">Figure S4</a>) and predicted CD and CSD domains, consistent with a polymorphic full length gene or an incorrect base call.</p

PAML analysis results of genes that occur in <i>D. melanogaster</i>.

<p><i>mel</i> = <i>D. melanogaster</i>, <i>sim</i> = <i>D. simulans</i>, <i>sec</i> = <i>D. sechellia</i>, <i>yak</i> = <i>D. yakuba</i>, <i>ere</i> = <i>D. erecta</i>, <i>tak</i> = <i>D. takahashii</i>, <i>bia</i> = <i>D. biarmipes</i>, <i>ele</i> = <i>D. elegans</i>, <i>fic</i> = <i>D. ficusphilia</i>.</p

Results of the population genetic analyses of HP1 genes that occur in <i>D. melanogaster</i>.

<p>ENC = effective number of codons, Syn = synonymous, Nonsyn = nonsynonymous, %ile = percentile, ave n = average # alleles, NS = nonsynonymous, S = synonymous, π ratio = NSπ/Sπ, poly = #polymorphisms, fix = # fixations, NI = neutrality index, FETpval = Fisher's Exact Test probability value.</p