Small RNA library construction and sequencing

We initially constructed small RNA libraries using size selection (gel excision) of total RNA extracted from fresh A. vaga biomass. One of the samples was irradiated to check whether the small RNA populations would be affected under conditions causing significant DNA damage. Sequencing on the Illumina GA IIx platform, after filtering to remove adaptors and low-quality sequences, yielded ∼9.6 M and ∼10.3 M sequences (a total of ∼20-M reads) for the wild-type (WT) and IR samples, respectively.

We observed a high level of A. vaga ribosomal RNA representation in our libraries, especially in the IR sample [15% ribosomal RNA (rRNA) in WT and 62% rRNA in IR], which is likely due to overabundance of rRNA breakdown products after irradiation, as was previously observed by others (Lee et al. 2009). In addition, size-selected libraries also contained transfer RNA (tRNA) breakdown products. After filtering of rRNA and tRNA sequences, the combined dataset was reduced to ∼6.6-M reads, of which ∼3 M could be mapped to the reference genome.

The patterns of length distribution in small RNA populations were found to be different in comparisons between WT and IR samples, especially in the 22- to 24-nt fraction, which was significantly increased in the IR sample, as was previously observed in other irradiated species, where siRNAs were shown to be involved in DNA repair (Wei et al. 2012). Of particular interest, however, was the pronounced shoulder formed by longer RNA species, which we inferred to correspond to piRNAs (Figure S1A). Typically, piRNAs are 25–30 nt in length and display a strong 5′-uridine bias (Vagin et al. 2006; Aravin et al. 2007a,b; Brennecke et al. 2007; Gunawardane et al. 2007; Siomi et al. 2008, 2011; Iwasaki et al. 2015). In agreement with expectations, the initial annotations met the correlation between TEs and potential piRNAs, in accordance with the role of piRNAs as central players in transposon silencing. Indeed, the 25- to 31-nt shoulder (∼1.3-M reads) exhibited a strong 5′-end uridine bias, especially in the 29-nt area of the IR sample (Figure S1B). However, the overwhelming majority of reads was represented by a large peak composed of siRNAs and miRNAs (Figure S1A), and it was difficult to achieve full separation of the presumed piRNA fraction in the shoulder to perform meaningful analyses. In addition, TE-mapped reads displayed highly nonuniform distribution across TE length (Figure S1D), which was likely due to low representation of piRNAs in the samples dominated by rRNAs and other types of small RNAs and the resulting amplification bias.

While our initial experiments agreed with the involvement of small RNAs in repair of DNA damage (Francia et al. 2012; Wei et al. 2012), we were particularly interested in examining the population of piRNAs, which suppresses TE activity in the germ line and may act to establish trans-generational silencing. To investigate piRNAs in A. vaga, we had to select a different method that would result in better representation of the piRNA population, since the commercially available ribo-depletion kits proved ineffective in this species. In model organisms, piRNAs can be purified through their association with Piwi proteins following immunoprecipitation with the corresponding antibodies. However, the sheer number of Piwi variants in A. vaga (eight different proteins) does not make it feasible to employ anti-Piwi antibodies for piRNA purification. We therefore decided to use Hitrap-Q chromatography columns for isolation of protein-bound small RNAs and simultaneous elimination of rRNA/tRNA, which results in significant enrichment with Argonaute complex proteins and small RNAs bound to them (Lau et al. 2006). The Illumina HiSeq library construction and sequencing of protein-bound sRNA molecules resulted in 21.7 million reads (two wild-type replicas). By using a mapping procedure to report only the alignments in the best alignment stratum (i.e., those having the least number of mismatches) and choosing the unique mapping option (one read — one mapped location), we matched a total of ∼19.3-M sequencing reads (91.13%) to the A. vaga reference genome. Since this method did not yield any differences between WT and IR samples, presumably due to near saturation of protein complexes by preexisting small RNAs, we proceeded with the analysis of WT samples to obtain a comprehensive picture of piRNA distribution in the A. vaga genome.

