A putative RNA-interference-based immune system in prokaryotes: computational analysis of the predicted enzymatic machinery, functional analogies with eukaryotic RNAi, and hypothetical mechanisms of action

Identification, classification and evolutionary analysis of cas genes

In the original study of the cas gene neighborhoods, which was performed with ~40 genomes, we identified ~20 protein families that were tightly or more loosely associated with the system we now call CASS. A recent update by Haft and coworkers with >200 genomes yielded a diverse "guild" of ~45 Cas protein families [24]. Many of the Cas proteins show very low sequence conservation which makes identification of homologous relationship between them a non-trivial task. We employed an iterative approach to the exhaustive analysis of the Cas protein sequences. The protein sequences of each Cas family were compared to the protein sequences from all available prokaryotic genomes using PSI-BLAST [30, 31], the proteins encoded in the neighborhoods of all identified candidates were used as queries for further searches, and the process was iterated until convergence. The sequences of Cas proteins were, in addition, carefully compared to each other, in an attempt to identify possible traces of common origin of some of these genes that have so far eluded detection [for a complete list of Gene Identification (GI) numbers of the detected Cas proteins, see Additional file 1].

COG1518, a universal marker of CASS

In agreement with the previous observations [23], we found that Cas1 (COG1518 in the Clusters of Orthologous Group of proteins classification system [32]) is the best marker of the CRISPR-associated systems (CASS). This gene encodes a highly conserved protein and that is represented in all cas neighborhoods, with the single exception of Pyrococcus abyssii. A PSI-BLAST search for COG1518 members in the completely sequenced prokaryotic genomes revealed the presence of at least one representative of this COG in 77 of the 177 analyzed genomes.

Given that COG1518 is represented in all versions of CASS, the distribution of this gene among prokaryotic lineages serves as a proxy for the distribution of CASS itself (P. abyssii was included as having CASS, the absence of COG1518 notwithstanding); this distribution is, obviously, non-uniform (Fig. 1). All sequenced archaeal genomes, including, notably, the tiny genome of Nanoarchaeum equitans, encode this protein. In contrast, in each of the bacterial lineages with a representative set of sequenced genomes, half or more species do not have COG1518. We came across several cases of differences in the presence of CASS among closely related bacterial species and even between strains of the same species. Thus, Corynebacterium diphtheriae has a COG1518 gene but Corynebacterium efficiens does not; similarly, CASS was detected in Mycobacterium tuberculosis but not in Mycobacterium bovis or Mycobacterium leprae; significant differences in the composition of CASS were noticed even between strains of several bacterial species, e.g., Thermus thermophilus [33].

Phylogenetic analysis of COG1518 combined with gene neighborhood comparisons yields a clear classification of CASS based on gene composition and (predicted) operon organization (Fig. 2a) which is, generally, in agreement with the recent study of Haft et al [24]. Altogether, we delineated 7 distinct versions of CASS which show substantial differences in gene composition and predicted operon organization. Only four CASS variants (CASS1-4) have a more or less stable operon organization. While a distinct operon organization motif is recognizable in the other CASS variants as well (Fig. 2a), these versions have undergone numerous, complex rearrangements in different prokaryotic genomes including gene duplications, fusion with other CRISPR-related operons, insertion of additional genes, and others. Thus, in line with the observation that CASS is often eliminated from genomes during evolution, this plasticity of CASS operon organization suggests that the CRISPR-associated genomic regions are "hot spots" of recombination and ensuing genome rearrangement. Taken together, all these observations demonstrate the extraordinary evolutionary mobility of CASS.

The recent work of Haft et al has identified a "guild" of over 45 protein families comprising the CASS [24]. Here we refine this classification by unifying many of the families into superfamilies and also expand the list of CASS-linked genes to include those that are found in cas operons less commonly and comprise the "cloud" surrounding the core CASS components. Only COG1518 and COG1343 (see below) are invariably present in all versions of CASS (with the exception of P. abyssii and several bacterial species that carry disrupted COG1518 gene). A few genomes that have a cassette of CRISPR-associated genes without COG1518 also possess another CRISPR-associated operon in a different genomic location that does include this gene. Altogether, we end up with ~25 gene families that are tightly associated with CASS (Table 1)

	Family	SubfamilyA	Phyletic distributionB	Comments
