We assembled a list of sixty coding genes involved in TE control (Fig. 1A, Additional file 1: Table S1). KZFPs were not included because i) they were previously analyzed and ii) their fast evolution in terms of gene copies makes it difficult to reconstruct orthology relationships [22, 28, 31]. We first sought to determine whether the protein products on these genes are engaged in a genetic conflict with TEs. If this were the case, evidence of positive selection might be expected as a signature of fast evolution in response to TE challenge. We thus focused on primates and retrieved sequence information of the coding sequences from public databases, with the aim to analyze evolutionary patterns (Additional file 1: Table S1). To test the hypothesis that positive selection has been driving the evolution of TE control genes, we applied likelihood ratio tests (LRTs) implemented in the PAML suite [32, 33]. All genes were screened for recombination and split into different regions if necessary. The neutral models (M7 and M8a) were rejected in favor of the positive selection model (M8) for 14 genes, corresponding to a considerably high fraction of 23.3% (Fig. 1A and 1B, Additional file 1: Table S2). The positively selected genes included three Piwi proteins and other molecules that participate to the piRNA pathway (SPOCD1 and TEX15), as well as components of the HUSH and NuRD complexes. In line with previous findings, we detected evidence of positive selection in MORC2 [30] and in the N-terminal region of DNMT3A [29] (Fig. 1B).

We were also interested in identifying the positively selected sites in these genes. To this aim, we used four methods: the BEB (Bayes empirical Bayes) analysis from M8, MEME (Mixed Effects Model of Evolution), FUBAR (Fast, Unconstrained Bayesian AppRoximation), and FEL (Fixed Effects Likelihood). To be conservative, a site was called as selected if it was detected by at least two methods. A total of 144 selected sites were identified (Additional file 1: Table S2). The average fraction of selected sites per protein resulted equal to 0.84%, with the highest proportion being observed for SPOCD1 (3.29%). In general, components of the piRNA pathway were found to display a high number of selection signals (Fig. 1B). Among these, we detected three sites that were independent targets of selection in the PAZ (PIWI-Argonaute-Zwille) domains of the two paralogous genes PIWIL3 (I310 and R311) and PIWIL4 (Q282). These sites define the same structural patch in the 3D structures of the corresponding Piwi proteins (Fig. 1C), suggesting that changes in this region are functionally relevant and adaptive. Although a number of positively selected sites was found to be located in the MID and PAZ domains, none of them is involved in the binding of guide RNAs.

Given their relevance for TE silencing and the signatures of positive selection we detected, we next focused on Piwi/Ago proteins and sought to determine their distribution in eukaryotes [48]. To this aim, we used the HMM (hidden Markov model) profile of the PIWI domain, which is the defining feature of this protein family, to search the reference proteomes of representative eukaryotic phyla. After exclusion of reference proteomes with a BUSCO (Benchmarking Universal Single-Copy Orthologs) completeness score lower than 80%, we analyzed 1,815 proteomes and identified 11,642 hits (Additional file 1: Table S4).

We found an extreme variability in the number of PIWI-domain containing proteins in eukaryotes (Fig. 5A). Complete loss was observed in 7.7% of the organisms we analyzed; most of these are fungi, apicomplexa protozoa, and green algae (Chlorophyta). The highest number of Piwi proteins was instead detected in two species of flowering plants in the genus Gossypium (G. barbadense and G. anomalum, phylum Streptophyta), with more than 100 proteins each. Other plant proteomes also showed a remarkable expansion of Piwi proteins, whereas some comprised relatively few (Fig. 5A). This is in line with the notion that the Piwi/Ago family has recurrently expanded and contracted during plant evolution [49].

In most metazoa, the number of Piwi proteins had a relatively narrow range, with some exceptions accounted for by species in the phyla Nematoda and Rotifera, plus the single proteomes of Tardigrada (Ramazzottius varieornatus, 18 Piwi hits) and Platyhelminths (Macrostomum lignano, 23 hits) (Fig. 5A). The wide variability in PIWI module content observed in nematodes is in line with the observation that some lineages encode an extended array of these proteins [50], whereas others have lost part of the piRNA pathway [51]. With respect to rotifers, these animals were recently shown to host a wide diversity of active TEs, which most likely triggered expansion of the Piwi pathway [52].