Transposable elements and small RNAs

The length distribution plot of the reads mapped to the A. vaga reference genome shows that they are represented mostly by 25- to 32-nt-long sRNAs, peaking at 29 nt (Figure 1A). Sequences longer than 24 nt exhibit a strong 5′-U bias (e.g., 96% of 29-nt reads have uridine at the 5′-end; Figure 1B). After filtering reads not assigned to any annotated features such as TEs or gene models, small RNAs showed a strong preference toward annotated transposons, and the majority of reads were mapped to TEs in reverse (antisense) orientation: while 71% of all reads were derived from TEs (sense and antisense strands combined), 66% of all reads were in the antisense direction. Since these pi-like RNAs meet the typical piRNA characteristics (25- to 32-nt size range, peaking at 29 nt; strong 5′-U bias; most sequences are TE derived), we further consider these sequences as piRNAs, although technically they were not isolated using anti-Piwi antibodies. However, to keep in mind that a small number of the studied sRNA population may represent noncoding small RNAs other than piRNAs, we collectively designate them as sRNA, even though the vast majority fits the piRNA category.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Characterization of protein-bound small RNAs in A. vaga. (A) Length distribution of sRNAs mapped to annotated genomic features (TEs, CDS, UTRs, and introns). The y-axis shows combined read counts for both replicates. (B) Percentage of 5′-terminal nucleotides (A, T, G, and C) in each sRNA size class. (C) The share of the assembly (in kilobases) occupied by major TE groups (top) and distribution of sRNA reads between TE groups (bottom) in the sense (blue) or antisense (red) orientation. Crypton TEs are not shown due to insignificant contribution to both plots. (D) Distribution of small RNAs across the AI gene categories: annotated CDS (blue), UTR (orange), and intron (yellow) features in the antisense (left) and sense (right) orientation. (E) The 5′-terminal nucleotide composition of small RNAs mapped to A. vaga genes of putatively foreign origin (AI ≥ 0); nucleotides are colored as in B.

It was of interest to investigate the piRNA population for signatures of possible amplification mechanisms. One of the models of piRNA biogenesis, which applies to Drosophila and many other species, proposes that secondary piRNAs are generated through iterative Piwi-mediated cleavage of transcripts with complementary sequence, known as the ping-pong amplification cycle (Brennecke et al. 2007; Gunawardane et al. 2007). We investigated the overlap (in nucleotides) between stranded and reverse reads corresponding to TE annotations. According to the ping-pong model, piRNAs that map to opposite genomic strands tend to have a constant overlap (e.g., 10 nt), in which Argonaute proteins bound to sense-strand piRNAs catalyze antisense-strand cleavage at the tenth-nucleotide base pair that generates the 5′-end of antisense piRNAs. However, in a plot of overlap values for TEs, there is no obvious ping-pong amplification signature at any overlap size (Figure S2A). This result, together with the lack of internal nucleotide bias in sense reads (Figure S2C), leaves open the question of how secondary piRNAs are generated and amplified in A. vaga, with numerous RdRPs (20 copies; Flot et al. 2013) presumably playing a prominent role.

The initial TE annotation in A. vaga covered ∼3% of the genome (Flot et al. 2013). Even after TE reannotation, which included previously unknown TE types (I. R. Arkhipova, F. Rodriguez, and I. A. Yushenova, unpublished results) and used the highly sensitive RepeatMasker tool yielding a number of false positives, TE annotations still cover only 4.4% of the genome, which is much less than TE content reported for other metazoans (Arkhipova and Rodriguez 2013). The top panel in Figure 1C depicts the share of the genome (in kilobases) occupied by each of the major TE subclasses. We also reported the high diversity of TE families, which are represented by only one or two full-length copies (209 of 255 families) and about 10 times as many fragmented copies (Flot et al. 2013). Figure S3A illustrates that the low ratio of full-length copies to partial fragments is more or less maintained in each TE subclass.