1	COG1518	COG1518 (cas1)	All	Putative novel nuclease/integrase; Mostly α-helical protein
2	COG1343	COG1343 (cas2), COG3512, ygbF-like; MTH324-like; y1723_N-like;	All	Small protein related to VapD, fused to helicase (COG1203) in y1723-like proteins
3	COG1203	COG1203 (cas3)	All	DNA helicase; Most proteins have fusion to HD nuclease
4	RecB-like nuclease	COG1468 (cas4), COG4343	All	RecB-like nuclease; Contains three-cysteine C-terminal cluster
5	RAMP	COG1688, COG1769, COG1583, COG1567, COG1336, COG1367, COG1604, COG1337, COG1332, COG5551, BH0337-like, MJ0978-like, YgcH-like, y1726-like, y1727-like	All	Belong to "RAMP" superfamily, possibly RNA-binding protein, structurally related to a duplicated ferredoxin fold (PDB: 1WJ9)
6	COG1857	COG1857, COG3649, YgcJ-like, y1725-like	All	α/β protein; probable enzymatic activity, possibly, a nuclease
7	HD-like nuclease	COG1203 (N-terminus), COG2254	All	HD-like nuclease
8	BH0338	BH0338-like MTH1090-like	All, mostly archaea and FIRM	Large Zn-finger-containing proteins, possibly, nucleases (nuclease activity has been reported for MTH1090 [75].
9	ygcL	ygcL	Bacteria, mostly PROTEO	Large Zn-finger containing proteins;
10	COG1353	COG1353, MTH326-like, alr1562, slr7011	All, mostly Archaea	Putative novel polymerase; Multidomain protein with permuted HD nuclease domain, palm domain, polymerase-thumb-like domain and Zn-ribbon; MTH326-like has inactivated polymerase catalytic domain; alr1562 and slr7011 – predicted only on the basis of size, presence of HD domain, and location with RAMPs in one operon
11	COG1517/HTH	COG2462	Archaea	Former COG2462; Fusion of COG1517-like domain to HTH-type transcriptional regulator; Possible regulator of the system expression in archaea
12	COG1421	COG1421	All, mostly Archaea	~150 aa protein; Has a few motifs similar to ygcK-like; mostly α-helical protein
13	ygcK	ygcK-like	Bacteria, mostly PROTEO	~180 aa protein; has a few motifs similar to COG1421; mostly α-helical protein
14	COG3337	COG3337	All, mostly Archaea	~110 aa; mostly α-helical protein
15	COG1517	COG1517 COG4006	All, mostly Archaea	Some are fused to HTH domain (see COG1517/HTH), some proteins have the domain duplication; structure is available (1XMX); domain appears to have a Rossmann-like fold.
16	COG3513	COG3513	Bacteria, mostly PROTEO	Huge protein; contains McrA/HNH-nuclease related domain and RuvC-like nuclease domain
17	PH0918	PH0918-like	All, mostly Archaea	Specific for Pyrococcus and Thermococcus. The pair ST0031/SSO1401 and AF1873, most likely, belong to the same family because have similar length and located in the identical place in an operon but due to low conservation are not alignable
18	AF1870	AF1870-like	Archaea	Former COG3574; ~150 aa protein.
19	AF0070	AF0070-like	Archaea	~420 aa protein, no prediction
20	y1724	y1724-like	Bacteria, mostly PROTEO	~450 aa protein, no prediction
21	SPy1049	Spy1049-like	Bacteria, mostly FIRM	~220 aa protein, no prediction
22	TTE2665	TTE2665-like	Bacteria, mostly CHLOR	~130 aa protein, no prediction
23	LA3191	LA3191-like	Few bacteria	~650 aa, no prediction

^A – Subfamilies are named by the corresponding COG number or by a protein ID ^B – All indicates that the family is widespread in all major prokaryotic lineages; PROTEO – proteobacteria; FIRM – firmicutes; CHLOR – Chlorobia

COG1857 – another putative enzyme found in most CASS versions

Another conserved CASS family that resists straightforward functional prediction is COG1857. This gene is present in most CASS, with the exception of the "minimal" variants, CASS 4, 6, and 7. The COG1857 proteins consist of ~350 aa and typically have an N-terminal β-strand followed by a loop containing a conserved asparagine (see Additional file 3). Although no other motifs containing potential catalytic residues seem to be conserved throughout COG1857, perhaps, due to the extensive divergence between the subfamilies of this family, the observed conservation pattern suggests that this protein is an enzyme, perhaps, yet another nuclease. COG1857 proteins are typically encoded immediately upstream of another widespread CASS gene, COG1688 (named Cas5 by Haft et al. [24]). We have previously identified COG1688 as a member of the Repair-Associated Mysterious Protein (RAMP) superfamily [23] that has been renamed to Repeat-Associated Mysterious Protein (preserving the acronym) by Haft and coworkers [24]. In CASS3, COG1688 is replaced by the y1726 family which we found to include another set of RAMPs distantly related to COG1688 (Figs. 2a and 3A).

Ancillary CASS components

The CASS-core and the pol-cassette together comprise the extended central componentry of CASS that consists of six COGs plus the RAMP superfamily (Table 1). In addition, ~20 other families were consistently found in CASS, even if in a minority of the CASS-containing species (Table 1). It has to be taken into account that most of the CASS proteins have highly diverged sequences, and further unification as well as expansion and the ensuing rise in status of some of the families remains a possibility. For example, COG1517 is such a growing family. Typically, the COG1517 genes are located on the periphery of pol-cassettes and, in many cases, contain fused helix-turn-helix (HTH) domains. The core domain of COG1517 is ~150 aa; several proteins, e.g. TM1812, contain a duplication of this domain, and several fusions, in addition to the one with HTH, were detected. The sequences of COG1517 proteins tend to be highly diverged, and some of these were detected in genomes with no CASS, indicating that the association of this gene with CASS is not as tight as it is for the extended core of CASS. For one such protein from Vibrio cholerae, the 3D structure has been solved, again, in one of the structural genomics projects (PDB: 1xmx). The structure is composed of three domains. The N-terminal domain (A-1-A137, β1 through β7, Fig. 4) is a modified Rossmann-like fold, with the next to the last β-strand (β6) being antiparallel to the rest of the strands. The middle domain (A161-A249, α6 through α10) is a winged helix-turn-helix (HTH) with an additional α-helix on each terminus (α6 and α10, respectively). The C-terminal domain (composite of A250-A383 and A138-A156, α11 through α14 plus α5 and β8) belongs to the restriction endonuclease superfamily. In the 1xmx structure, the endonuclease domain is positioned between the HTH and the Rossmann-like domain (Fig. 4). The lowest VAST P-values of the matches between each of the 1xmx domains and the respective domains with known folds are highly statistically significant (<10^-5) which is strongly suggestive of homologous relationships.

The functions of several other CASS gene families remain obscure. For instance, CASS1, 5, 7 contain genes upstream of COG1857 that encode large (500–600 amino acids), homologous proteins; the best conserved family in this set of proteins is represented by BH0338 and its orthologs. Some members of the BH0338 family contain a Zn-ribbon in the middle of the sequence but otherwise have no recognizable domains or motifs. Among several conserved motifs of these proteins, are two conserved aspartates and a distal conserved glycine, a combination that resembles the motifs seen in the PALM polymerase domain (not shown). Although we were unable to obtain additional evidence of the potential connection of this family with polymerases, it is tempting to speculate that these proteins might contain an extremely diverged version of the PALM domain. Similar, albeit less pronounced, motifs are detectable in the MTH1090 family proteins which are present in CASS5 and CASS7. This similarity and the fact that the respective genes occupy the same position in the corresponding operons suggest that the BH0338-family and the MTH1090-family proteins are highly diverged homologs. CASS2 also includes a large protein (YgcL-family), some of which contain Zn-clusters; however, the conserved motifs of these proteins and their position in the respective operons are different from those of the MTH1090 and BH0338 families. CASS4 contains another huge protein (COG3513, ~1150–1400 aa) with two recognizable domains, a McrA/HNH nuclease and a RuvC-like nuclease (RNAseH fold). These observations emphasize the striking diversity of still poorly characterized CASS components, particularly, the plethora of predicted nucleases of various classes and potential novel ones.

Hypothesis: CASS is a prokaryotic defense system that functions on the RNAi principle

Based on the properties of CRISPR and Cas proteins, we speculate that this system is a functional analog of the eukaryotic siRNA systems and propose possible mechanisms of the putative prokaryotic small RNA interference. The crucial observation reported independently by Mojica et al [27], Pourcel et al [29], and Bolotin et al [28] is that a certain fraction (~10% according to [27]) of the unique inserts in CRISPR units are homologous to fragments of viral (bacteriophage) or plasmid genomes. Only a miniscule fraction of the existing phage and plasmid sequences is currently available, whereas the total diversity of prokaryotic mobile elements is humungous and apparently exceeds the diversity of prokaryotes at least by an order of magnitude [41, 42]. Thus, it is not far fetched to propose that most, if not all, CRISPR inserts are derived from mobile replicons [29].

Should that be the case, it seems more or less obvious that CASS is a prokaryotic defense system against foreign replicons that functions on the antisense RNA principle. More specifically, it seems likely that the inserts are transcribed and silence the cognate phage or plasmid genes via the formation of a duplex between the prokaryotic small interfering (psiRNA) and the target mRNA followed by cleavage of the duplex or translation repression. Indeed, Mojica et al [27] mention the analogy between the CRISPR and eukaryotic RNA interference systems but propose no specific mechanisms for the action of the putative defense systems and, crucially, do not explore the connection between the putative psiRNA and the predicted activities of Cas proteins. Important supporting evidence has been independently obtained through the analysis of the small non-messenger RNA expression in the euryarchaeon Archaeoglobus fulgidus which showed that CRISPR are transcribed (from a leader-promoter sequence), apparently, in the form of a multiunit precursor that is subsequently cleaved into CRISPR monomers and oligomers [18]; very similar observations have been subsequently reported for the crenarchaeon Sulfolobus solfataricus [19]. Furthermore, as noticed by Pourcel and coworkers [29], one of the unique CRISPR inserts in the MIGAS strain of the bacterium Streptococcus pyogenes is homologous to a prophage present in other strains of the same bacterium that, conversely, do not carry the CRISPR. This is compatible with the possibility that the insert makes the bacterium immune to the given phage.

Figure 5a presents our current hypothesis on the basic mechanism of the functioning of the prokaryotic RNAi system. We speculate that CRISPR regions are transcribed from a promoter located in the AT-rich CRISPR leader. Transcription might be regulated with the participation of one of the Cas proteins and, conceivably, would be stimulated by stress, such as phage infection, although the results of expression analysis in A. fulgidus suggest some level of constitutive transcription [18]. The work of Tang and coworkers further indicates that the primary transcript is likely to encompass the entire CRISPR repeat region. This transcript would be cleaved into 70–100 nt pre-psiRNA, conceivably, by the putative p-dicer, the COG1203 protein (Fig. 5a). We further postulate that p-dicer catalyzes the second, perhaps, slower (judging by the results of Tang et al) processing step that releases mature psiRNA species (Fig. 5a). The psiRNA molecules would bind RAMPs in a size-specific manner and anneal to the target mRNA. The resulting complex would recruit p-slicer, forming the minimal form of the prokaryotic analog of RISC (pRISC) that would cleave the mRNA and could be recycled to attack the next target molecule, thus, silencing the respective gene (Fig. 5a).

Figure 5b shows a version of this pathway that involves the activity of the CASS polymerase, by analogy with the eukaryotic RdRp, which participates in some RNAi pathways in most eukaryotes, but apparently has been lost in arthropods and chordates [45–47]. The initial steps in this scheme are the same as in the basic one (Fig. 5a) – transcription of the CRISPR, processing of the psiRNA precursor, and annealing of psiRNA to the target mRNA mediated by the pRISC – but, at the next step, psiRNA is postulated to serve as the primer for elongation by the CASS polymerase, yielding an extended double-stranded form of the target (Fig. 5b). This form would be cleaved by p-dicer analogously to the cleavage of viral and transposon dsRNAs by the eukaryotic dicers. The p-dicer might function as a complex with the respective RAMP to form a distinct version of pRISC. This could be the endpoint of the pathway, or else, the dsRNA degradation products could be utilized as new psiRNAs, resulting in amplification of the silencing effect (Fig. 5b). The CASS polymerase is most common in thermophiles, and it is tempting to speculate that the prevalence of this form of the psiRNA pathway has to do with the instability of the psiRNA-target duplex under the high ambient temperatures of these organisms.

The most complex and uncertain aspect of the putative prokaryotic RNAi system discussed here is the formation of new CRISPR units containing unique psiRNA genes specific for new targets encountered by the organism (Fig. 5c). The path to the creation of new psiRNAs would begin just like the response pathway, i.e., with transcription of the CRISPR locus and the first processing step yielding the 70–100 nt psiRNA precursors (compare Fig. 5c with Fig. 5a). At the next step, however, there must be a mechanism to replace the unique insert within the pre-psiRNA with a new fragment of foreign (e.g., phage) RNA. The nature of this mechanism remains unclear. In principle, two possibilities can be envisaged: i) reverse transcription with copy choice whereby a reverse transcriptase, most likely, the predicted CASS polymerase (COG1353) switches from using the pre-psiRNA as a template to using a phage mRNA, and then back, and ii) direct, non-homologous RNA recombination between a pre-psiRNA and a foreign mRNA, followed by reverse transcription of the resulting recombinant RNA (Fig. 5c). Both mechanisms are non-trivial in their molecular choreography and are unlikely to occur with high efficiency. Nevertheless, there are precedents for both in the molecular biology of retroviruses and other RNA viruses. In particular, reverse transcriptase switches templates in each cycle of retrovirus first-strand cDNA synthesis although, in this case, copy-choice is facilitated by the spatial juxtaposition of the two templates within the virus particle; a similar mechanism is responsible for recombination in retroviruses [48]. In addition, and probably, more relevantly to the psiRNA case, reverse transcription with copy-choice is thought to be involved in the incorporation of copies of cellular genes, such as oncogenes, into retroviral genomes [49, 50]. The alternative, namely, direct recombination between RNA molecules might seem far fetched, but such a process has been demonstrated, by several groups independently, to occur in RNA viruses, apparently, via a protein-independent mechanism [51, 52]. During the formation of new psiRNA species, these low-frequency processes might be facilitated by the high abundance of the phage mRNAs involved. Indeed, it has been shown that the unique inserts in CRISPR most often correspond to fragments of essential, highly conserved phage genes that are typically expressed at a high level in infected bacteria[27]. Once the dsDNA molecule consisting of a CRISPR unit with the new, unique insert is produced, by one of the mechanisms outlined here or, perhaps, via a different pathway, it must insert into the CRISPR array via homologous recombination (Fig. 5c). We suspect that this process is mediated by the COG1518 protein, the universal marker of CASS containing conserved motifs resembling those of different nucleases [23]. It seems likely that this protein functions as the CRISPR integrase/recombinase, perhaps, in cooperation with the COG1343 protein, another universal component of CASS.