In simpler eukaryotes, a substantial expansion in the number of PIWI-containing proteins was observed in Ciliophora (Fig. 5A). Interestingly, in these organisms, Piwis have taken over novel functions and contribute to programmed DNA elimination, a process whereby somatic nuclear development occurs by removal of TEs and other repeat types [53]. Finally, the situation was very variegated for fungal proteomes (Fig. 5A, Additional file 1: Table S4). A remarkable expansion of Piwi proteins was observed in some organisms in the phylum Mucoromycota. Notably, all these are fungi in the Glomeromycotina subphylum of arbuscular mycorrhizal symbionts, which establish mutualistic interactions with land plants and display high TE proportions in their genomes. Remarkably, these asexual fungi and use TEs to generate genomic diversity while avoiding TE overproliferation through the Piwi pathway and other mechanisms [54].

Eukaryotic Piwi/Ago proteins are known to commonly display a similar architecture, which comprises the signature PIWI domain, plus PAZ, MID, and AgoN (Argonaute N-terminal) domains, as well as linkers and terminal IDRs (Fig. 1B). Prokaryotic Agos have instead a more diverse domain composition [55, 56]. Given the large number of PIWI-containing proteins we obtained with the HMM search, we set out to explore the diversity of eukaryotic PIWI-containing proteins in terms of domain architecture (Additional file 1: Table S4). We thus retrieved from UniProt all domains that are annotated in the eukaryotic PIWI-containing proteins. We detected a total of 72 different domains, with similar domain diversity in fungi, plants, and metazoans (Fig. 5B). Conversely, the few proteins in the other kingdoms showed that PIWI modules only co-occur with kinase domains in Alveolata. As expected, the PAZ domain was the most frequently detected, although it should be noted that its frequency may be underestimated in our analyses as we only relied on available annotations. The same applies to other domains (MID, AgoN, and Argonaute linkers) commonly found in Piwi/Ago proteins (Fig. 5B). A number of domains annotated in PIWI-containing proteins were found to function in the maintenance of genome integrity (e.g., BRCT, HORMA, MCM), in nucleic acid binding (e.g., BZIP, DRBM, HTH, MADF, RRM, SAMD1 and different types of zinc fingers) and in protein-protein interaction (DH, F box, PLAT, PUL, and SAM) (Fig. 5B). We also detected some remarkable domain associations. In three instances from fungi and green algae, the PIWI module was acquired independently by dicer-like proteins, suggesting that such proteins bring together the recognition and cleavage of mobile elements operated by dicer with the effector functions of the PIWI domain [57] (Fig. 5B and Fig. 5C). Another interesting example is represented by three related proteins encoded by fungi in the genus Colletotrichum, where the PIWI module appears in association with cytidine/dCMP deaminase (CDA) domains (Fig. 5B and Fig. 5C). The latter domains play different roles in eukaryotes, including defense from viruses or mobile elements [58]. It is thus possible that the CDA domains in these fungal Piwi/Ago proteins contribute a further effector function in the silencing of TEs or viruses. Similar considerations may apply to metazoan proteins that carry a PIWI module together with domains that are commonly found in antiviral proteins, such as DEAD box helicase and PARP (poly(ADP-ribose) polymerase) [59, 60] (Fig. 5B and Fig. 5C). Notably, we also found two plant proteins that carry both PIWI and reverse transcriptase (RT) domains (Fig. 5B and Fig. 5C). The latter have been co-opted to serve functions in cellular organisms and, in prokaryotes, they often have antiviral roles [61, 62]. Whether this is also the case in eukaryotes remains to be evaluated, but the association of PIWI and RT in these proteins suggests that it may be the case. Finally, we identified three histone acetyltransferases (HAT) in metazoa that recruited PIWI and PAZ modules (Fig. 5B and Fig. 5C). Because histone acetylation is usually associated with an active chromatin conformation, the PIWI domains in these proteins might function in gene regulation rather than TE repression [63].

Structural homology searches identify a bacterial close relative of eukaryotic Piwi/Ago proteins with conserved IDR location

