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.

Introduction

The Piwi-interacting RNA (piRNA) pathway is a small RNA silencing system, which functions in germline specification, gametogenesis, transposon silencing, genome integrity, and stem cell maintenance across the animal phylogeny. piRNAs are endogenous small noncoding RNA molecules ranging from 26 to 31 nucleotides in length. Animals lacking piRNAs exhibit activation of transposable elements (TEs) and defects in gametogenesis (Carmell et al. 2007). Piwil1 (Piwi-like 1) and Piwil2 (Piwi-like 2) are core proteins interacted with piRNAs (Aravin et al. 2006; Girard et al. 2006; Houwing et al. 2007, 2008) and act to repress TEs both transcriptionally and posttranscriptionally. Other components involved in piRNA pathway have also been identified in both vertebrates and Drosophila (Ishizu et al. 2012). TDRDs proteins (Tudor domain-containing proteins) can associate with both Piwil1 and Piwil2 through their N-terminal symmetrical dimethyl arginines (sDMAs) in both mouse and zebrafish (Chuma et al. 2006; Chen et al. 2009; Reuter et al. 2009; Huang, Houwing, et al. 2011). The protein arginine methyltransferase 5 (Prmt5) catalyzes methylation of the arginines of the sDMAs (Kirino et al. 2009; Vagin et al. 2009). Mitochondrial protein Pld6 (Phospholipase D member 6, as known as MitoPLD) is also involved in primary piRNA biogenesis in the mouse germline (Huang, Gao, et al. 2011; Watanabe et al. 2011). Genetic analyses have revealed that the mouse Ddx4 (DEAD box polypeptide 4, a.k.a. Vasa, Mvh) is essential for the piRNA ping-pong cycling (Kuramochi-Miyagawa et al. 2010). Asz1 (or Gasz) (Ankyrin repeat, SAM, and basic leucine zipper domain containing 1) plays a role in retrotransposon silencing by stabilizing Piwil2 (MILI) in nuage (Ma et al. 2009). Both Mov10l1 (Moloney leukemia virus 10-like 1) and Fkbp6 (FK506 binding protein 6) are involved in delivering piRNAs to the mouse Piwi proteins (Frost et al. 2010; Zheng et al. 2010; Xiol et al. 2012). Kinesin motor protein Kif17 (Kinesin family member 17) interacts with the mouse Piwil1 and acts in the nucleocytoplasmic transport in the chromatoid bodies (Kotaja et al. 2006). In addition, Hen1 (HEN1 methyltransferase homolog 1) is responsible for piRNA stability, through methylation of piRNA 3′-end of both mouse and zebrafish (Kirino and Mourelatos 2007; Kamminga et al. 2010). These piRNA pathway genes, which are evolutionarily conserved across the animal phylogeny, act jointly to protect the genome from invasive TEs in the germline (Ishizu et al. 2012).

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 (Duvernell et al. 2004), 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, 2008). However, zebrafish piRNAs have been observed to match to both retrotransposons and DNA transposons (Houwing et al. 2007, 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.

Materials and Methods

Data Preparation

The coding sequences of 18 genes (Piwil1, Piwil2, Tdrd1, Tdrkh, Pld6, Ddx4, Asz1, Mov10l1, Fkbp6, Kif17, Prmt5, Hen1, Ago1, Ago2, Ago3a, Ago3b, Ago4, and β-actin) in 29 vertebrates (21 mammals and 8 fishes) were obtained from Ensembl (www.ensembl.org, last accessed May 31, 2014) and National Center for Biotechnology Information (NCBI; http://www.ncbi.nlm.nih.gov, last accessed May 31, 2014) (supplementary table S1, Supplementary Material online). Five miRNA/siRNA pathway genes (Ago1, Ago2, Ago3a, Ago3b, and Ago4), belonging to Ago subfamily in the Argonaute family, and β-actin were used as control genes. Two analyses were performed to ensure orthology among genes. First, we blasted these genes in GenBank (http://www.ncbi.nlm.nih.gov, last accessed May 31, 2014) to ensure the high hits. Second, ClustalW is used to align the orthologous genes (Chenna et al. 2003) and to check their conservation among vertebrate species. A few coding sequences in some species (supplementary table S1, Supplementary Material online) are not available in Ensembl or NCBI because of unfinished genome sequencing. Some partial sequences with large gaps or sequencing errors were discarded in this study.

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 (dS) and nonsynonymous (dN) were calculated using the modified Nei–Gojobori method (Zhang et al. 1998) implemented in MEGA 5 (Tamura et al. 2011). The ω ratios (dN/dS) 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 ω ratio (dN/dS) 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 ω ratio. Two-ratio model was then used (a background ratio ω0 and a different ratio ω1 on the foreground branch). The significance for the likelihood ratio tests (LRTs) was tested using χ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 ω ratios (0 < ω0 < 1, ω1 = 1, ω2 > 1) are set for foreground and two ω ratios (0 < ω0 < 1, ω1 = 1) for background. Null model parameters are the same as in model A, whereas ω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 ω values (dN/dS) across the Piwil1 sequence, we performed sliding window analysis of the ω 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.

Homology Modeling

Three-dimensional (3D) protein structures of four zebrafish proteins (Piwil1-PAZ domain, Piwil2-PAZ domain, Tdrd1-the third Tudor domain, Tdrkh-Tudor domain) were built using the SWISS-MODEL server (http://swissmodel.expasy.org, last accessed May 31, 2014) (Schwede et al. 2003) with solved 3D homologous proteins as templates. The crystal structures of templates were obtained from the Protein Data Bank (PDB). The PDB ID (identity) are as follows: 3O7V (human Piwil1-PAZ) (Tian et al. 2011), 3O7X (human Piwil2-PAZ) (Tian et al. 2011), 4B9W (mouse Tdrd1-the third Tudor domain) (Mathioudakis et al. 2012), 3FDR (human Tdrkh-Tudor) (Chen et al. 2009). The models were visualized by Swiss-PDB Viewer v3.7 and rendered with POV-Ray v3.6 (Persistence of Vision Ray Tracer, USA). The positively selected sites distributed on these domains were highlighted as red spheres to show their locations in the protein structures.

Degenerate Polymerase Chain Reaction, Rapid Amplification of cDNA Ends, and Reverse Transcription Polymerase Chain Reaction

The swamp eels (M. albus) were obtained from markets in Wuhan, China. Their sexes were confirmed by microscopic analysis of gonad sections (ovary, ovotestis, and testis). Total RNAs were isolated from adult tissues using TRIZOL (Invitrogen, CA). SMART cDNAs were reversely transcribed from the RNAs using the SMART cDNA synthesis kit following protocol (Clontech, CA). Degenerate PCR was used to amplify the conserved regions of ten genes (Piwil1, Piwil2, Tdrd1, Tdrkh, Prmt5, Ago1, Ago2, Ago3a, Ago3b, and Ago4) using multiple degenerate primers (supplementary table S2, Supplementary Material online). PCR was performed in a 20 µl reaction mix containing 10 mM Tris–HCl, pH 8.3, 1.5 mM MgCl2, 50 mM KCl, 150 mM dNTP, 0.2 µM each primer, and 1 unit of Taq DNA polymerase. Based on the conserved regions sequences, 3′-RACE of these genes was performed using common primer CDSIII: 5′-ATTCTAGAGGCCGAGGCGGCCGACATGd(T)30N_1N-3′ (N = A, G, C, or T; N_1 = A, G, or C) and gene-specific primers (supplementary table S3, Supplementary Material online). 5′-RACE of these genes were performed using common SMARTIII primer: 5′-AAGCAGGGTTCAACGCAGAGTGGCCATTACGGCCGGG-3′ and gene-specific primers (supplementary table S3, Supplementary Material online). All sequences were cloned and sequenced. Semiquantitative RT-PCR was performed for five Ago genes (Ago1, Ago2, Ago3a, Ago3b, and Ago4) and six piRNA pathway genes (Piwil1, Piwil2, Tdrd1, Tdrkh, Ddx4, and Prmt5). Hprt was used as a control. The cDNAs from adult tissues were amplified using gene-specific PCR primers (supplementary table S4, Supplementary Material online).

TE Analysis

To analyze diversity of TEs in fish genomes, we combined homology-based and de novo approaches to identify TEs in genomes of six teleost fish species. Zebrafish (Danio rerio), tetraodon (Tetraodon nigroviridis), fugu (Takifugu rubripes), medaka (Oryzias latipes), and stickleback (Gasterosteus aculeatus) genome data were downloaded from Ensembl release 63. Swamp eel (M. albus) genome (DDBJ/EMBL/GenBank under the accession AONE00000000) was sequenced in our laboratory recently. The homology-based approach utilized database Repbase (release 16.02) with RepeatMasker (version 3.2.9) (Tarailo-Graovac and Chen 2009). The de novo approach utilized two prediction programs (RepeatModeler version 1.0.4 [Price et al. 2005] and LTR-FINDER version 1.0.5) to build the de novo repeat libraries based on the genome sequences. The multicopy genes and contaminations were removed from the libraries. Then, the RepeatMasker was used again to find repeats in these repetitive sequence libraries. Finally, we combined all the results generated by these methods and analyzed the diversity of TEs in fish genomes.

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.

Results

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 dN 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 dN/dS values for each gene in two groups (fig. 2B). We observed that dN/dS 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 dN/dS 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 dN/dS ratio of Ago3b in teleosts is higher than in mammals.

F

Phylogenetic tree of Piwil1 in 30 vertebrates. Phylogenetic analysis was performed with PhyML and the lamprey (an order of jawless fish, also a very ancient lineage of vertebrates) Piwil1 was used as an outgroup. Bootstrap values (percentage) were estimated using bootstrap analysis with 100 replicates and the values above 50% are shown next to the nodes. The number of amino acid substitutions per site was used to assess evolutionary distances. The distance scale (0.1) is shown. The accession numbers of all sequences in 30 vertebrates are listed in supplementary table S1, Supplementary Material online.

F

Nonsynonymous (dN) and synonymous (dS) nucleotide distance, and dN/dS ratio for 18 genes between teleosts and mammals. The first 12 genes (Piwil1, Piwil2, Tdrd1, Tdrkh, Pld6, Ddx4, Asz1, Mov10l1, Fkbp6, Kif17, Prmt5, and Hen1) are involved in piRNA pathways. Five Ago genes (Ago1, Ago2, Ago3a, Ago3b, and Ago4), which are involved in miRNA/siRNA pathways, and β-actin were used as control genes. (A) Nonsynonymous (dN) and synonymous (dS) nucleotide distance between teleosts and mammals. Straight lines between two circles (dS) or two triangles (dN) show distance between teleosts and mammals. All values are mean ± SE. (B) dN/dS ratio between teleosts and mammals. Black and gray columns show dN/dS ratios of different genes in teleosts and mammals, respectively. Gene full names are as follows: Piwil1 (Piwi-like 1), Piwil2 (Piwi-like 2), Tdrd1 (Tudor domain containing 1), Tdrkh (Tudor and KH domain containing protein), Pld6 (Phospholipase D member 6), Ddx4 (DEAD box polypeptide 4), Asz1 (Ankyrin repeat, SAM and basic leucine zipper domain containing 1), Mov10l1 (Moloney leukemia virus 10-like 1), Fkbp6 (FK506 binding protein 6), Kif17 (Kinesin family member 17), Prmt5 (Protein arginine methyltransferase 5), Hen1 (HEN1 methyltransferase homolog 1), Ago1 (Argonaute RISC catalytic component 1), Ago2 (Argonaute RISC catalytic component 2), Ago3a (Argonaute RISC catalytic component 3a), Ago3b (Argonaute RISC catalytic component 3b), and Ago4 (Argonaute RISC catalytic component 4). The dN/dS ratios of piRNA pathway genes are significantly larger than those of miRNA/siRNA pathway genes (P < 0.01, t-test).

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 ω1 (dN/dS) values for all piRNA pathway genes except Hen1 in the teleost lineage were significantly larger than ω0 in mammals (table 1). In contrast, ω1 for Ago2, Ago3a, and Ago4 in the teleost lineage were significantly smaller than ω0 in the mammal lineages. The ω 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.

Table 1

Branch model tests for piRNA pathway genes and Ago genes

Gene	2ΔLa	P Valueb	ω0c	ω1d
Piwil1	712.356	6.150E-157	0.022	0.165
Piwil2	48.940	2.639E-12	0.081	0.137
Tdrd1	8.499	0.003	0.253	0.300
Tdrkh	33.835	5.999E-09	0.126	0.210
Pld6	52.876	3.552E-13	0.036	0.131
Ddx4	76.679	2.011E-18	0.074	0.175
Asz1	35.661	2.348E-09	0.135	0.283
Mov10l1	60.794	6.336E-15	0.099	0.179
Fkbp6	31.633	1.862E-08	0.094	0.185
Kif17	28.705	8.428E-08	0.063	0.105
Prmt5	30.149	4.001E-08	0.031	0.059
Hen1	25.956	3.492E-07	0.360	0.230
Ago1	0.402	0.526	0.021	0.019
Ago2	27.078	1.954E-07	0.022	0.009
Ago3a	45.255	1.729E-11	0.025	0.006
Ago3b	9.874	0.001	0.025	0.038e
Ago4	79.917	3.905E-19	0.018	0.003
β-actin	3.359	0.066	0.013	0.018

Note.—Higher ω values are in bold.

aTwice the difference of likelihood values between two nested models (one-ratio model vs. two-ratio model).

bP values significant at the 1% threshold after Bonferroni correction in bold.

cω0 indicates a ratio of the mammalian clade (background).

dω1 indicates a different ratio of the teleost fish clade (foreground).

eAgo3b is a fish specific copy of Ago3.

Adaptive Evolution in Functional Domains of piRNA Pathway Genes

Although ω 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).

Table 2

Detection of positive selection for piRNA pathway genes and Ago genes

Gene	2ΔLa	P valueb	Estimates of the Parameters in the Modified Model Ac	Positively Selected Sitesd
Piwil1	12.997	0.0003	ω0 = 0.068; ω2 = 40.906	73K, 158A, 259S, 285R, 297C, 357V, 363V, 448S, 517G, 531Y, 535Q, 555V, 684L, 692S, 850N, 852D
			P0 = 0.861; P1 = 0.052
			P2a = 0.080; P2b = 0.004
Piwil2	17.351	3.106E-05	ω0 = 0.083; ω2 = 10.057	305P, 340T, 345V, 386P, 553S, 556G, 602T, 647C, 666V, 672S, 692E, 713I, 735S, 737P, 804T, 827S, 839S, 983Y
			P0 = 0.826; P1 = 0.076
			P2a = 0.089; P2b = 0.008
Tdrd1	28.065	1.173E-07	ω0 = 0.183; ω2 = 86.504	5F, 14L, 31P, 91N, 161Y, 163R, 352Q, 438N, 453R, 474K, 496E, 518P, 662E, 776S, 787E, 790G, 791Q, 793V, 922P, 1000P, 1011F, 1054L
			P0 = 0.588; P1 = 0.199
			P2a = 0.157; P2b = 0.053
Tdrkh	20.404	6.269E-06	ω0 = 0.116; ω2 = 60.256	69H, 123Q, 126V, 131Q, 138L, 142C, 152E, 157V, 159D, 177K, 190S, 198E, 213Q, 241Q, 251P, 348S, 373E, 385Y, 387D, 402S, 404L, 418L, 453T, 480S, 505K, 515S, 542G
			P0 = 0.598; P1 = 0.146
			P2a = 0.204; P2b = 0.050
Pld6	4.124	0.042	ω0 = 0.065; ω2 = 6.764	59H, 68P, 82S, 94S, 120V, 146S, 166G, 178L, 179T, 196L
			P0 = 0.784; P1 = 0.044
			P2a = 0.161; P2b = 0.009
Ddx4	9.197	0.002	ω0 = 0.068; ω2 = 335.073	51S, 104F, 177R, 268S, 337Q, 357A
			P0 = 0.766; P1 = 0.177
			P2a = 0.045; P2b = 0.010
Asz1	10.732	0.001	ω0 = 0.155; ω2 = 40.485	154C, 190I, 396K, 414Q, 419C
			P0 = 0.700; P1 = 0.146
			P2a = 0.126; P2b = 0.026
Mov10l1	17.078	3.587E-05	ω0 = 0.085; ω2 = 11.126	302K, 311S, 341P, 366H, 390P, 409V, 438A, 560I, 561G, 737Q, 762V, 775H, 778F, 779H, 787S, 866N, 875G, 886D, 892K, 908C, 912P, 913R, 920C, 921Q
			P0 = 0.694; P1 = 0.130
			P2a = 0.146; P2b = 0.027
Fkbp6	6.917	0.009	ω0 = 0.109; ω2 = 46.815	79T, 96A, 118Q, 158S, 164T
			P0 = 0.721; P1 = 0.049
			P2a = 0.213; P2b = 0.014
Kif17	8.569	0.003	ω0 = 0.045; ω2 = 4.658	122F, 222S, 226A, 366G, 560L, 575S, 581S, 595H, 640G, 641Q, 645G, 651V, 678T, 684V, 730E
			P0 = 0.798; P1 = 0.099
			P2a = 0.090; P2b = 0.011
Prmt5	19.085	1.250E-05	ω0 = 0.036; ω2 = 97.992	5S, 51E, 95E, 110C, 124G, 125P, 135L, 187S, 199C, 212T, 216K, 269S, 407D, 463S, 488K, 492C, 523D, 530Q, 531C, 538C, 552T, 562K, 587S, 590D, 602G
			P0 = 0.874; P1 = 0.030
			P2a = 0.091; P2b = 0.003
Hen1	4.658	0.031	ω0 = 0.160; ω2 = 99.653	179C
			P0 = 0.602; P1 = 0.351
			P2a = 0.029; P2b = 0.016
Ago1	6.078	0.013	ω0 = 0.017; ω2 = 9.523	131R, 716S
			P0 = 0.979; P1 = 0.008
			P2a = 0.011; P2b = 0.0001
Ago2	2E-06	0.999	ω0 = 0.012; ω2 = 1
			P0 = 0.964; P1 = 0.004
			P2a = 0.030; P2b = 0.0001
Ago3a	0	1	ω0 = 0.014; ω2 = 6.090
			P0 = 0.999; P1 = 0.0003
			P2a = 0; P2b = 0
Ago3b	0	1	ω0 = 0.030; ω2 = 12.715
			P0 = 0.992; P1 = 0.007
			P2a = 0; P2b = 0
Ago4	0	1	ω0 = 0.007; ω2 = 10.684
			P0 = 1; P1 = 0
			P2a = 0; P2b = 0
β-actin	4.324	0.037	ω0 = 0.012; ω2 = 998.998
			P0 = 0.982; P1 = 0.009
			P2a = 0.008; P2b = 0.00007

aThe 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 ω2 fixed at 1 is the null model.

bP values significant at the 5% threshold after Bonferroni correction in bold.

cP0 is the proportion of codons that have ω0 in all branches, P1 is the proportion of codons that have ω1 = 1 in all branches, P2a is the proportion of codons that have ω0 in the background branches but ω2 in the foreground branches, and P2b is the proportion of codons that have ω1 in the background branches but ω2 in the foreground branches.

dSites 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′-/5′-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 ω 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.

F

Distribution of the positively selected sites in functional domains. (A) Putatively positively selected sites with BEB PPs higher than 90% were shown. Functional domains in x axis: PAZ–MID (Piwi/Ago/Zwille-Middle), PIWI (P-element induced wimpy testis), Tudor, KH (K homology), PLDc (Phospholipase D family), DEAD (Asp-Glu-Ala-Asp box), Ank (Ankyrin repeat), DEXD (DEAD-like helicase superfamily ATP binding domain), AAA12 (ATPases Associated with a wide variety of cellular Activities superfamily), FKBPc (FKBP-type peptidyl–prolyl cis–trans isomerase), Kinesin (Kinesin motor domain), BBP1c (C-terminal domain of BBP1), Mtase (Arginine-N-methyltransferase). The numbers above stacked columns show the ratios of the sites distributed on functional domains/total positively selected sites for each gene. (B) Sliding window analysis shows ω value distribution (dN/dS) along the coding sequence of Piwil1 between teleosts (red line) and mammals (blue line). Functional domains of Piwil1: PAZ, MID, and PIWI. sDMA represents the binding sites for Tudor domain of Tdrd1 or Tdrkh.

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 (Chen et al. 2009; Mathioudakis et al. 2012). In addition, 44% positively selected sites were distributed in Tdrkh-KH domain, which is an RNA-binding 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′-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 protein–protein/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.

F

Structural modeling suggests functional importance of positively selected sites. NJ trees recapitulate vertebrate phylogeny. The amino acid alignments of nine representative species are shown in the middle panels. Homology modeling of zebrafish Piwil1 PAZ domain (from 275E to 387T) (A), zebrafish Piwil2 PAZ domain (from 459R to 574T) (B), the third Tudor domain (from 616L to 807G) of zebrafish Tdrd1 (C), and zebrafish Tdrkh Tudor domain (from 348S to 440C) (D). The binding pockets (piRNA 3′-end binds to PAZ; sDMA of Piwi binds to Tudor) are shown with dashed ellipse. Arrows (in amino acid alignments) or red spheres (in 3D models) indicate the positively selected sites. The modeling domains are colored by structural motifs: α helix, blue; β sheet, cyan; coil, purple. The modeling structures in Piwil1, Piwil2, and Tdrkh are not a whole domain because of incompleteness of corresponding template.

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, 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.

F

Branch model analysis of the ω values for 12 genes in the teleost lineage. (A) Cloning and expression levels of five Ago genes and six piRNA pathway genes in the swamp eel detected by RT-PCR. Hprt was used as an inner control. Ago genes are expressed ubiquitously, whereas the piRNA pathway genes are mainly expressed in three kinds of gonads except Tdrkh. (B) The branch model analysis (one-ratio model against two-ratio model). Black columns indicate ω values of the genes in the swamp eel (foreground branch), and gray columns indicate ω values of the other fish species (background branch). *LRT, P < 0.05; **P < 0.01. Ago genes and β-actin were used as controls.

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.

F

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 superfamilies in human and mouse, which are inactive in mammals. The red dots indicate the branches with rapid evolution. (B) Diverse TE loci of zebrafish piRNAs. Nonredundant piRNA sequences of zebrafish from Ketting’s study (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.

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 host–parasite 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 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 inthe 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.

Supplementary Material

Supplementary tables S1–S4 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).