Additional lines of evidence relevant for the predicted RNAi function of CASS

A connection between CRISPR and RAMPs

We searched for CRISPR in the 89 complete prokaryotic genomes containing cas genes and found remarkable diversity, in both the number of repeat units and the spacer size, even among closely related species and strains, which is in agreement with previous observations [25, 27–29]. Given our hypothesis that RAMPs are size-specific psiRNA-binding proteins, we examined possible connections between the number of RAMP genes, the number of CRISPR units (which is the same as the number of unique inserts) and the length heterogeneity of the unique inserts in prokaryotic genomes. The number of RAMPs showed strong positive correlations both with the number of CRISPR units (Fig. 6a) and with the variance of the insert lengths (Fig. 6b) which seems to be compatible with our hypothesis on the psiRNA-binding function of RAMPs. Significant correlations were observed also between the number of CRISPR and the variance of the insert length, and between each of these variables and the number of COG1518 proteins, COG1517 proteins or all Cas proteins minus the RAMPs. However, the correlations seen for the RAMPs and the CRISPR inserts were the strongest, suggesting that these were the biologically most relevant links (Table 4). It should be emphasized that, given the extreme sequence divergence of the RAMPs, it is most likely that the correlation between the number of encoded RAMPs and the number of CRISPR units present in a genome is due to the RAMPs discriminating between inserts (by size, under our hypothesis) rather than recognition of the repeats themselves.