Small RNAs cover 70.2% of cumulative annotated TE length (not considering different coverage depths or strandedness), thus showing a significant improvement over the size-selected library prepared from total RNA (4.9%). For the vast majority of active TE copies, sRNA coverage is sufficient to target them for silencing, since a TE does not have to be fully covered by sRNAs over its entire length to be subjected to silencing. Those copies that do not yield significant sRNA coverage (<10 counts) usually belong to transcriptionally inactive and decaying families. The distribution of sRNA reads is rather uniform over the entire TE length, not displaying any bias toward specific TE regions such as coding or repeated sequences. While sRNA mapping in most cases agrees with preexisting TE annotations, in some cases we were able to extend host–TE boundaries for previously unknown TE types by relying on small RNA coverage (El Baidouri et al. 2015; I. R. Arkhipova, F. Rodriguez, and I. A. Yushenova, unpublished results). Nevertheless, piRNAs show certain differences in numbers per class or superfamily. The bottom part of Figure 1C shows the distribution of small RNAs, in sense or reverse orientation, across the major A. vaga TE subclasses. The terminal inverted repeat (TIR)- containing and Athena elements yield a substantial share of reverse sRNAs (1.63 M and 1.18 M reads, respectively), which is expected because both subclasses occupy a bigger proportion of the genome than the others, as may be seen in the upper panel. However, LTR and non-LTR retrotransposons, which show relatively moderate genome coverage when compared to TIR or Athena, also comprise a substantial share of reverse reads (0.88 M and 1.39 M, respectively), which is indicative of higher density of sRNA coverage in these retrotransposons or perhaps more frequent collapse of several copies into a single contig.

To assess the transcriptional activity of TEs, we used transcriptome profiling. To this end, 76-bp Illumina reads from cDNA libraries were mapped to the A. vaga genome as described in Flot et al. (2013) and used to compare distributions of aligned RNA-seq reads and sRNA reads in annotated TE families. Figure S4 shows RNA-seq and sRNA count plots for previously described individual A. vaga retrotransposon (LTR, non-LTR, Penelope) and DNA TE (TIR and Helitron) families (Arkhipova and Meselson 2005b; Gladyshev and Arkhipova 2010; Arkhipova et al. 2013; Flot et al. 2013). In most cases, TE families with higher antisense sRNA counts also display lower RNA-seq counts, with a few exceptions. For instance, the high transcript levels in Zator1 are confined to the N-terminal part of the single ORF containing an extended coiled-coil motif and yield a truncated inactive transposase, which may potentially act as a repressor (Majumdar and Rio 2015). In other cases where TEs represent recent arrivals and may be undergoing an initial increase, such as Hebe or Tc/mariners (Arkhipova and Meselson 2005b; Gladyshev and Arkhipova 2010), the distribution of RNA-seq counts is biased toward several copies, which may have been collapsed together. A curious exception is the rotifer-specific Soliton clade of non-LTR retrotransposons: several Soliton families display atypically high RNA-seq/sRNA count ratios (Figure S4B), but are nevertheless characterized by particularly low copy numbers and by virtual absence of fragmented copies. Apparently, for these TEs, higher transcriptional activity is not correlated with higher proliferative capacity, perhaps being limited to a certain somatic tissue.

Small RNAs and genes of foreign origin

In the sRNA population that is not derived from TEs, a significant fraction (30%) originates from predicted CDS (Figure 1A), and 26% of all non-TE hits matched genes with AI > 0 (Figure S3B), many of which may have been acquired by HGT. Characteristically, sRNA coverage of such genes generally exhibits patterns that are similar to those observed in TEs. Specifically, 81% of sRNAs mapping to putatively foreign genes are represented by 25- to 32-nt reads, with 29-nt reads being the most abundant (either in reverse or sense orientation), and >79% of reads in the 25- to 32-nt size range are mapped to annotated CDS regions.