We next sought to explore the deep evolutionary history of Piwi/Ago proteins. Growing evidence has indicated that a number of eukaryotic innate immunity modules, some of which also control TE replication, originated from prokaryotic antiphage systems [64]. These include Piwis/Agos, other components of the same pathway (e.g., PLD6), as well as MOV10/MOV10L1 [64]. Recently, using an HMM-based approach, Bastiaanssen and coworkers showed that Asgard achaea (or Asgardarchaeota) encode a wide diversity of Piwi/Ago genes and identified a protein (HrAgo1) encoded by the Lokiarchaeon ‘Candidatus Harpocratesius repetitus’ that shares a common origin with eukaryotic Piwi and Ago proteins [24]. An alternative strategy to identify distantly related proteins is to rely on structural conservation. The recent development of Foldseek allows the rapid screening of protein structure databases with high sensitivity [65]. We thus used the structures (or Alphafold models) of human Ago1-4 (hAgo1-4) and Piwi1-4 (hPiwi1-4) to search for proteins with similar folds in bacteria and archaea. Surprisingly, in all cases the best hit, with E values < 1.49 x 10–53, was an Argonaute protein encoded by a Planctomycetota bacterium (PlaAgo, A0A7C6F4G8) sampled from hot springs in the USA (strain SpSt-977) (Additional file 1: Table S5). PlaAgo shows higher sequence identity with hAgos than HrAgo1, making it the closest relative of human argonautes identified in bacteria or archaea. Sequence identity with hPiwis is instead lower (Fig. 6A). Pairwise structural alignments of PlaAgo with hAgo4 and hPiwi2 using TM-align showed very high scores and structural similarity covering all domains (AgoN, L1, PAZ, L2, MID, and PIWI) (Fig. 6B and 6C). Based on the comparison with hAgos, the catalytic tetrad of PlaAgo is formed by residues Asp728, Glu774, Asp806, and His946.

To explore in further detail the relationships among PlaAgo and other prokaryotic and eukaryotic Piwis/Agos, we generated a phylogenetic tree using the MID-PIWI domains of several Piwi/Ago proteins (Additional file 1: Table S6). These included both short (solely consisting of the PIWI and MID domains) and long (also including the PAZ and AgoN domains) prokaryotic Agos, either catalytically active (long AgoA) or not (long AgoB and short Agos). Results indicated that plaAgo is sister to eukaryotic Agos, whereas, in line with previous data, HrAgo1 is sister to eukaryotic Agos and Piwis (Fig. 6D) [24]. A phylogenetic tree generated from a FoldMason structural alignment of a subset of well characterized Piwis/Agos confirmed these findings (Fig. 6E, Additional file 1: Table S7).

Interestingly, in contrast to HrAgo1 (which only comprises a short disordered N-terminal region), PlaAgo presents IDRs at both the N- and C-termini (Fig. 6C), reminiscent of the domain organization of hAgos and hPiwis. We thus focused on the IDRs in PlaAgo, as protein disorder is known to be relatively uncommon in prokaryotes [66–68]. Indeed, out of 126 prokaryotic Agos we included in the phylogenetic tree, only 9 have IDRs (most of these internal rather than terminal). We first asked whether some sequence conservation between PlaAgo and hPiwis/hAgos was detectable within the IDRs. We thus inspected local similarity using the PRSS program, which computes the significance of a pairwise alignment by shuffling the amino acids in one of the two sequences. Remarkably, results indicated that, when PlaAgo is compared to hAgo1-4, significant sequence similarity extends through the C-terminal IDR. With respect to the N-terminal PlaAgo IDR, it shows significant similarity to the IDRs in hPiwi1 and hPiwi3 (Fig. 6F).

Sequence conservation of IDRs across huge evolutionary distances is quite unusual. We thus retrieved the N- and C-terminal IDRs of Piwi/Ago proteins of representative organisms from different realms of life (Additional file 1: Table S8). We found higher sequence conservation for the C-terminal IDR than for the N-terminal one, which also widely differed in length, with the former being shorter (Fig. 7A). Previous studies indicated that whereas sequence and size conservation is low in IDRs, specific features may be instead conserved across orthologs [37, 38, 69–71]. We thus used the SVR predictor [35] to calculate length-independent parameters. In addition to Sconf/N and ν, we computed SHD and NCPR. Results indicated that N-and C- terminal IDRs are relatively homogeneous in terms of conformational parameters and sequence features (Fig. 7B). Despite relatively high NCPR, those at the N-terminus have low entropy and are compact, a feature likely determined by the clustering of hydrophobic residues [42, 72] (Fig. 7B). The C-terminal IDRs are extended and have high Sconf/N; with two exceptions, they tend to be negatively charged with low SHD (Fig. 7B). Overall, these results indicated that, despite the huge phylogenetic distance and poor sequence/size similarity, the IDRs in Piwi/Ago proteins maintain some patterning features that in turn reflect into conserved conformational features.