Values selected for correlation	Rank Correlation coefficienta
Number of spacers vs S.D. of spacer lengths	0.398
Number of spacers vs number of COG1518 proteins	0.475
Number of spacers vs number of RAMPs	0.668
Number of spacers vs number of COG1517 proteins	0.524
Number of spacers vs total number of Cas proteins minus RAMPs	0.580
Number of spacers vs total number of Cas proteins minus RAMPs and COG1517 proteins	0.556
Number of RAMPs vs S.D. of spacer lengths	0.651

^aAll correlation coefficients are highly statistically significant P < 10^-3

Putative psiRNAs: relationships with phage and plasmid genes and secondary structure

Importantly, the CRISPR insert sequences from even closely related bacterial strains are unrelated to each other, with only a few exceptions among enterobacterial strains (data not shown). This suggests that the inserts are replaced rapidly on the evolutionary scale, perhaps, via the mechanism outlined in Fig. 5c. A broader implication is that the dominant phages and plasmids encountered by even closely related bacteria are different leading to the rapid generation of distinct repertoires of psiRNAs.

In terms of their ultimate function, the putative psiRNAs appear to be conceptually most similar to eukaryotic siRNAs in that they are homologous to foreign, rather than endogenous, genes and are predicted to function as part of a defense system. However, the psiRNAs also share a major feature with the miRNAs in that they are embedded in and, probably, cleaved from a specific precursor molecule, rather than from a long dsRNA like the siRNAs (Fig. 5). In an attempt to elucidate common features of the psiRNA precursors, we predicted the secondary structures of all available CRISPR units. The results were somewhat ambiguous as shown in Fig. 7. The folding free energies of GC-rich CRISPR units were substantially lower than those for the shuffled CRISPR units; notably, the CRISPR folding energies were distributed similarly to those of miRNA precursors (Fig. 7a), with numerous CRISPR units capable of folding into highly stable secondary structures. These observations suggest the possibility that psiRNA precursors possess a distinct secondary structure that could be important for recognition by p-dicer, RAMPs or other proteins. However, the folding energy distribution for AT-rich CRISPR shows only a slight shift compared to the distribution for shuffled sequences and is completely different from the distribution for miRNA precursors of the same base composition (Fig. 7b). Thus, AT-rich psiRNA precursors, which comprise approximately half of the CRISPR, on average, do not seem to be capable of folding into stable secondary structures. Thus, the existence of a functionally relevant consensus structure of psiRNA precursors remains uncertain.

Figure 8 shows two predicted secondary structures of putative psiRNA precursors in which the inserts are homologous to phage or plasmid genes. The relatively GC-rich sequence from Xanthomonas axonopodis forms a stable secondary structure, whereas the AT-rich sequence from Streptococcus thermophilus can fold only into a weak stem-loop structure. Curiously, in both cases, the unique insert contributes significantly to stem formation. Clearly, further analysis of the structures of putative psiRNA precursors by computational and experimental means is required to identify their functionally important features.

