Rapid Evolution of piRNA Pathway in the Teleost Fish: Implication for an Adaptation to Transposon Diversity

The Piwi-interacting RNA (piRNA) pathway is responsible for germline specification, gametogenesis, transposon silencing, and genome integrity. Transposable elements can disrupt genome and its functions. However, piRNA pathway evolution and its adaptation to transposon diversity in the teleost fish remain unknown. This article unveils evolutionary scene of piRNA pathway and its association with diverse transposons by systematically comparative analysis on diverse teleost fish genomes. Selective pressure analysis on piRNA pathway and miRNA/siRNA (microRNA/small interfering RNA) pathway genes between teleosts and mammals showed an accelerated evolution of piRNA pathway genes in the teleost lineages, and positive selection on functional PAZ (Piwi/Ago/Zwille) and Tudor domains involved in the Piwi–piRNA/Tudor interaction, suggesting that the amino acid substitutions are adaptive to their functions in piRNA pathway in the teleost fish species. Notably five piRNA pathway genes evolved faster in the swamp eel, a kind of protogynous hermaphrodite fish, than the other teleosts, indicating a differential evolution of piRNA pathway between the swamp eel and other gonochoristic fishes. In addition, genome-wide analysis showed higher diversity of transposons in the teleost fish species compared with mammals. Our results suggest that rapidly evolved piRNA pathway in the teleost fish is likely to be involved in the adaption to transposon diversity.