We obtained primate orthologous coding sequences from the NCBI database (http://www.ncbi.nlm.gov). Primate genes carrying premature stop codons or with a low sequence coverage were excluded. A list of primate species analyzed for each gene is reported in Additional File 1: Table S1.

The RevTrans 2.0 utility was used to generate multiple sequence alignments (MSAs) using MAFFT v6.240 as an aligner [99, 100]. The resulting alignments were manually trimmed and adjusted to contain only well-aligning regions. MSAs were analyzed for the presence of recombination signals using Genetic Algorithm Recombination Detection (GARD, parameters: general discrete model, three rate classes) [101]. GARD is a genetic algorithm based on phylogenetic incongruence among alignment segments which detects the best-fit number and location of recombination breakpoints. When evidence of recombination was detected, the coding alignment was split based on the recombination breakpoints and sub-regions (of at least 50 nucleotides) were used as input for subsequent molecular evolution analyses. Recombination breakpoints were identified for MOV10L1, PIWIL3, and TNRC18 alignments (Additional File 1: Table S2).

Phylogenetic trees were reconstructed using the phyML program (version 3.1) with a General Time Reversible (GTR) model plus gamma-distributed rates and 4 substitution rate categories with a fixed proportion of invariable sites [102]. The MSAs and phyML trees were used as inputs to perform evolutionary analysis. To detect positive selection, the codon-based codeml program implemented in the PAML (Phylogenetic Analysis by Maximum Likelihood) suite was applied [32], using the F3x4 codon frequencies model (codon frequencies estimated from the nucleotide frequencies in the data at each codon site) [32, 33]. A model (M8, positive selection model) that allows a class of sites to evolve with dN/dS > 1 was compared to two models (M7 and M8a, neutral models) that do not allow dN/dS > 1. To assess statistical significance, twice the difference of the likelihood (ΔlnL) for each model (M8a vs M8 and M7 vs M8) is compared to a χ2 distribution (1 and 2 degrees of freedom for M8a vs M8 and M7 vs M8 comparisons, respectively). To be conservative and to obtain robust results, we called a gene as a target of positive selection only if it was detected by both M7 versus M8 and M8a versus M8 comparisons.

In order to identify specific sites subject to positive selection, we applied four different methods: 1) the Bayes Empirical Bayes (BEB) analysis (with a cutoff of 0.90), which calculates the posterior probability that each codon is from the positive selection site class (under M8 model) [103]; 2) Fast Unbiased Bayesian AppRoximation (FUBAR) (with a cutoff of 0.90), an approximate hierarchical Bayesian method that generates an unconstrained distribution of selection parameters to estimate the posterior probability of positive diversifying selection at each site in a given alignment [104]; 3) Mixed Effects Model of Evolution (MEME) (with a p-value cutoff < 0.1), which allows the distribution of dN/dS to vary from site to site and from branch to branch at a site [105]; 4) Fixed Effects Likelihood (FEL) (with a p-value cutoff < 0.1), a maximum-likelihood (ML) approach to infer dN/dS on a per-site basis, assuming that the selection pressure for each site is constant along the entire phylogeny [106]. Again, to be conservative and to limit false positives, only sites detected using at least two methods were considered as positive selection targets.

Episodic positive selection was detected using the adaptive branch-site REL test for episodic diversification (aBSREL) [107], using a subset of primate species for which data of TE genomic content is available [4]. aBSREL tests which branches in a phylogeny have a proportion of sites that have evolved under positive selection. In this case, for each gene alignment, we tested the tips of a primate phylogenetic tree derived from the zoonomia database [4].

The single-likelihood ancestor counting (SLAC) method was applied to identify sites under negative/purifying selection [106]. GARD, FUBAR, MEME, FEL, SLAC, and aBSREL analyses were run locally through the HyPhy suite V2.5.29 [108].