Genome sequences, databases and sequence analysis

The complete bacterial and archaeal genome sequences were retrieved from National Center for Biotechnology Information (NCBI, NIH, Bethesda) FTP site. The non-redundant database of protein sequences at the National Center for Biotechnology Information (NIH, Bethesda) was iteratively searched using the PSI-BLAST program [31]. The cut-off of E < 0.01 was normally employed for inclusion of sequences in the position-specific weight matrices. Each retrieved sequence was used as the query for additional searches until no new sequences could be detected. For detecting subtle sequence conservation, the PSI-BLAST search results were visually examined and sequences with greater E-values, but containing signature motifs of a given protein family were included into profiles on a case by case basis [30, 60]. Multiple alignments of protein sequences were constructed using the MUSCLE program [61] and corrected on the basis of PSI-BLAST results. Protein secondary structure was predicted using the JPRED program [62]. Protein structure comparisons were performed using the DALI [63] and VAST [40] programs, and ribbon diagrams of protein structures were generated made using the program BOBSCRIPT [64].

Phylogenetic analysis

Distance trees were constructed from multiple protein sequence alignments after excluding poorly aligned positions, by using the least-square method as implemented in the FITCH program of the PHYLIP package [65, 66]. Maximum likelihood trees were constructed using the ProtML program of the MOLPHY package, with the JTT-F model of amino acid substitutions, by optimizing the least-square trees with local rearrangements [67, 68].

Identification and analysis of CRISPR repeats

Search for repeats was performed as follows: first, all exactly matching 20-mer anchor substrings were identified in the nucleotide sequence of a bacterial or archaeal genome. Alignments around these anchors were expanded in both directions, to include all adjacent positions with information content above or equal to 1.5 bits [69]. All identified high-similarity fragments were used as queries in nucleotide BLAST [70] search (word size 7; mismatch penalty -1; both gap opening and extension costs -1; E-value threshold 0.001) to detect more diverged versions of the repeats as well as the instances of the repeat in the opposite strand. To determine the repeat boundaries more precisely, sequences of all loci with short (20 nt) flanks added were collected and aligned using the MUSCLE program [61]. Alignments were trimmed from the 5' and 3' termini up to columns with information content exceeding 1.3 bits. Repeat families that shared the spatial arrangement typical of CRISPR (median repeat length of 20–50 nt; median spacer length of 15–60 nt) were identified as CRISPR candidates and further examined for chromosomal proximity to cas genes. The custom PERL scripts used for this analysis are available upon request.

Sequences of eukaryotic miRNAs precursors, CRISPRs, and randomly shuffled CRISPR sequences were computationally folded, and the free energy of the most stable secondary structure was calculated using a a dynamic programming algorithm that employs nearest neighbor parameters to evaluate free energy [71]. Energy minimization was performed by dynamic programming method that finds the secondary structures with the minimum free energy by summing up the contributions from stacking, loop length, and other structural features, using improved thermodynamic parameters [72].

Similarity of inter-CRISPR spacers to other sequences

Nucleotide sequences of inter-CRISPR spacers were used as a query in MEGABLAST [70] searches (word size 11; e-value threshold 0.01) against GenBank; hits to virus or plasmid sequences and to distantly related prokaryotes were counted separately for each source organism.

Reviewer's report 1

Eric Bapteste, Department of Biochemistry and Molecular Biology, Dalhousie University, Halifax, Nova Scotia, Canada

Main comments:

The scientific quality of this paper and its methodology is certain. There has been a lot of good and interesting work done here. As indicated by its title, it thus provides multiple information and one hypothesis about RNA-interference in prokaryotes. Because of this broad scope, the manuscript is quite large. In fact, several of its parts could be read on their own, depending on the reader-specific interest, this is notably the case for the part dealing with the hypothesis of a RNAi prokaryotic immune system. I would thus suggest that a shorter version of the paper, centered around this very interesting hypothesis could be proposed online (with the first part turning into Supp. Mat.), because I feel that this part is going to receive more attention anyway, and it would be unfortunate if some readers did not consider this aspect because they are scared by the overall size of the paper. But this is simply a suggestion, and the authors are more than welcome to disregard my opinion.

Author response: While we fully understand the sentiment and agree that the hypothesis on prokaryotic RNAi is of greater general interest than the detailed presentation of the protein sequence-structure analysis, we strongly feel that the latter provides a badly needed foundation for the hypothesis as well as important information in its own right. Furthermore, we tend to believe that the general spirit of online publication is to present complete results of a study (of course, there are exceptions). The reader can easily navigate between sections, so the length of a paper does not represent a particularly severe problem. Furthermore, we made certain modifications to the protein analysis part in response to similar but more specific comments of Martijn Huynen, in particular, introduced additional subheadings which, hopefully, makes this part of the paper more reader-friendly.

The seductive hypothesis of a RNAi based immune system is presented as an analogy with the eukaryotic RNAi system. The use of analogy is potentially challenging: on the one hand, it allows a powerful and elegant presentation of many complex genomic results, but on the other hand, it is questionable, since the analogy may impose an a priori model to interpret biological features, and if this model is incorrect, if the analogy does not hold, there is a risk that the genomics data receive a fairly biased interpretation. In this respect, it would be interesting if the authors discuss whether an homologous immune system could have been possible in eukaryotes and in prokaryotes, and why it is not found. Indeed, such an homologous system would be a more natural reference to interpret the data than an analogy.

Author response: We understand the epistemological concerns regarding the role of analogy in this study. However, as indicated by the reviewer himself, the analogy is strong. Moreover, this analogy is manifest at two different levels: i) the presence of inserts homologous to phage and plasmid genes in CRISPR units and ii) presence of predicted activities compatible with a siRNA-like system among cas gene products, in particular, the dicer analog. Had the analogy been false, there would be no reason whatever for this congruence. We briefly comment to that effect in the revised manuscript. The idea of homologous immune systems in prokaryotes and eukaryotes seems a little far-fetched. Nevertheless, this comment prompted us to incorporate a brief comparison of the evolutionary histories of RNAi and classical immune systems.