To characterize small RNAs matching foreign genes, we associated each GeneID with the AI value, distinguishing between reads corresponding to sense (stranded) direction or antisense (reverse) direction in each gene feature (UTRs, introns, and CDS). Small RNAs from the reverse strand of transcripts comprise the majority of reads (77%) in AI ≥ 0 genes, for which CDS features are precisely targeted (67%) by antisense reads. For sense reads in the 25- to 32-nt size range, 39% targeted CDS features, while a noticeable input was provided by intron regions (41%).

We also observed that small RNAs associated with foreign genes in A. vaga have a strong bias for 5′-U (Figure 1E), which is the terminal nucleotide in 87% of reads in the 25- to 32-nt size range. Again, as in pi-like small RNAs mapping to TEs, there is no obvious 5′-nucleotide bias in sense reads (data not shown). Nor could we observe a ping-pong amplification signature with a constant overlap between reverse and sense reads (Figure S2B).

To analyze sRNAs associated with the gene dataset, we divided it into four categories: AI ≥ 30; 0 ≤ AI < 30; −30 < AI < 0; and AI ≤ −30 (Figure S3, B and C; Materials and Methods). All sRNAs mapped to annotated genes were subdivided into UTR-, intron-, or CDS-matching reads. At the same time, a distinction was set based on read orientation (sense or antisense). Notably, while the bona fide foreign genes (AI ≥ 30) account for ∼10% of genes with BLAST hits (Figure S3B), they constitute nearly a quarter of all genes with sRNA coverage. The proportion of genes with piRNAs mapping to CDS for each AI category is shown in Figure S3C, where the number of genes with piRNA coverage (from two biological replicates with a minimum count of 10 per gene) is plotted for each AI group, showing that the proportion of CDS-matching antisense reads in AI ≥ 30 genes does not depend on a few genes with high sRNA count.

Notably, Figure 1D shows that across the AI categories, i.e., from nonmetazoan to host genes, a transition is clearly observed in antisense reads, featuring a shift from CDS (over 70% in foreign AI ≥ 30) to UTR (∼60% in metazoan AI ≤ −30). Thus, while our results for metazoan genes agree with previous observations that genic piRNAs are mostly 3′-UTR directed (Robine et al. 2009; Saito et al. 2009; Siomi et al. 2011; Le Thomas et al. 2014), the putatively foreign genes largely display piRNA coverage patterns that are similar to TEs in uniformly targeting the gene body, rather than the 3′ UTR.

Overall, we find that the patterns of piRNAs matched to numerous alien gene families resemble the patterns derived from TEs in read distribution and orientation, indicating that the RNA-mediated silencing machinery may be acting to regulate expression of putatively foreign genes.

Small RNA profiles in foreign multigene families

To further understand the interactions between sRNA machinery and foreign genes and to place them into phylogenetic context, we sought to investigate the immediate genomic environments and phylogenetic relationships in sufficiently expanded gene families and their nonmetazoan counterparts and to compare RNA profiles (RNA-seq and sRNA) in each group. Analysis of multigene family members in the phylogenetic context is particularly informative, since it allows one to discriminate between relatively recent insertion/duplication events and more ancient gene acquisitions residing in the host genome for a long time and to determine whether the more recent copies are more likely to be subjected to piRNA response. The presence of repetitive and/or modular structures in many expanded gene families complicates their proper taxonomic assignment and phylogenetic analysis, and therefore such families could not be considered in the phylogenetic context. Here, we focus on two foreign multigene families devoid of repetitive and modular structures: β-lactamases (bla) and NAD:arginine ADP-ribosyltransferases (ART).

Serine β-lactamases are known to provide resistance to β-lactam antibiotics such as penicillin or cephalosporin (Ambler 1980). They are found in the majority of principal bacterial groups or taxonomic subdivisions. We extracted all A. vaga CDS homologous to serine β-lactamases (36 copies, with AI between 22 and 304; File S1, sheet 2) and subjected them to phylogenetic analysis (Figure 2A), including representatives of known bacterial β-lactamase classes (Hall and Barlow 2004). Each node giving rise to an A. vaga β-lactamase clade that also includes a bacterial counterpart with significant clade support value apparently corresponds to an independent horizontal transfer, yielding at least nine such events (and possibly more; see below).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

