Evolutionary history of the vertebrate Piwi gene family

Sequence acquisition

We used Ensembl v.101 (Yates et al., 2020) and NCBI release 239 (Sharma et al., 2018) to collect known Piwi and Piwi-like coding sequences (CDS) from representative species of all major groups of vertebrates. Specifically, our sampling included cyclostomes (jawless fishes; lamprey) and gnathostomes (jawed vertebrates) from chondrichthyes (cartilaginous fishes), ray-finned fishes, lobe-finned fishes and tetrapods (amphibians, reptiles, birds and mammals). To help understand the ancestral state and focus on changes that occurred within vertebrates, we included Piwis from closely related deuterostomes as outgroups (Table S1, Fig. S1). The longest CDS was selected if there were multiple transcripts for a Piwi gene. Unidentified CDSs that shared similarity to Piwis were queried against the NCBI non redundant protein database via BLAST (Altschul et al., 1990) to confirm homology. If a Piwi paralog was absent or poorly annotated in Ensembl, we queried the NCBI for the paralog using human Piwi sequences as BLAST queries.

Phylogenetic reconstruction

All Piwi sequences were translated to amino acids and aligned using the LINSI strategy in MAFFT v7 (Katoh & Standley, 2013). We constructed phylogenetic trees from the amino acid alignment with two methods. First, we used IQ-Tree2 (Minh et al., 2020) and let IQ-Tree2 find the best fitting substitution model (Kalyaanamoorthy et al., 2017). Ten independent phylogenetic analyses were run simultaneously to explore tree space and we selected the tree with the highest likelihood score. Node support for the best scoring tree was evaluated with the ultrafast bootstrap method (Hoang et al., 2018) and the SH-like approximate likelihood ratio test (SH-aLRT) (Guindon et al., 2010). In both cases, we used 1,000 pseudoreplicates. Second, we used ExaBayes v1.5.1 (Aberer, Kobert & Stamatakis, 2014) to conduct Bayesian analyses. We ran four simultaneous chains for 1,000,000 generations sampling every 500 generations, using the same substitution model chosen by IQ-Tree2 but without additional model parameters. Chains were considered as having converged once the average standard deviation of split frequency among runs was less than 0.01. We summarized results with a majority-rule consensus tree from the set of sampled trees after the first 25% were discarded. We repeated these tree construction analyses on a reduced alignment generated after removing poorly aligned regions using Gblocks and allowed up to five contiguous non-conserved positions and columns to contain gaps; -b3 = 5 -b5 = a (Castresana, 2000) to test the impact of alignment quality on our results.

Synteny analyses

To understand the conservation of synteny as well as gene gain and loss, we examined the position of protein coding genes up and downstream of Piwi genes in representative vertebrates. We used gene coordinate information from Ensembl v.101 as well as estimates of orthology and paralogy from the Compara database (Herrero et al. 2016). If a Piwi gene was presumed missing, we used BLAST in Ensembl with human Piwi queries to locate any Piwi derived fragments near expected positions. If no fragments were found, we concluded the gene was missing, but if fragments were found and a full annotation was lacking, we presumed the Piwi to be pseudogenized.

Estimating d_N/d_S rates and selection tests

For each Piwi paralog in each major lineage with multiple samples (birds, mammals, and ray-finned fishes), we made independent codon alignments by translating DNA sequences to amino acid sequences, generating LINSI amino acid alignments and ‘reverse translating’ alignments into codon alignments using a custom Python script (available at github.com/mike2vandy/PiwiGeneFamily along with other scripts used). We specifically conducted site-tests for diversifying selection (Goldman & Yang, 1994) among alignments using the codeml module of PAML v4.9. We calculated the likelihood of models that allow d_N/d_S to vary among codon positions (M1a, M2a, M8a and M8). We used the likelihood ratio test (LRT) to test for significant differences between nested models that do not allow selection to those that do (M1a vs. M2a, M8a vs. M8). We also estimated a single d_N/d_S rate for each alignment using codeml’s M0 one-ratio model and pairwise d_N/d_S distances using the modified Nei-Gojobori (proportional) method in MEGA-X (Kumar et al., 2018). We noticed an intronic sequence that included stop codons was retained in the manatee Piwil4 sequence. We suspect this was an annotation artifact and because of it, this sequence was not included in selection tests and d_N/d_S estimates.

To view changes in d_N/d_S along the mammalian Piwil1 and Piwil3 branches we constructed a phylogenetic tree of mammalian Piwil1 and Piwil3 and used a free-ratio model to estimate d_N/d_S rates on all branches (Yang, 1998). We also used the adaptive branch-site random effects likelihood (aBSREL; (Smith et al., 2015)) approach implemented in Datamonkey v2.2 (Weaver et al., 2018) which detects episodes of positive selection on all designated foreground branches. We performed a LRT between the null model (d_N/d_S = 1) against the alternative where a branch was undergoing some form of selection (d_N/d_S ≠ 1). All Piwil3 branches were labeled as foreground and tested.