This being said, I do not feel that the use of the analogy was a problem here, as it is convincingly presented and argued by the authors. We could eventually question more if all the so-called CASS genes are really involved in the prokaryotic immune systemand deserve their label: some may just be present in the genomic proximity of CRISP, yet having nothing to do with the RNA interference. They might just be mobile «travelling» genes. A further study of the genomic distribution of the homologs of the CASS genes in bacterial genomes may help to clarify which genes are strongly and exclusively CRISP related and which ones can be found also in alternative locations in different genomes. Then, perhaps the «striking diversity of still poorly characterized CASS components» described on page 16 would appear less striking if some CASS categories are simply not relevantly defined, and include unrelated proteins, since maybe, the use of the analogy would had led to too relaxed definitions of CASS. In another situation, by contrast, the use of the analogy could perhaps be too strict. On page 18, the authors wonder how to identify «the slicer counterpart (p-slicer)» in prokaryotes. They explain that this identification is «less straightforward because of the diversity of predicted nucleases within the CASS». But, after all, why should there be only one p-slicer, as in eukaryotes? It is possible that prokaryotes have multiple «p-slicers».

Author response: Yes, we agree, the possibility of multiple slicers exists, and we modified the text to acknowledge this. With regard to the rest of this comment, however, we feel that the current diversity of prokaryotic genomes is already sufficient to make conclusions on the strength of the association of individual genes with CASS, and the genes are classified here accordingly, as true CASS components and loosely associated "satellites". As far as the latter are concerned, inferences on involvement with CASS functions are made only for those genes whose activities appear clearly relevant, like the RT or Argonaut.

Finally, the strong suspicion about the analogous and intricate functions of COG1518 and COG1343 (cf. page 21) could be similarly toned down. Maybe these genes do play the essential analogous role of CRISPR integrase/recombinase consistently with the analogy, but maybe they fulfill several different tasks. Perhaps the authors would like to comment more on some of these minor points.

Author response: Actually, this prediction is not based on analogy with eukaryotic RNAi systems but rather on the mutualistic association of these genes to CRISPR and the features of the proteins themselves. We dropped the "strong" suspicion but, generally, we strongly believe that this is the best possible prediction. Reference to multitasking might not be particularly productive unless there are good ideas regarding what these multiple functions might be (this is different from the above possibility that there are multiple slicers which is, indeed, compatible with the data).

Also, to go back to the CASS gene evolution, the authors mention, page 8, « the extraordinary evolutionary mobility of CASS». It is unclear to me how this statement has been tested, and how the authors have established that CASS genes are more mobile than any average gene randomly picked in the same collection of prokaryotic genomes. For this reason, I am not sure if, as claimed by the authors, on page 14 «the pol-cassette comprises a distinct evolutionary unit that is often transferred horizontally independently of the CASS-core». Does the CASS-core really have an established vertical mode of inheritanceor, as the authors stated before, a «non-uniform» distribution (cf. page 8)? This might be more strongly argued.

Author response: Several distinct issues are addressed here. With regard to the 'extraordinary mobility' of CASS, this is demonstrated by the trees in Fig. 2(more trees have been published previously in our own 2002 paper and by Haft et al.), but even more convincingly, by the persistent pattern of presence-absence of CASS in closely related species and even strains of bacteria. We consolidated the argument such that this becomes clear the first time "extraordinary" mobility comes up. It is true that we did not compare the mobility of the CASS components with that of garden-variety prokaryotic genes in a rigorous, quantitative manner. While this is doable, in principle, all methods we are aware of are open to debate, and we feel that the exercise is beyond the scope of this paper. Given the above argument, we believe that, qualitatively, it is clear that CASS is unusual in this respect. With regard to the pol-casette, we believe that the discrepancy between the topologies of the two trees in Fig. 2is quite sufficient for the statement on independent HGT. As for the "vertical mode of inheritance" of CASS, there seems to be a semantic issue here. We do not really claim vertical inheritance for CASS but neither is such a pattern necessary to detect horizontal mobility. What is required is a predominant pattern of vertical inheritance among other genes that allows us to use a species tree to detect HGT. Of course, we realize that there are substantial arguments for abandoning "tree thinking"[73]altogether but, on balance, we still believe that a species tree conceptualized as a central trend in the evolution of gene ensembles [74]is, at least, a useful tool for analysis of genome evolution.

These few questions show a strength of the present work which interestingly opens perspectives and suggests that some additionnal analyses should now be conducted, because the topic deserves consideration. Maybe the authors would feel like addressing some of the points below in a revised version of the current paper, or in future analyses.

Author response: These are, indeed, very interesting questions, we appreciate them. Some are for future studies but we can provide certain answers now.

Further study could iclude the following:

- Do other genomic regions harboring concentrations of nucleasescomparable to the ones around CRISP exist elsewhere in the genomes?

Author response: Hardly. As indicated in the paper, the CRISPR neighborhood is the second most prominent (i.e., the one that ranks second in the number of genes) neighborhood in prokaryotic genomes after the ribosomal superoperons, so it is quite outstanding. However, there are other, considerably small constellations of nucleases, such as the classical recBCD operon encoding repair proteins, some restriction-modification systems, and, perhaps, others that are still poorly understood and deserve investigation.

- If yes, is there more than one prokaryotic immune system definable on this analogous ground? Notably, did bacteria without CRISPR evolve a totally different immune system?

Author response: There is no evidence of that. Furthermore, as repeatedly emphasized in this paper, CASS shows extreme evolutionary volatility, apparently being lost quite easily, in a short time, on evolutionary scale. It is hardly imaginable that these bacteria evolved a distinct immune system in the short time elapsed since the loss of CASS. Of course, purely hypothetically, one could perceive the possibility that another immune system is disseminated horizontally, like CASS, and prokaryotes having both, could differentially lose one of them. However, we are unaware of any support for such a scenario. Another prominent prokaryotic defense mechanism is restriction-modification; it would be interesting to examine the relationship between RM systems and CASS, that could be a subject for a future study.