Small RNAs in the A. vaga β-lactamase family. (A) Phylogenetic relationships of A. vaga serine β-lactamase homologs and their bacterial counterparts. The codon-based maximum likelihood phylogenetic tree of serine β-lactamases includes bacterial classes A, C, and D (in black) and the A. vaga homologs (in red). Red stars denote A. vaga genes with small RNA coverage. The node with an asterisk (*) is expanded in B to show RNA profiles. The key for bacterial subdivisions is shown in the top left corner. (B) Dendrogram and RNA profile histogram of A. vaga β-lactamase homologs. The top histogram shows antisense sRNA counts and RNA-seq counts for the clade in A marked with *. The dendrogram below shows phylogenetic relationships of bla genes and nucleotide sequence identity (%) between nearest neighbors, which can be grouped into allelic pairs (yellow boxes) by synteny analysis. The ab initio gene prediction model SNAPT00018705001 was included as a member of the allelic pair. (C) Genomic environments of A. vaga bla genes with sRNA coverage. The alignment is centered on the seven leftmost bla genes from B, organized in two high-identity groups, with the unpaired scaffold 32 in the middle. Predicted genes, yellow block arrows; TE insertions, light blue bars; and telomeric repeats, small red arrowheads. Regions of homology are gray shaded, and the intensity of shading corresponds to nucleotide sequence identity (%) as indicated. BLASTN output was used for pairwise comparison of scaffolds, plotting similarities with e-values <0.001. The figure was produced using Easyfig 2.2.2 (Sullivan et al. 2011).

Not every copy of β-lactamase (bla) in A. vaga exhibits the same RNA profile, either in sRNA or RNA-seq profiles. Figure 2B provides a dendrogram representation of the clade bracketed in the β-lactamase tree of Figure 2A. Nucleotide sequence identity is shown for nearest neighbors, which can be grouped into allelic pairs after checking for collinearity of adjacent genes. The sRNA profiles are clearly observed in two groups in the left part of the chart, represented by seven gene models that form the most recent branches in the phylogenetic tree. Each of the two groups contains allelic pairs with high nucleotide sequence identity (>98%) located in collinear gene environments and an additional duplicated copy. The additional annotated bla copies in each group may represent assembly variants due to polymorphic TE insertions; however, they still draw sRNA and RNA-seq counts regardless of their activity status. GSADVT00009239001 does not share the local environment with three other members of the clade and does not have an allelic partner, but its coverage is approximately double in comparison with other members of this clade, and could therefore represent an allelic pair formed by gene conversion (Flot et al. 2013). The sRNA and RNA-seq profiles from each member of a pair display similar counts, indicating that identical reads are randomly assigned to either of the high-identity copies.

Importantly, all bla clades with sRNA coverage are likely to represent relatively recent insertions into subtelomeric regions, as may be judged from their phylogenetic placement and the presence of short stretches of telomeric repeats surrounding these bla copies (Figure 2C). The immediate genomic environment of bla genes in each group is different, indicating independent origin of each segment framed by telomeric repeats. Telomeric repeat signatures are not seen in the more basal branch, perhaps because of gradual erosion over time (Figure 2C, middle scaffold 32), and all of the earlier-branching clades, arranged in allelic pairs with syntenic environments (Figure 2B, right half), lack interstitial telomeric repeat stretches and do not display sRNA coverage, which may indicate full incorporation of these genes into the host environment.