Expression analyses

To assess historical changes in Piwi expression through vertebrate evolution, we collected RNASeq data from representative tissues (brain, heart, kidney, liver, ovary and testes) and species (elephant shark, spotted gar, clawed frog, spiny toad, chicken, opossum, cow, and human) from NCBI’s short read archive (SRA; Leinonen, Sugawara & Shumway, 2011). Accession numbers for species and tissues can be found in Table S2. For all species, we used appropriate Ensembl’s CDSs downloaded from BioMart as a mapping reference. We only used the longest CDS per gene. Prior to mapping, RNASeq reads were cleaned with Trimmomatic v0.38 (Bolger, Lohse & Usadel, 2014) using the trimming and adapter removal parameters: HEADCROP:5, SLIDINGWINDOW:5:30, MINLEN:50, and ILLUMINACLIP:2:30:10. RNASeq reads were mapped to reference CDS libraries using default parameters of RSEM v1.3.1 (Li & Dewey, 2011) and Bowtie v1.2.2 (Langmead et al., 2009) to estimate gene expression in units of transcripts per million (TPM).

Phylogenetic analysis

Our first goal was to reconstruct the Piwi phylogeny from vertebrate representatives and deuterostome outgroups using both maximum likelihood (ML) and Bayesian strategies. We collected 168 sequences from 63 species (Table S1). We generated trees in three different strategies: a tree of select representative species (Fig. 1), a tree from all collected sequences (Fig. S2), and because the N-terminus was difficult to align, we also estimated a tree from the most alignable regions of the amino acid alignment which included the PAZ and PIWI domains (Fig. S3). All trees were relatively similar and sequences fell into the two major clades, the Piwil1 group and Piwil2 group (Fig. 1). There have been no additional duplications of Piwil2 among deuterostomes and the subfamily is monophyletic. By contrast, Piwil4 and the mammalian specific Piwil3 were derived from independent duplications of a vertebrate Piwil1 (Fig. 1). These additional duplications within the Piwil1 clade indicate that the gnathostome Piwil1 is not a 1:1 ortholog of the ancestral Piwil1 gene. However, we have retained the traditional nomenclature for simplicity.

The Piwil4 gene of vertebrates was likely derived from an ancient duplication of Piwil1 in the ancestor of gnathostomes, however a lamprey Piwil1-like sequence was found sister to gnathostome Piwil4 in both ML and Bayesian trees. Therefore, we can interpret the Piwil4 duplication in two different ways. The Piwil4 duplication occurred in the ancestor of vertebrates, but Piwil1 was lost in lamprey and the lamprey sequence is an ortholog of Piwil4. Alternatively, the Piwil4 duplication occurred in gnathostomes, but the lamprey Piwil1-like gene is artificially sister to Piwil4. Piwil4 was present in cartilaginous and lobe-finned fishes, mammals and reptiles, but was absent in ray-finned fishes and birds (Fig. 1; Fig. S1).

Our analyses suggest Piwil3 is the result of a duplication of Piwil1 in the common ancestor of therian mammals given its presence in marsupial and placental mammals and its absence from Monotreme genomes (Fig. 1). However, not all placental mammals have retained Piwil3. Piwil3 was present in marsupials (opossum and koala), Laurasiatherians (bat, dog, and cow) and most Euarchontoglires (primates, rodents and relatives), but we were unable to find it in Afrotherians (elephant and manatee), Xenarthrans (armadillo) and mouse-like rodents which included the kangaroo rat, deer mouse and house mouse (Fig. S1).

An additional copy of Piwil1 (Piwil1a and Piwil1b) had been previously described in the clawed frog (Wilczynska et al., 2009), but our analyses also revealed this additional paralog was present in the spiny toad (Fig. S2). Interestingly, a third intact Piwil1 copy was identified in the spiny toad (Fig. S2). This copy is most likely a retrogene, diagnosed by a lack of introns in the coding sequence (Fig. S4A).

Synteny analyses

We used genomic coordinates to identify genes up- and down-stream of Piwi genes to resolve or confirm homology and infer the duplication mechanism behind the vertebrate Piwi expansion. Synteny analyses included representative species with contiguous regions for at least 10 genes up and downstream of each Piwi. In most cases, this included the elephant shark, spotted gar, zebra fish, clawed frog and chicken and human, but varied among Piwis. The lamprey assembly on Ensembl lacked contiguity around the Piwi genes, therefore we queried sea lamprey coordinate data from a GFF file at https://genomes.stowers.org/sealamprey (Smith et al. 2018). Unfortunately, there was a lack of synteny among Piwi genes between cyclostomes and gnathostomes, but synteny around Piwil1, Piwil2 and Piwil4 has generally been conserved among gnathostomes (Fig. 2; Fig. S5). Synteny around Piwil2 was less conserved than Piwil1 and observed three different sets of upstream genes, however Slc39a14 and Ppp3cc were consistently downstream of Piwil2 in all lineages (Fig. S5).