A substantial amount of piRNAs is derived from TEs, which play a crucial role in silencing TEs (Aravin et al. 2007). TEs are selfish genetic elements that can mobilize within a genome, potentially leading to serious genome damage. piRNAs act as guardians of the genome and protect it from invasive TEs in the germline cell (Kalmykova et al. 2005). In the process of gametogenesis, derepression of even a single active TE may cause infertility, a disaster from a Darwinian point of view.
The types of transposons vary dramatically in vertebrate species, ranging from a few to over a hundred (Malone and Hannon 2009). In Drosophila, there are approximately 150 different TE types, which have different expression, replication, and mobilization strategies. There are much more TE types in teleost fishes than that in mammals, whereas some types are absent from the mammalian genomes (Aparicio et al. 2002;Volff et al. 2003). Around 20 clades of retrotransposons are present in pufferfish genomes, whereas only six clades were detected in human and mouse genomes (Volff et al. 2003). Notably, fish genomes contain over 30 types of LINE-1 retrotransposons, each of which is active, whereas mammalian genomes possess only one (Furano et al. 2004). About 4.8% transcriptome of the platyfish is derived from TE sequences which represent 40 different families, suggesting that many TEs are still active (Schartl et al. 2013). DNA transposons in fish species are often active, whereas they are inactive in mammals (Bohne et al. 2012). Although TE family numbers are lower in mammals than in fish species, TE copy numbers are generally much higher in mammals than in fishes (Volff 2005). Therefore, TEs in fish genomes appear to have undergone a higher turnover than in the mammalian genomes , which is similar to that observed in Drosophila (Eickbush and Furano 2002).
In mouse, piRNAs are mainly derived from three classes of retrotransposons: LTR, LINE, and SINE (Aravin et al. 2007(Aravin et al. , 2008. However, zebrafish piRNAs have been observed to match to both retrotransposons and DNA transposons (Houwing et al. 2007(Houwing et al. , 2008 and are divided into more than 20 groups (Zhou et al. 2010). Diverse piRNAs may provide an adaptive defense to repress diversified transposons in teleosts. These observations prompt us to test whether high diversity and turnover rate of TEs in fish have driven the evolution of piRNA pathway genes, which is involved in TE silencing.
On the basis of their variety and abundance, TEs present an imposing threat to the germline cells and their genomes. With regard to the germline cells, they have a special need to protect their genome, because genetic information must be faithfully transmitted to offspring. Thus, it is crucial to repress TE activity to ensure the integrity of the germline genomes. As a result, piRNA machinery was suggested to have undergone adaptive evolution, indicating that TE activity has driven the evolution of piRNA pathway (Malone and Hannon 2009). Thus it seems that as a result from a coevolutionary arms race, there is a process of reciprocal adaptation between piRNA machinery and TE activities. Several recent studies observed evidences of selection in the piRNA pathway in Drosophila (Vermaak et al. 2005;Castillo et al. 2011;Kolaczkowski et al. 2011;Simkin et al. 2013). So far, studies on adaptive evolution of piRNA pathway have been largely confined to the Drosophila species. As such, adaptive evolution of piRNA pathway remains unknown in vertebrates.
To test a possible association between piRNA pathway evolution and TE diversity in the teleost fish species, we conducted a comparatively evolutionary analysis for piRNA pathway genes between the teleost fish species and mammals, and compared these genes with Ago genes that function in miRNA/siRNA pathway. The swamp eel (Monopterus albus) that belongs to the order Synbranchiformes is a protogynous fish, which can transform its sex from female via intersex into male (Zhou et al. 2001;Cheng et al. 2003). We tested whether the selective pressure of piRNA pathway differs between swamp eel (a kind of hermaphroditic fish) and the other fish species. We found an accelerated evolution of piRNA pathway genes in the teleost lineage. Notably, we detected many positively selected sites distributed on functional domains of piRNA pathway genes on the ancestral branch of teleosts, suggesting adaptive evolution of piRNA pathway genes as organisms evolve. In addition, we observed that five piRNA pathway genes in the swamp eel evolved more rapidly than other fish species, and their high expression level in three types of gonads (ovary, ovotestis, and testis) during sex reversal suggests the crucial role of piRNA pathway genes during the sex reversal. To further test whether high diversity and turnover rate of TEs exist in fish, we observed various types of fish TEs in comparison with mammals, suggesting a link of TEs diversity to the rapid evolution of piRNA pathway in fish species. Thus, this study highlights that rapidly evolved piRNA pathway is likely to adapt to transposon diversity in the teleost fish lineages.

Phylogenetic Analysis
The protein sequences of 30 vertebrates were aligned with ClustalW program (Chenna et al. 2003) and manual adjustments. The maximum-likelihood tree was constructed using PhyML 3.0 (Guindon and Gascuel 2003;Guindon et al. 2010) with LG model (Le and Gascuel 2008). To root the tree, lamprey Piwil1 was used as an outgroup. The phylogeny was analyzed using 100 replicates. The neighbor joining (NJ) trees using bootstrap analysis with 1,000 replicates were constructed using MEGA 5 (Tamura et al. 2011) with Poisson model.

Synonymous and Nonsynonymous Distance
The coding sequences of 18 genes were aligned separately with ClustalW (Chenna et al. 2003) and manual adjustments. The sequences of each gene from both mammals and teleosts were divided into two groups (one group: 21 mammals, and the other: 9 teleost fish species). Synonymous (d S ) and nonsynonymous (d N ) were calculated using the modified Nei-Gojobori method (Zhang et al. 1998) implemented in MEGA 5 (Tamura et al. 2011). The o ratios (d N /d S ) were calculated to compare the evolutionary distance for each gene between two groups.

Branch Model Test
To test various selection pressure acting on different clades, we used branch model based on the o ratio (d N /d S ) in the codeml in PAML 4 package (Yang 2007). A species tree of 30 vertebrates (from Ensembl and teleost phylogeny; Chen and Mayden 2010) was used to determine the evolutionary relationship. We first tested the one-ratio model as the null hypothesis, in which the entire tree has one o ratio. Two-ratio model was then used (a background ratio o 0 and a different ratio o 1 on the foreground branch). The significance for the likelihood ratio tests (LRTs) was tested using w 2 analysis. 2Áln L was used (twice the difference of log likelihood between two models) in this analysis.

Detection of Positive Selection
To detect positively selected sites, the branch-site model in the program codeml in PAML 4 was used (Yang 2007). We set the ancestral branch of teleosts as foreground and the other lineages as background branch. The test 2 in PAML 4 was used to compare two models. In model A, three o ratios Null model parameters are the same as in model A, whereas o 2 ¼ 1. LRTs were calculated to test the significance of differences between two models. If the test shows positive selection, the sites under selection were analyzed according to high Bayes empirical Bayes posterior probabilities (BEB PPs) (Yang et al. 2005). All positively selected sites were aligned to the protein sequences and functional domains in the Pfam database (Coggill et al. 2008;Punta et al. 2012).

Sliding Window Analysis
To analyze distribution of all sites with the o values (d N /d S ) across the Piwil1 sequence, we performed sliding window analysis of the o values of each interval estimated by SWAAP (Pride 2000) with window size 60 bp and step size 9 bp. The Piwil1 sequences from 20 mammals and 8 teleosts were used in the analysis.

Annotation of Repeat-Associated piRNAs
The zebrafish piRNA sequences were from Ketting's study (Houwing et al. 2007) together with our data (Zhou et al. 2010). Nonredundant piRNA sequences were used to search against TE superfamilies using Repeat-Masking (Kohany et al. 2006) with default parameters in Repbase (Jurka 2000). Percentage for each TE superfamily was calculated with respect to total TE-associated piRNAs. A total of 97,666 zebrafish piRNAs were analyzed using the annotated TE database in Repbase.

piRNA Pathway Genes Evolve More Rapidly in Teleosts than in Mammals
To gain insights into evolutionary dynamics of piRNA pathway in vertebrate phylogeny, we first performed phylogenetic analysis of Piwil1 protein, which is a core component of piRNA machinery in vertebrates. Comparative analysis of amino acid substitution rates among 30 vertebrate species, including 21 mammals and 9 teleost fish species, showed longer branch length of the Piwil1 in teleosts than in mammalian species (similar results in other piRNA pathway proteins were also observed), suggesting higher substitution rate in fish branches ( fig. 1). To further investigate accelerated evolution of piRNA pathway in teleosts, using the modified Nei-Gojobori method (Zhang et al. 1998), we compared the synonymous and nonsynonymous nucleotide substitution distances of 12 piRNA pathway genes together with 5 miRNA/siRNA pathway genes between two groups (mammals vs. teleosts). The d N value differences between two groups were markedly higher in piRNA genes (except for Hen1) than in Ago genes ( fig. 2A). This suggested that piRNA pathway genes in teleosts undertook more nonsynonymous substitutions than in mammals. To analyze evolutionary rates in piRNA pathway genes between mammals and teleosts, we computed the d N /d S values for each gene in two groups ( fig. 2B). We observed that d N /d S ratios of piRNA pathway genes were higher than that of Ago genes (P < 0.01), suggesting piRNA pathway genes evolved more rapidly than miRNA/ siRNA pathway genes. In addition, the d N /d S ratios of 12 piRNA pathway genes except for Hen1 in teleosts were markedly higher than that in mammals. These results suggested that evolutionary rates of piRNA pathway genes were higher in teleosts than in mammals. Ago3b, a teleost-specific duplicated gene of Ago3 (McFarlane et al. 2011), is an exception. Gene duplication leads to a relaxation of selection and consequently an elevated rate of molecular evolution for the duplicated genes (Kondrashov et al. 2002). Therefore, the d N /d S ratio of Ago3b in teleosts is higher than in mammals.
To confirm these results obtained from distance analysis, we further performed branch model analysis by the codeml algorithm from the PAML 4 package (Yang 2007).
Comparative analysis between the two-ratio model and one-ratio model showed that o 1 (d N /d S ) values for all piRNA pathway genes except Hen1 in the teleost lineage were significantly larger than o 0 in mammals (table 1). In contrast, o 1 for Ago2, Ago3a, and Ago4 in the teleost lineage were significantly smaller than o 0 in the mammal lineages. The o of Ago1 and -actin have no significant difference between mammals and fishes. These results from branch model analysis were consistent with synonymous and nonsynonymous distance analysis. These data suggested that piRNA pathway genes evolved more rapidly in teleosts than in mammals.

Adaptive Evolution in Functional Domains of piRNA Pathway Genes
Although o ratios of piRNA pathway genes were smaller than 1 at whole gene level, it cannot be judged that positive selection did not occur in piRNA pathway genes. To examine the possible positive selection occurred in piRNA pathway, we used the branch-site model (Yang 2007) to test whether there was positive selection on the ancestral branch of the teleost lineages. Branch-site model test showed evidence of positive selection for most piRNA pathway genes (table 2). There were more positively selected sites in Piwi proteins and TDRDs proteins than other piRNA pathway proteins, for example, 13, 12, 9, and 7 positively selected sites with BEB PPs higher than 95% were detected in Piwil1, Piwil2, Tdrd1, and Tdrkh, respectively (table 2).
The ancestral branch leading to the teleosts was assigned as the foreground branch. Twice the difference between the likelihood values of the null model and the alternative model is indicated as 2ÁL. The modified model A is used as the alternative model and the modified model A with o 2 fixed at 1 is the null model. b P values significant at the 5% threshold after Bonferroni correction in bold. c P 0 is the proportion of codons that have o 0 in all branches, P 1 is the proportion of codons that have o 1 ¼ 1 in all branches, P 2a is the proportion of codons that have o 0 in the background branches but o 2 in the foreground branches, and P 2b is the proportion of codons that have o 1 in the background branches but o 2 in the foreground branches. d Sites with the BEB PPs higher than 90% are shown. Sites higher than 95% are single underlined and sites higher than 99% are double underlined. The sites of Mov10l1 and Fkbp6 are indexed by the amino acids at the sites in the tetraodon, whereas the sites of the other genes are indexed by the amino acids at the sites in the zebrafish.
We further analyzed locations of the positively selected sites in these genes. Using the protein topologies, we observed that their distributions were clearly concentrated on functional domains ( fig. 3A). As for Piwi proteins, 50% (Piwil1) and 56% (Piwil2) positively selected sites were distributed in PAZ-MID (Piwi/Ago/Zwille-Middle) domain, which is responsible for recognition of piRNA 3 0 -/5 0 -end (Tian et al. 2011). In addition, 31% (Piwil1) and 22% (Piwil2) positively selected sites were distributed in PIWI domain, which is similar to ribonuclease H (Song et al. 2004;Tolia and Joshua-Tor 2007). Branch-site results were further supported by sliding window analysis, which showed variation in o values along the Piwil1 gene ( fig. 3B). These data suggested that positive selection was strongest at the interface between Piwi proteins and target piRNA molecules, implicating their possible roles in modulating interaction of Piwi proteins with diverse piRNAs in teleosts.
As for two TDRDs proteins, 50% and 30% positively selected sites were distributed in Tudor domain of Tdrd1 and Tdrkh, respectively. Tudor domain can interact with sDMA in the N-terminal region of Piwi proteins Mathioudakis et al. 2012). In addition, 44% positively selected sites were distributed in Tdrkh-KH domain, which is an RNAbinding domain (Siomi et al. 1994;Valverde et al. 2008). The adaptive substitutions in the Tudor domain suggested their possible roles in modulating protein-protein interaction in the piRNA pathway. In contrast, LRTs were not significant and no sites have been detected in Ago2, Ago3a, Ago3b, and Ago4, whereas only two positively selected sites have been detected in Ago1.
3D structure homology modeling further showed that some positively selected sites in Piwil1-PAZ and Piwil2-PAZ were distributed on pocket of Piwi-PAZ for binding piRNA 3 0 -end, which may facilitate piRNA recognition ( fig. 4A and  B). Some positively selected sites in Tdrd1-third Tudor domain and Tdrkh-Tudor domain were localized in the pocket of Tudor for binding Piwi-sDMA ( fig. 4C and D) to modulate Piwi-Tudor protein interaction in the piRNA pathway. These analysis highlighted importance of amino acid changes in PAZ domain and Tudor domain, which are involved in proteinprotein/RNA interaction. Overall, these results suggested that evolution of functional domains in piRNA pathway genes is likely to be adaptive to their functions in piRNA pathway in the teleost fish species.

Higher Substitution Rate in piRNA Pathway Genes in the Swamp Eel
The longest branch in the swamp eel (hermaphroditic fish) Piwil1 as observed in the teleost lineages ( fig. 1) indicated a higher substitution in swamp eel than in other fish species. To test possible accelerated evolution in this species, we first cloned five core piRNA pathway genes (Piwil1, Piwil2, Tdrd1, FIG. 6.-Association between rapid evolution of piRNA pathway and TEs diversity in the teleost fish species. (A) Comparison of TE superfamilies with copy number >1,000 between teleost fish species and mammals. Blue columns, retrotransposon superfamily numbers (Class I); red columns, DNA transposon superfamily numbers (Class II). The data of human and mouse are from a previous study results (Lander et al. 2001;Mouse Genome Sequencing Consortium et al. 2002). As many as 12-25 superfamilies of retrotransposons have been observed in the teleost fish species, whereas only 10-12 superfamilies are detected in the genomes of both human and mouse. For Class II, 7-18 superfamilies are detected in the fish species, whereas there are 2-3 Tdrkh, and Prmt5) and five Ago genes (Ago1, Ago2, Ago3a, Ago3b, and Ago4) from swamp eel using degenerate PCR and RACE technique. RT-PCR showed that piRNA pathway genes were mainly expressed in three kinds of gonads (ovary, ovotestis, and testis), whereas Ago genes were expressed in all adult tissues detected, indicating a role of piRNA pathway genes in the germline ( fig. 5A). We then examined evolutionary rate of piRNA pathway genes in the species. The branch model analysis (one-ratio model vs. two-ratio model) showed that five piRNA pathway genes (Piwil1, Piwil2, Tdrd1, Ddx4, and Prmt5) in the swamp eel evolved more rapidly than that in other fish species significantly (fig. 5B). These results suggested that piRNA pathway genes in the swamp eel evolved at an elevated rate compared with the other gonochoristic fishes.

Diversity of TEs and TE-Derived piRNAs in Teleosts
To connect rapid evolution of piRNA pathway genes with TEs diversity in the teleost fish species, using newly released whole genome sequence data, together with new genome (swamp eel), TEs in six teleost fish species were analyzed. There were 12-25 superfamilies of retrotransposons (Class I TEs) detected in different fish species ( fig. 6A), whereas only 10-12 superfamilies in human and mouse genome. For Class II, 7-18 superfamilies were detected in the fish species, whereas there were 2-3 superfamilies in both human and mouse (Lander et al. 2001;Mouse Genome Sequencing Consortium et al. 2002). Thus, the teleost fish genomes have much higher TE diversity compared with mammalian genomes, suggesting a correlation between TEs diversity and the rapid evolution of piRNA pathway. Further analysis of loci distribution of piRNAs in zebrafish genome showed that TE-derived piRNAs were corresponding to 16 superfamilies of retrotransposons and 13 superfamilies of DNA transposons ( fig. 6B). However, most TE-derived piRNAs were originated from about three superfamilies of retrotransposons in mouse (Aravin et al. 2007). In addition, DNA trsansposons (Class II) are inactive in mammals (Huang et al. 2012). These data suggested an association between rapid evolution of piRNA pathway and TEs diversity in the teleost fish species.

Discussion
Advantageous mutation accumulation through positive selection may result in adaptive evolution. Adaptive evolution often drives rapid evolution (Swanson et al. 2001). Genes involved in some biological processes often evolve rapidly. For example, immune-related genes are generally thought as a class of rapidly evolving genes, because they are often involved in hostparasite coevolutionary arms race. In fact, many reproductive genes, such as those involved in sex determination (Whitfield et al. 1993;Zhang 2004) and gamete recognition (Swanson et al. 2001), also evolve rapidly as a result of adaptive evolution. In this report, we present that adaptive evolution may have shaped piRNA pathway in the teleost fish lineage, and the rapid evolution of piRNA pathway is consistent with the diversity of TEs, which support a correlation between piRNA pathway evolution and transposon diversity in the teleost fish. Based on these studies, a model of piRNA pathway in silencing diverse TEs in teleost fish genome is proposed ( fig. 6C). Diverse piRNAs transcribed from corresponding TE loci are loaded onto Piwi proteins, which have undergone an adaptive evolution in the lineage of the teleost fish. Besides Piwi proteins, other factors, such as Tdrd1 and Tdrkh, are also involved in piRNA pathway. Finally, piRNA pathway represses TE activation through DNA methylation to protect genome in germline cells.
Studies in Drosophila have indicated that RNAi genes evolve rapidly (Obbard et al. 2006;Kolaczkowski et al. 2011). Because RNAi processing is associated with TEs silencing, the synergetic effects between RNAi and TEs should be important in maintenance of genome stability through adaptive evolution between small RNA pathway and TEs. Our observations in the teleost fish lineages further suggest that rapid evolution of fish piRNA pathway is likely to be involved in the adaption to transposon diversity. Rapid evolution of piRNA pathway in fish may silence TEs more efficiently, which is consistent with observations in the Drosophila (Fablet et al. 2014). We find that positive selection has occurred in the teleost lineages with strongest at the interface between silencing machinery and target piRNAs. These results extend our understanding of piRNA pathway adaptation to ensure genome stability.
piRNA pathway genes play crucial roles in the transposon repression and the germline functions. Rapid evolution of  (Houwing et al. 2007) together with our data (Zhou et al. 2010) were used to search against TE superfamilies using Repeat-Masking with default parameters in Repbase. The number indicates percentage for each superfamily with respect to total TE-associated piRNAs. Different colors in pie diagram represent different kinds of the TE superfamily. TE superfamilies with percentages less than 1% are included in other TEs. DIRS, named for the founding member from Dictyostelium discoideum; ERV, endogenous retrovirus; LTR, long terminal repeat; L1, long interspersed element 1; L2, long interspersed element 2; non-LTR, nonlong terminal repeat; SINE, short interspersed elements. (C) A model of adaptive evolution of piRNA pathway to silence diverse TEs in fish genome. Diverse piRNAs transcribed from corresponding TE loci are loaded onto Piwi proteins, which have undergone an adaptive evolution in the lineage of the teleost fish. Besides Piwi proteins, other factors, such as Tdrd1 and Tdrkh, are also involved in piRNA pathway. Finally, piRNA pathway represses TE activation through DNA methylation to protect genome in germline cells. piRNA pathway genes in the teleost fish lineages is consistent with observations in the Drosophila species (Simkin et al. 2013), supporting an evolutionary "arms race" between TEs and the piRNA pathway. Moreover, key functional domains of Piwi proteins we detected may serve as hotspots of adaptive evolution. The adaptation at MID and PAZ domains of Piwil1 and Piwil2 provides access to their diverse target piRNA molecules, which is consistent with much more diversified TEs in the teleost genomes. As an epigenetic regulation, the piRNA/TE system is important for transposon silencing to ensure integrity of animal genomes. In addition, other epigenetic modes are also implicated in TE regulation, such as chromatin structure and DNA methylation (Puszyk et al. 2013). All these epigenetic regulations play crucial roles in genome defense.
A distinct biological feature that may shape susceptibility to TEs is fertilization mode (Huang et al. 2012). The mammals reproduce by internal fertilization, an adaptation to live on land rather than in water as teleosts do. In contrast, most teleosts usually reproduce by external fertilization. As both sperm and eggs are mobile in a water environment, the fertilization in water is subjected to pressure from harsh aquatic environment, climate change, and water pollution. In order to ensure sufficient fertilized eggs for species survival, external fertilization animals definitely have to shed much more sperms and eggs into water than internal fertilization animals do. Thus, gametogenesis is very active so as to maintain sustained reproductive capacity, giving rise to higher risk of transposon jumping in genomes at the same time. The piRNA pathway genes, as guardians of the genome playing crucial role in silencing transposons, are supposed to be subjected to strong and persistent selective pressure in the teleosts. A typical example of the rapid evolution in the piRNA pathway is observed in the hermaphroditic fish, swamp eel; piRNA pathway genes evolved at an elevated rate compared with the other gonochoristic fish species. In response to selective pressures from diversified TEs and reproductive mode of external fertilization, piRNA pathway genes in the teleosts are supposed to evolve faster than their orthologous genes in the mammals.

Conclusions
The piRNA pathway is responsible for germline specification, gametogenesis, transposon silencing, and genome integrity. Mobile elements can disrupt genome and its functions. piRNA pathway evolution and its adaptation to transposon diversity in the teleost fish remain unknown. In this report, we present that adaptive evolution is likely to shape piRNA pathway in the teleost fish lineages, and the rapid evolution of piRNA pathway is associated with the diversity of TEs, which suggested a link of adaptation evolution of piRNA pathway to transposon diversity in the teleost fish.