Another expanded foreign multigene family of relatively simple structure is ART. While its bacterial counterparts are primarily annotated as uracil-DNA glycosylases (UDG), most of these genes in A. vaga lack the UDG moiety, and contain only the C-terminal ADP-ribosylating domain of certain bacterial UDG-like enzymes, which is most similar to the VIP2 family of protist and bacterial actin-ADP-ribosylating insecticidal binary toxins. Of 156 A. vaga genes of probable nonmetazoan origin (AI > 0) with ART domain, 69 copies formed an ortholog group. As in β-lactamases, a variety of RNA profiles is observed in A. vaga ART homologs (Figure 3B). Two of these show a particularly high sRNA peak. After inspection of their genomic context, they can be considered as an allelic pair with collinearity in the region of scaffold overlap (Figure 3A) and, being 95% identical, share most of the RNA reads. Both have an oppositely oriented insertion of the same TIR fragment (Avmar1a) at their 3′ flanks, which may have played a role in establishing piRNA production from this gene. In addition, GSADVT00063672001 has a unique P-element insertion at the 5′-flank, which, however, shows little piRNA coverage. The decaying pseudogene retained in one of the haplotypes has the footprints of microhomology-mediated deletion and is likely on the way out. Again, the presence of interstitial telomeric repeats framing the ART-containing segment (Figure 3A, red triangles) supports its relocation to the subtelomeric region.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

Small RNAs in the A. vaga ART family. (A) Synteny analysis and sRNA plots in two ART-containing scaffolds. Annotation and genomic context of gene models GSADVT00063672001 on scaffold 680 (1–13,600) and GSADVT00013577001 on scaffold 54 (13,000–1) are shown. Features are designated as in Figure 2C; Ψ, pseudogene. Plots above each region show sRNA read counts (blue, forward strand; red, reverse strand). It is not known whether synteny continues beyond the 5′-end of the scaffolds, as both start with the region of interest. (B) The top histogram shows RNA profiles for A. vaga ART genes (AI > 0), which clustered together as an ortholog group; the dendrogram below shows their phylogenetic relationships; labels are as in Figure 2B. Transcripts with high sRNA counts, assigned to an allelic pair after inspection of their genomic environment, are compared in A.

Finally, an example of a foreign gene family with very few members also illustrates the greater likelihood of piRNA response in genes recently relocated to subtelomeres. The A. vaga RVT genes are domesticated single-copy reverse transcriptases of fungal/microsporidian origin (Gladyshev and Arkhipova 2011). Figure 4A shows that the two RVT_A copies lack introns, are embedded in a genomic environment rich in TEs and other foreign genes, are surrounded by telomeric repeats, and elicit a pronounced piRNA response, which is possibly driven by an LTR remnant at the 3′-end present in both alleles. In contrast, the two alleles in the transcriptionally active RVT_B lineage are located in a metazoan gene-rich environment, lack telomeric repeats in the vicinity, have acquired introns, and are virtually devoid of piRNA coverage (Figure 4B).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

Small RNAs in the A. vaga RVT family. The sRNA coverage plots are shown for allelic scaffolds s546(18,500–50,000)/s1012(26,700–1) (A) and s7(1–14,500)/s24(20,700–35,300) (B), comprising the rvtA and rvtB gene pairs, respectively. Features are designated as in Figure 2C. Scaffold 1012 has an assembly gap (mostly filled with scaffold 6459) within the rvt gene, which was previously sequenced in full for both members of the allelic pair from a fosmid library (JN235987.1 and JN235988.1; Gladyshev and Arkhipova 2011). An additional pseudogenized pair (rvtC) does not yield transcripts or sRNAs and is not shown.

Host genes targeted by small RNAs

While the majority of genic sRNAs target CDS of foreign origin, a number of metazoan CDS are also serving as sRNA targets. It should be pointed out that not all HGTs can be detected by the AI metric—only those derived from nonmetazoan species. We have identified a subset of host genes with AI ≤ −30 (of likely metazoan origin) displaying sRNA coverage over the entire CDS region (a group of 327 genes, see Materials and Methods). As mentioned above, TEs and foreign genes tend to colocalize in subtelomeric regions of the A. vaga genome (Flot et al. 2013). It was therefore of interest to investigate possible colocalization between TEs and host genes with sRNA coverage. This co-occurrence was tested statistically, using the annotated TEs and host genes with AI ≤ −30. First, we counted the number of TEs in windows of 500, 1000, and 2000 bp around each of the selected host genes (a group of 327 genes with sRNA counts) and around host genes without sRNA coverage. Then, we compared the distribution of TE density around host genes with sRNA coverage and around the remaining host genes, using a t-test in three different window sizes. According to the t-test, the density of TEs is significantly higher around the selected group (AI ≤ −30 and sRNA coverage) for genomic windows as small as 500 bp (Table 1).