How did the psiRNA pathway arise in thermophiles(cf. page 20)? Does it result from a transfer? Was it ancestral?

Author response: Very interesting, fundamental questions, indeed. In response, we expanded the discussion of these and other aspects of evolution of CASS. The specific preponderance of CASS in thermophiles, noticed already in the 2002 paper, when we thought that this was a thermophile-specific repair system, remains a mystery. Whatever the nature of this association, it seems likely that CASS is ancestral in thermophiles (at least in hyperthermophiles).

- Could we imagine that multiples promoters exist, both sense and antisense, which would activate the transcription of CRISPR, generating even more RNAi(cf.p 24)?

Author response: In principle, existence of multiple promoters cannot be ruled out. However, the leader sequence seems to be the only natural candidate for the promoter function. The rest of the CRISPR cassette is homogeneous (repetitive), so it is unclear where an alternative promoter would be located. Further, in the two archaeal systems that have been studied experimentally (Archaeoglobus and Sulfolobus) all transcription of CRISPR loci appear to be unidirectional.

- Finally, it might be challenging, though interesting to test in vitro on bacterial cultures if, as proposed by the authors, the presence of CRISP and CASS, has really an impact on the fitness of prokaryotes in presence of viruses.

Author response: We certainly hope that the computational analyses and predictions described in this paper stimulate a lot of experimentation aimed at elucidation of the biological functions of CASS and roles of its individual components.

We greatly appreciate these insightful and stimulating comments.

Minor comments/questions:

On page 7: «functionally analogous» is redundant.

Author response: We see the point but do not really agree. The word "functionally" seems to add clarity.

On page 8: the sentence «the distribution of COG1518 and, by implication, CASS among prokaryotic lineages...» is too «bold» for me: even if the conclusion is correct, I am not sure one can generalize as suggested here from the case of one protein only.

Author response: Indeed, we can. Rephrased to clarify and emphasize this.

On page 15: «several other CASS gene families remain mysterious» is a mysterious sentence. I am not sure what this does really mean.

Author response: That there is no clue as to the possible functions of these proteins; modified to clarify.

On page 21: I miss the idea of the sentence starting by «In addition, and probably, more relevantly etc.» to «retroviral genomes». Could you rephrase it to explicit it a little bit more?

Author response: Rephrased – hopefully, to clarify.

On page 23: what is the criterion retained for homology between the plasmid genes, fragments of phages and the CRISPR sequences?

Author response: The following quote from the Methods addresses this issue:

"Nucleotide sequences of inter-CRISPR spacers were used as a query in MEGABLAST[70]searches (word size 11; e-value threshold 0.01) against GenBank; hits to virus or plasmid sequences and to distantly related prokaryotes were counted separately for each source organism."

On page 48: To me, the multiple positive correlations evoke multiple causalities and the possibility of some hidden correlations. Would you say that all the relevant combinations have been considered here?

No, we won't claim that. More complex multiple regression analysis would be required to separate correlations that reflect true causality; for the purposes of this paper, we felt it was sufficient to note the strongest correlations.

Reviewer's report 2

Patrick Forterre, Biologie Moléculaire du Gène chez les Extrêmophiles (BMGE) Institut de Génétique et Microbiologie (IGM), Université Paris-Sud, Centre d'Orsay, 91405 Orsay Cedex, France, and Biologie Moléculaire du Gène chez les Extrêmophiles (BMGE), Département de Microbiologie Fondamentale et Médicale, Institut Pasteur, Paris, France

In this very important paper, Makarova and coworkers propose a detailed mechanism for a putative procaryotic antiviral immunity system mediated by CRISPS sequences and their associated Cas proteins (the CAS system, CASS sensu the authors). Their model is based on the hypothesis that these elements represent a prokaryotic-specific antiviral mechanism analogous to the eukaryotic RNAi system. In procaryotes (Bacteria and Archaea) there is no homologs of the proteins involved in the eucaryotic RNAi system. Untill recently, it was therefore widely believed that restriction-modification mechanisms were the only defense available to procaryotes to fight viral infections. However, it has been proposed last year by several groups that procaryotic CASS could also play a significant role in fighting viral aggression in archaea and bacteria (Mojica et al. 2005, Pourcel et al., 2005, Bolotin et al., 2005). CRISPR sequences, which are transcribed but non-coding, are formed by the tandem repetition of units containing both a conserved element (similar all along a given CRISPR) and a variable element, the spacer, different from one unit to the other. The spacer sequences have strikingly no homologous sequence in databases, except for viral or plasmid sequences. Both Mojica et al. (2005) and Bolotin et al., (2005) have suggested that transcription of CRISPR sequences produce anti-sense RNA that can inhibit transcription of incoming viral (plasmid) sequences and Mojica et al. (2005) mentioned the analogy of such system with eukaryotic RNAi. However, these authors did not elaborate on the specific mechanism involved and how the cas proteins could be involved in the processing of viral RNA.

In this work Makarova and co-workers have first performed an updated analysis of cas proteins using genomic context analysis and sensitive methods (iteration approaches) to detect low level of similarity and to classify cas proteins in families and superfamilies. They were able to identify several new putative cas proteins and to define 25 superfamilies of cas proteins and 7 different types of CASS organization (named CASS1 to 7). They have also analyzed all available CRISPR repeated sequences and their putative secondary structures. More importantly, they try to predict the biological function of the cas proteins and their mechanism of action in the framework of the RNAi hypothesis. Previously, it has been suggested that cas proteins were involved in the formation and spreading of the CRISPR. For instance, Bolotin et al. Have predicted that cas proteins are acting at the DNA level by promoting cleavage, recombination and ligation. Makarova and al are the first to suggest that several cas proteins should instead interact at the RNA level, by promoting RNA degradation and RNA-RNA hybridization. They specifically suggest the existence of procaryotic homologs of eucaryotic dicer (helicase-nuclease) and splicer (nuclease). They also propose that a previously suspected DNA polymerase could be an RNA dependent RNA polymerase used to stabilize RNA/RNA hybrid by extending iRNA hybridized to their viral mRNA target. They also suggest the involvement of a reverse transcriptase in the formation of the linker sequences from viral (plasmidic) mRNA. In my opinion, all these proposals are reasonnable and very convaincing. Another prediction is that RAMP proteins recognize linker sequences of different sizes. This is supported by a correlation between the number of linker sequences and the number of RAMPs encoding genes (Fig. 6). In that case, it's not clear to me why this could not be due to the binding of RAMPs to the repeated units, since these units exhibit conserved sequences and their number (identical to the number of linkers) should be also correlated with the number of RAMPs.