The sequences of reviewed (Swiss-Prot) canonical proteins (20,420) of the reference human proteome (UP000005640) were retrieved from Uniprot ( https://www.uniprot.org/ ). Intrinsically disordered regions (IDRs) were estimated using the Metapredict V2 tool [109, 110], that defines which residues from a protein sequence are disordered by applying a deep-learning algorithm based on a consensus score calculated from eight different predictors [109]. Metapredict V2 was run using default parameters and consecutive disordered stretches equal to or longer than 30 residues were labeled as IDRs. For primate genes, IDRs were kept only if the same orthologous region was predicted in at least 50% of the species analyzed.

The conformational entropy per residue (Sconf/N) and the Flory scaling exponent (ν) were calculated for all analyzed IDRs. Sconf/N is a measure of the landscape of different structures accessible to an IDR. ν derives from the scaling laws of polymers that describe how chain dimensions vary as a function of chain length [111]. Both parameter were estimated using a Colab notebook (https://colab.research.google.com/github/KULL-Centre/_2023_Tesei_IDRome/blob/main/IDR_SVR_predictor.ipynb) [35, 112], which uses a support vector regression model trained on simulations performed using the CALVADOS model [112, 113]. The same SVR predictor was used to derive other measures: SHD (sequence hydropathy decoration) [42], FCR (the fraction of charged residue), and NCPR (net charge per residue).

We conducted our analysis on a dataset of variant call format (VCF) files comprising genetic variation data of autosomal chromosomes from 54 distinct human populations and 828 individuals from the Human Genome Diversity Panel [114]. VCF files were further indexed and filtered for minimum quality parameters (-minQ/-minMapQ < 20) using bcftools and samtools [115, 116].

For each population, we calculated two summary statistics of genetic diversity: Tajima’s D and Fay’s H. These statistics were calculated in sliding windows, with a window size of 50 kilobases (kb) with a step size of 10 kb, using the software ANGSD [117]. We ranked the values for each summary statistics, and each population separately, then assigned the minimum rank value to each gene, considering all the windows encompassing its genomic coordinates.

We found no significant difference in lengths between candidate genes and the rest of the human genes. To test for signals of selection, we ranked the calculated summary statistics within each population. For each gene, we identified windows that fell into the extreme 5% of ranked values for Tajima’s D and Fay and Wu’s H (Additional File 1: Table S9). We then tagged the TE control genes that fell in those lower values and visualized the frequency of which those genes are found across the multiple populations. The analysis was conducted using RStudio and R version 4.2.2 and utilized ggplot2 package [118] for visualization.

The HMM profile of the PIWI domain (PF02171) was downloaded from InterPro (https://www.ebi.ac.uk/interpro/). The profile was used as the input for hmmsearch (version 2.39.0) to search the reference proteomes of representative eukaryotic phyla [119]. Default parameters for hit identification were applied (Significance E-value < 0.01 and < 0.03 for sequences and hits). BUSCO scores [120] were obtained from UniProt and only proteomes with a completeness score higher than 80% were retained.

Protein structures or models were obtained from the RCSB PDB (www.rcsb.org) or Alphafold database [121, 122]. IDs are available in figure legends or in Tables S5 and S7 (Additional File 1).

Foldseek (https://search.foldseek.com/) [65] was run using either the solved structures or Alphafold models of hAgo1-4 and hPiwi1-4 (Additional File 1: Table S5) limiting the search to Eubacteria and Archaea.

For pairwise structure alignment, TM-align [123] was run from the dedicated website (https://zhanggroup.org/TM-align/) using the model of PlaAgo (AF-A0A7C6F4G8-F1-model_v4) and the resolved structures of hAgo4 (6OON) or hPiwi2 (7YFX).

A structural alignment of PlaAgo and well characterized Piwis/Agos was obtained using FoldMason [124]. The multiple structural alignment (MSTA) file was used as the input for IQ-TREE [125]

A phylogenetic tree of prokaryote and eukaryote Piwis/Agos was generated by aligning representative sequences (including some from Asgard archaea in Bastiaanssen and coworkers) [24] with MAFFT [100] (Additional File 1: Table S6). The aligned sequences were used as an input for IQ-TREE.

PRSS was accessed through the University of Virginia website (https://fasta.bioch.virginia.edu/) [126] and was run using default parameters (Blosum50 scoring matrix, gap open and extension penalty − 8 and − 2, 200 uniform shuffles).

Linear Models were run using a phylogenetic generalized least squares (PGLS) model. In particular, for each primate species, we counted the number of genes showing evidence of episodic positive selection and we correlated these counts with the frequency of recently acquired transposable elements as defined by [4]. To account for the phylogenetic relationship among species, a phylogenetic tree retrieved from the zoonomia database [127] was implemented in the model under a Brownian motion process of evolution. PGLS models were run using the R package ape [128].

To evaluate whether positively selected sites are enriched in IDRs, we applied a binomial test, using the counts of sites falling in IDRs as number of successes, the number of total selected sites as trials, and the ratio between all IDR length divided by the length of all analyzed regions as the probability of success.

Wilcoxon rank sum and Kendall’s correlation tests were performed in the R v.4.0.5 environment.

All data are available in the manuscript and its supplementary files.