Another way of determining whether the selected group of host genes is located in the genomic environment similar to that of foreign genes is to employ the number of interstitial telomeric repeats in the vicinity of a gene as a proxy for localization in subtelomeric regions. We counted the number of telomeric repeats in windows of different sizes around annotated TEs, putatively foreign genes, selected metazoan genes with sRNA coverage, and the remaining bulk of metazoan genes. Table 2 shows that, in contrast to the majority of metazoan genes, the sRNA-covered host genes are about as likely to reside in the vicinity of telomeric repeats as are genes of putatively foreign origin.

Host multigene families with small RNA coverage

To find out whether the limited subset of metazoan CDS yielding sRNA coverage can also be regarded as relatively recent, we focused on the sufficiently large multigene family that permits informative phylogenetic analysis, which was identified as an overrepresented category in functional GO annotations (File S1, sheet 3). The CYP450 enzymes catalyze oxidative transformation, which contributes to a broad array of biological functions in living organisms, such as activation or inactivation of many endogenous and exogenous bioactive compounds, particularly xenobiotic detoxification, or even modification of nonribosomal peptides (Guengerich 1991; Haslinger et al. 2015). CYP enzymes are found in all kingdoms and domains of life (metazoans, plants, fungi, protists, bacteria, and archaea) and even in a virus. The number of putatively functional full-length CYP genes in metazoan genomes can reach a few dozen copies (57 in humans, 102 in mice, 84 in D. melanogaster, 74 in Caenorhabditis elegans) (Nelson 2009). After ortholog clustering analysis, which helps to distinguish orthologs and paralogs in large metazoan gene families, we identified 92 A. vaga CYP forming an ortholog group (File S1, sheet 4). Despite the fact that none of the copies appear to yield truncated proteins, they differ in transcriptional activity (RNA-seq counts), and, importantly, their origins are quite dissimilar. BLASTN and phylogenetic analysis (Figure 5A) showed AvCYP genes to be related to different CYP clans previously identified in other taxa. The AI values also indicate that some of them could have been acquired by independent HGT events.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5

Small RNAs in the A. vaga cytochrome P450 (CYP) gene family. (A) The maximum likelihood phylogenetic tree includes A. vaga CYP homologs (in red) and CYP genes from other domains of life. CYP gene sequences and names were taken from the CytP450 homepage (Nelson 2009). Red stars denote genes with small RNA coverage. Genes marked with blue dots indicate probable foreign origin (AI > 0). The node marked with an asterisk (*), which includes most genes with significant antisense read counts (>100), is shown in B. (B) The histogram on the top shows RNA profiles from the AvCYP node marked in A, with antisense sRNA counts and RNA-seq rpk counts. The phylogenetic relatedness and nucleotide sequence identity (%) between nearest neighbors is shown in the dendrogram below. Two pairs grouped with 87.6% identity do not exhibit synteny in their genomic environments. Genes with sRNA reads, marked by red stars in A, are shown in red.

Analysis of small RNA reads and their mapping patterns reveals a patchy distribution (Figure 5A), with one of the expanded clades displaying higher levels of sRNA coverage. This expanded node, marked with an asterisk (16 CYP genes, 12 of which display sRNA coverage), is shown in Figure 5B as a dendrogram with small RNA counts and RNA-seq reads per kilobase (rpk) values per CYP copy. An allelic pair with the highest transcriptional activity lacks sRNA coverage, while in the remaining allelic pairs the transcript levels are much lower, in agreement with increased levels of sRNA coverage. Interestingly, in this group, sRNA coverage is correlated with the presence of an oppositely oriented transcriptional unit in the 3′-flanking region (File S5, sheet 5).