Piwil1 was flanked by Rimbp2, Stx2, and Ran on the 3′end in all gnathostomes, but synteny was only conserved upstream of Piwil1 among bony and jawed vertebrates where Piwil1 is generally flanked by Fzd10, Tmem132d and Glt1d1 (Fig. 2A). By contrast, synteny around Piwil4 has been conserved throughout all gnathostome evolution. Although Piwil4 was lost in ray-finned fish and birds, the overall locus has remained intact (Fig. 3B). In the spotted gar, the distance between Fut4 and Amotl1 was 18,372 bases and there were no BLAST hits to Piwi genes in that region. The same was true for the chicken, where the distance between Fut4 and Amotl1 was 45,906 and lacked sequences similar to Piwil4. By contrast, the same region in the crocodile covers 157,135 bases and includes a Piwil4 gene (Fig. 3B). Interestingly, Piwil4 is intact in all mammals except the platypus. The Piwil4 locus in platypus does contain a lncRNA annotation (ENSOANT00000069592) between Fut4 and Amotl1 with similarity to Piwil4. However, a complete gene with an open reading frame could not be recovered. Therefore we suspect that Piwil4 has been pseudogenized in the platypus (Fig. 2B) but has not been fully purged from the genome.

Piwil1 and Piwil4 were not flanked by any co-duplicating gene families (Fig. 2), an indicative signal of whole genome or segmental duplications (Campanini et al., 2015; Opazo et al., 2015). Furthermore, both teleosts and salmonid fishes, which have experienced additional whole genome duplication events (Jaillon et al., 2004; Macqueen & Johnston, 2014; Lien et al., 2016), lacked any additional Piwi paralogs (Fig. S2).

Piwil3 is a direct neighbor of Piwil1 in marsupials which would suggest that the Piwil1 and Piwil3 genes of therian mammals are co-orthologs of the ancestral mammalian Piwil1 gene derived from a tandem duplication of Piwil1 in their common ancestor (Fig. 3A). Interestingly, the position of Piwil3 is conserved among marsupials, but Piwil3 is in a novel location in placental mammals where flanking genes vary among lineages (Fig. 3A). Taking these results together, the expansion of the Piwi gene family in vertebrates can best be explained by tandem duplication followed by a translocation due to lineage-specific genomic rearrangements.

Piwil3 and selection tests

We observed very long branches among Piwil3 paralogs (Fig. 1) and decided to test the role diversifying selection has played on the evolution of Piwil3. Using PAML’s free ratio model, we found that d_N/d_S rates among Piwil3 sequences were considerably faster than rates among therian Piwil1 sequences (Fig. 3B). Since the free-ratio model does not offer any additional information besides changes in d_N/d_S and is not statistically robust, we used site-tests among Piwil3 sequences and found that models allowing selection (M2a and M8) were significantly better fits to the data than models that disallow positive selection (M1a and M8a) (Table S3). Codons throughout the entire gene, including the conserved PAZ and PIWI domains (Fig. S6), exhibited signatures of positive selection (Fig. 3C; Fig. S7) according to Bayes Empirical Bayes (BEB) estimation of site posterior probabilities. A separate branch-site test, aBSREL in HyPhy, suggested 17 out of 22 tested Piwil3 branches have undergone episodic positive selection (Fig. S8).

We then compared Piwil3 patterns of evolution to the other paralogs within tetrapod groups (birds, mammals and ray-finned fishes). Site-tests of diversifying selection suggested no other Piwi paralog is evolving in a diversifying pattern. In all cases except Piwil3, models allowing positive selection were not improved over models restricting selection (Table S3). Further, one-ratio (Fig. 3D) and pairwise (Fig. 3E) d_N/d_S, although likely inflated due to saturation at synonymous sites, indicated Piwil3 is the fastest evolving Piwi among all lineages while Piwil1 in mammals is highly conserved.

Historical Piwi expression

Our last goal was to identify any differences in Piwi expression among lineages, and consistent with previous results, Piwi expression was largely localized to the gonads, specifically testis, although Piwis were sporadically expressed at low levels among somatic tissues (Fig. 4A). Aravin et al. (2008) reported that Piwil4 is only expressed for a short period of time during testis development in mice, but Piwil4 has since been detected in adult undifferentiated spermatogonia (Carrieri et al., 2017; Vasiliauskait et al., 2018). Consistently, Piwil4 was detected in the adult testes of all species examined, but lowly detected in the mouse (Fig. 4B). In addition, we detected weak expression of Piwil3 in the testis of the opossum, human, and cattle. Both Piwil1 paralogs were expressed in the testis of the clawed frog (Fig. 4B). We were also able to measure Piwi expression from the brain, ovary, and testes, in the spiny toad which had an additional Piwil1 retrogene. The spiny toad expression pattern was unique as Piwil4 and Piwil1a-retro were the only Piwis expressed in the testis, but Piwil1a and Piwil1a-retro were expressed in the ovary (Fig. S4B). Piwil1b nor Piwil2 were expressed in any library (Fig. S4B).