Author response: That RAMPs discriminate, one way or another, between CRISPR inserts, is strongly suggested by the extreme sequence divergence of RAMPs which is hardly compatible with recognition of identical repeats. To be explicit about it, we added a clarification at the end of this section.

The search for specific secondary structure associated to the repeated units did not give convincing results and suggest for me that the dyad symmetry observed in many repeat units could be due to the binding of proteins with repeated structure (possibly the duplicated ferredoxin-like fold present in RAMPs) and not the formation of secondary structures in the transcribed repeats.

Author response: It is hard to see how one excludes the other: it stands to reason that CRISPR do form distinct secondary structures which bind to symmetrical proteins.

The model proposed (including possible variation) thus implies many predictions that could be experimentally tested. Surprisingly, to my knowledge, only one cas protein has been studied at the bench up to now (ref 70 in the manuscript). This protein turns out to have DNAse activity in vitro, but I suspect that the authors have not tested a possible RNAse activity. This is surprising because the importance of these proteins was already highlighted in 2002 by two in silico papers that in one case suggested their participation to a "mysterious DNA repair system and in the other described their association with CRISPR sequences. The present paper, with much more specfic predictions, should hopefully strongly stimulate biochemists and molecular biologists to jump onto this really exciting story. As noticed by the authors in their conclusion, if their hypothesis turned out to be correct, this prokaryotic RNAi system could be exploited to silence any gene in organisms that encode CASS. Furthermore, the experimental study of this system should help us to get new critical insights on the dynamic relationships between viruses and archaeal/bacterial populations in nature.

Finally, I would like to know if the authors have some idea about the origin of this CAS system. Why is it present in all archaeal genomes sequenced so far? Is it possible that this system originated in Archaea and was later on introduced in bacteria by LGT?

Author response: Given the horizontal mobility of CASS, we can only speculate on the point of its origin. We expand such speculation in the revised conclusion including the possibility of archaeal origin.

– In some case, the authors should be more cautious in their statement. For instance, when they talk about the pol-cassette, it might led some reader to believe that the polymerase actvity of the COG1353 protein has been experimentally validated, which is not the case.

Author response: We added a few more "predicted". However, we did not want to abandon the term 'pol-cassette' as it is descriptive and succinct.

Reviewer's report 3

Martijn Huynen, Nijmegen Center for Molecular Life Sciences University Medical Center St. Radboud p/a Center for Molecular and Biomolecular Informatics, Nijmegen, Netherlands

This paper provides a highly interesting and well documented hypothesis about a cluster of genes that E. Koonin and co-workers have discovered some time ago. By combining biological knowledge with bioinformatics methods and creative thinking the authors propose that Archaea and to a bit lesser extent Bacteria posses an RNA-interference-based immune system involving CRISPR and cas genes, that is analogous the eukaryotic RNA interference systems. Although aspects of this hypothesis have been published before, specifically with respect to CRISPR, this paper is, as far as I can tell, the first that makes the analogy between the cas genes and the RNA interference system. The idea that prokaryotic genomes would internalize pieces of foreign DNA in order to be able to defend themselves against it, thus having an immune system with a memory, would be an interesting example of Lamarckian evolution.

I do have some questions and editorial comments that I think should be addressed.

1) Do the authors have any idea why this system has the phylogenetic distribution that it does, being present in such a small genome as the nanoarchaeon, but not in e.g. the majority of Firmicutes

Author response: No mechanistic idea, unfortunate as this might be. We added some additional discussion of the ultimate origin of CASS (see the response to Patrick Forterre).

2) concerning the feasibilty of the system proposed by the authors: Is there anything known about how many fiendly DNAs a prokaryote encounters in daily life, and how does that compare to the number of different elements in a CRISPR ?

Author response: Not enough for this particular comparison. However, it is well known that phages are extremely abundant, much more so than bacteria or archaea, and in the revised manuscript, we refer to this more specifically, with the corresponding references.

3) Regarding the Lamarckian scheme: That the unique element of the CRISPR correspond to highly conserved, essential elements of phage genes suggests that selection on genetic variation also plays role here. So the scheme would be partly Lamarckian.

Author response: Probably, so. The way we state it in the text "CASS seems to come closest to a true Lamarckian mode of evolution among all known systems of heredity" is compatible with this view.

4) I am not so convinced by the argument on page 3 that the results imply that even among closely related prokaryotes the most commonly encountered phages are different. First of all, it is more a corollary of the hypothesis, but second, it could also reflect the high turnover of phages over time, rather than niche.

Author response: This is a very good idea, we now mention this possibility both in the Abstract and in the Discussion.

5) I am puzzled on the involvement of more or less randomly selected pieces of DNA from foreign DNA/RNA in exactly the same location in the secondary structure of the psiRNA (the top of the hairpin). Does this pattern occur more often?

Author response: The situation when the insert forms a stem of varying stability with parts of the repeats is common but not universal. The positions of the inserts are not exactly the same although they are, indeed, very similar, and the stems in which the inserts are involved are imperfect. Of course, the exciting possibility exists that the CRISPR inserts are specifically selected for their ability to base-pair with the repeats, however, we do not have enough data to make that claim.