The independent loss model with ordered insertions for the evolution of CRISPR spacers

TL;DR

This paper models the evolution of CRISPR spacer arrays using an independent loss model with ordered insertions, deriving spacer distance distributions to estimate spacer loss rates from data.

Contribution

It introduces a mathematical model for CRISPR spacer evolution with ordered insertions and losses, providing explicit spacer distance distributions for small sample sizes.

Findings

Derived spacer distance distributions for n=2 and n=3 arrays.

Provided a method to estimate spacer loss rate from array data.

Abstract

Today, the CRISPR (clustered regularly interspaced short palindromic repeats) region within bacterial and archaeal genomes is known to encode an adaptive immune system. We rely on previous results on the evolution of the CRISPR arrays, which led to the ordered independent loss model, introduced by Kupczok and Bollback (2013). When focusing on the spacers (between the repeats), new elements enter a CRISPR array at rate $\theta$ at the leader end of the array, while all spacers present are lost at rate $\rho$ along the phylogeny relating the sample. Within this model, we compute the distribution of distances of spacers which are present in all arrays in samples of size $n=2$ and $n=3$. We use these results to estimate the loss rate $\rho$ from spacer array data.

Equations

$$(U,s_1,s_2,\dots )$$

$$(s_1,\dots ,s_{i-1},s_{i+1},s_{i+2},\dots )$$

$$\mathbf{P}(X_t\in A|X_r=x)=p_{t-r}(x,A)\text{ for all }r\le t.$$

$$S_{t,i}:=(S_t{)}_i.$$

$$E_0^{\ell ,\text{mathbbm}K}$$

$$E_{n+1}^{\ell ,\text{mathbbm}K}$$

$$V_i$$

$$\mathbb{P}(A=a,B=b)$$

$$(V_1,W_1),(V_2-V_1,W_2-W_1),(V_3-V_2,W_3-W_2),\dots$$

$$(V_1,W_1)\sim (C_1+A_1,C_2+B_1),(V_i-V_{i-1},W_i-W_{i-1})\sim (A_i,B_i),i=2,3,\dots$$

$$\sum _{b=0}^{\infty }\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)x^b=\frac{1}{(1-x{)}^{a+1}},i=0,1,2,\dots$$

$$x=\frac{1-e^{-\rho T}}{2-e^{-\rho T}}=1-\frac{1}{2-e^{-\rho T}},$$

$$\text{mathbbm}P(A=a)$$

$$\text{mathbbm}P($$

$$\displaystyle=\bigg{(}\Big{(}\frac{\theta}{\rho}(1-e^{-\rho T})\Big{)}^{v_{1}+w_{1}}e^{-2\frac{\theta}{\rho}(1-e^{-\rho T})}\,\frac{e^{-\rho T}}{2-e^{-\rho T}}$$

$$\displaystyle\qquad\cdot\sum_{v=0}^{v_{1}}\sum_{w=0}^{w_{1}}\binom{v+w}{v}\frac{1}{(v_{1}-v)!(w_{1}-w)!}\frac{1}{(\frac{\theta}{\rho}(2-e^{-\rho T}))^{v+w}}\bigg{)}$$

$$\displaystyle\qquad\qquad\cdot\bigg{(}\prod_{i=2}^{m}\binom{\Delta v_{i}+\Delta w_{i}}{\Delta v_{i}}\bigg{)}\cdot\bigg{(}\frac{e^{-\rho T}}{2-e^{-\rho T}}\bigg{)}^{m-1}\cdot\bigg{(}\frac{1-e^{-\rho T}}{2-e^{-\rho T}}\bigg{)}^{\sum_{i=2}^{m}\Delta v_{i}+\Delta w_{i}}$$

$$\text{mathbbm}P(V_1=v_1,W_1=w_1)$$

$$\cdot \sum _{a=0}^{v_1}\sum _{b=0}^{w_1}\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)\frac{1}{(v_1-a)!(w_1-b)!}\frac{1}{(\frac{\theta }{\rho }(2-e^{-\rho T}){)}^{a+b}}.$$

$$\text{mathbbm}P($$

$$=\sum _{a=0}^{v_1}\sum _{b=0}^{w_1}\text{mathbbm}P(A=a,B=b)\cdot \text{mathbbm}P(C_1=v_1-a)\cdot \text{mathbbm}P(C_2=w_1-b)$$

$$=\sum _{v=0}^{v_1}\sum _{w=0}^{w_1}\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)\frac{e^{-\rho T}}{2-e^{-\rho T}}{x}^{a+b}\frac{z^{v_1+w_1-(a+b)}}{(v_1-a)!(w_1-b)!}{e}^{-2z}$$

$$\displaystyle=z^{v_{1}+w_{1}}\frac{e^{-2z}e^{-\rho T}}{2-e^{-\rho T}}\sum_{a=0}^{v_{1}}\sum_{b=0}^{w_{1}}\binom{a+b}{a}\frac{1}{(v_{1}-a)!(w_{1}-b)!}\Big{(}\frac{x}{z}\Big{)}^{a+b}.$$

$${\int }_0^T\theta e^{-\rho t}dt=\frac{\theta }{\rho }(1-e^{-\rho T}).$$

$$\{S_{1,C_1+i}:i=1,2,\dots \}\cup \{S_{2,C_2+i}:i=1,2,\dots \}\subseteq S_{1\wedge 2}.$$

$$(\{s\in S_1\},\{s\in S_2\}{)}_{s\in S_{1\wedge 2}}$$

$$\text{mathbbm}P(s\in S_1)=\text{mathbbm}P(s\in S_2)=e^{-\rho T}$$

$$p_1$$

$$p_2$$

$$p_3$$

$$p_4$$

$$\text{mathbbm}P$$

$$=\sum _{k=a+b+1}^{\infty }\text{mathbbm}P(K=k)\text{mathbbm}P(A=a,B=b|K=k)$$

$$=\sum _{k=a+b+1}^{\infty }{p}_4(1-p_4{)}^{k-1}(\genfrac{}{}{0ex}{}{k-1}{a+b}\left)\right(\genfrac{}{}{0ex}{}{a+b}{a}){\left(\frac{p_1}{1-p_4}\right)}^a{\left(\frac{p_2}{1-p_4}\right)}^b{\left(\frac{p_3}{1-p_4}\right)}^{k-1-(a+b)}$$

$$=\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)p_1^a{p}_2^b{p}_4\sum _{k=a+b+1}^{\infty }\left(\genfrac{}{}{0ex}{}{k-1}{a+b}\right)p_3^{k-1-(a+b)}$$

$$=\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)p_1^a{p}_2^b{p}_4\sum _{k=0}^{\infty }\left(\genfrac{}{}{0ex}{}{a+b+k}{a+b}\right)p_3^k$$

$$\displaystyle=\binom{a+b}{a}\frac{p_{4}}{1-p_{3}}\Big{(}\frac{p_{1}}{1-p_{3}}\Big{)}^{a+b}$$

$$\displaystyle=\binom{a+b}{a}\frac{e^{-2\rho T}}{1-(1-e^{-\rho T})^{2}}\Big{(}\frac{e^{-\rho T}(1-e^{-\rho T})}{1-(1-e^{-\rho T})^{2}}\Big{)}^{a+b}$$

$$\displaystyle=\binom{a+b}{a}\frac{e^{-\rho T}}{2-e^{-\rho T}}\Big{(}\frac{1-e^{-\rho T}}{2-e^{-\rho T}}\Big{)}^{a+b}$$

$$X_i$$

Full text

††journal: arXiv

1 Introduction

The CRISPR Cas system is a widespread microbial adaptive defense mechanism against viruses and plasmids (Marraffini, 2015; Rath et al., 2015), that likely originated in archaea and spread to bacteria via horizontal transfer (Makarova et al., 2011a). The Clustered Regulary Interspaced Short Palindromic Repeats (CRISPR) have already been described in 1987 by Ishino et al. (1987). Later it turned out that the unique sequences between these repeats, so called spacers, are of foreign origin (Bolotin et al., 2005) and serve as an immunological memory passed to the offspring. New spacers are acquired and inserted at the leader end of the array (Barrangou et al., 2007), such that the order of spacers represents the chronological infection history of the bacterial population. Together with CRISPR associated (cas) genes these spacers can provide resistance against phages and plasmids by targeting molecular scissors to the corresponding sequences in the invading DNA (Barrangou et al., 2007).

The most prominent cas gene is Cas9, which recently led to a revolution in genome engineering (Doudna and Charpentier, 2014; Hsu et al., 2014). With the CRISPR-Cas9 enzyme mechanism an uncomplicated and cheap technology to alter the genome of potentially any organism is now available. Due to the precise targeting via engineered spacer sequences this system is speeding up the pace of research and gives rise to applications with incredible impact and opportunities in a variety of fields (Doudna and Charpentier, 2014). As just one among many examples, the concept of gene drive (Burt, 2003) in combination with the precision of CRISPR-Cas9 may enable us to alter the genetics of entire populations (Esvelt et al., 2014; Oye et al., 2014).

Here we will focus on the evolution of natural CRISPR-Cas systems and their spacer arrays in microbial genomes. CRISPR systems have been classified into different types and subtypes, with different sets of accompanying cas genes (Makarova et al., 2011b, 2015). A single genome can contain different types of CRISPR and the rates at which new spacers are inserted and old spacers are lost vary between the systems (Horvath et al., 2008). This suggests that different types may have different evolutionary dynamics and functions beyond defense, e.g. in regulation of gene expression (Westra et al., 2014). As another example, Lopez-Sanchez et al. (2012) suggest that CRISPR may control the diversity of mobile genetic elements in Streptococcus agalactiae.

So far we just got a glimpse of the ecological and evolutionary impact of CRISPR cas systems. In particular, the benefit of possessing a CRISPR system and the parameter regime where they are maintained (Levin, 2010; Weinberger et al., 2012), as well as the coevolutionary dynamics of bacteria containing CRISPR loci and phages (Koskella and Brockhurst, 2014; Han and Deem, 2017) have been considered. However many ecological and evolutionary aspects of CRISPR cas systems, as the frequent horizontal transfer of the whole system, are still not understood (Rath et al., 2015). Not only the evolution of the whole CRISPR system but the evolving and adapting spacer array itself is of interest. The spacer array represents snippets of previous phage/plasmid exposure that can help to disentangle the interplay between bacterial and viral populations (Childs et al., 2014; Sun et al., 2015). Since resistance is inherited by the offspring via the acquired spacer sequences, at least some CRISPR systems blur the distinction between Darwinian and Lamarckian modes of evolution (Koonin and Wolf, 2016). Modeling the evolution of spacer arrays will help to identify differences between CRISPR types and interpret the observed pattern of spacer insertion and deletion.

In 2013, Kupczok and Bollback introduced probabilistic models for the evolution of CRISPR spacer arrays. They concluded that a model with ordered spacer insertions and independent losses best describes the dynamics of spacer array evolution. In this model unique new spacers are inserted at the leader end and each spacer gets lost independently at a constant rate; see also Definition 2.1 below. Later, the model has been extended to conclude that there is no evidence for frequent recombination within the spacer arrays (Kupczok et al., 2015). In Kupczok and Bollback (2013) the constructed estimators for spacer insertion and deletion rates assume that no phylogenetic information is available. In contrast, in this paper we assume that the genealogy is known or has been reconstructed adequately, e.g. based on the cas genes in front of the spacer array. Given such a genealogy we look at distances between equal spacers, i.e. spacers that appear in more than one array. Kupczok and Bollback (2013) also considered an unordered independent loss model, where the order of spacers is irrelevant. In pangenome analysis (Mira et al., 2010; Vernikos et al., 2014), this model is known as the infinitely many genes model (Baumdicker et al., 2010) and has led to methods to jointly estimate gene gain and loss rates based on the frequency of genes in the sample (Baumdicker et al., 2012). These methods can directly be applied in our setting to jointly infer the rates of spacer insertion and deletion from (unordered) spacer frequencies. Here we compute the distribution of distances of (ordered) spacers, which are present in all arrays in samples of size $n=2$ and $n=3$. We show that including the order of spacer arrays by looking at equal spacer distances in a sample of arrays allows to decouple estimation of spacer insertion and spacer loss rate.

The paper is organized as follows: In Section 2, we introduce the ordered independent loss model. In Section 3, we compute the distribution of equal spacer distances of samples of size $n=2$ and $n=3$. As a by-product, we find sufficient statistics useful for the estimation of the loss rate. These are employed in Section 4 where we give maximum likelihood estimators and perform simulations in order to show their accuracy.

2 Model

The following model for the evolution of spacers in the CRISPR-system is based on work by Kupczok and Bollback (2013). Since this model was called the ordered model with independent losses, we follow this terminology here.

Definition 2.1** (The ordered independent loss model along a single line).**

Let $\mathcal{S}=(S_{t})_{t\geq 0}$ be a Markov jump process with state space $E=[0,1]^{\mathbb{N}}$ and the following dynamics: If $S(t)=s=(s_{1},s_{2},\dots)$, it jumps to

[TABLE]

The first type is called a gain-event (or insertion event), while the latter is called a loss-event (or deletion event). We refer to $S$ as the independent loss model with ordered gains, or as the ordered independent loss model, for short.

A simple result is the following, which is clear from the definition of $\mathcal{S}$.

Lemma 2.2** (Equilibrium of $S$).**

Let $(U_{1},U_{2},\dots)\in E$ be a sequence of independent, $U([0,1])$-distributed random variables. Then, the distribution $\pi$ of $(U_{1},U_{2},\dots)$ is an equilibrium of $\mathcal{S}$ from Definition 2.1.

Remark 2.3** (More equilibria).**

We note that $\pi$ is not the only equilibrium of the process. For example, let $N\sim\text{Poi}(\theta/\rho)$ and, conditioned on $N$, let $(V_{1},\dots,V_{N})$ be independent and $U([0,1])$-distributed, and $V_{N+1}=V_{N+2}=\cdots=1$. Then, the distribution of $(V_{1},V_{2},\dots)$ is stationary for $S$ as well.

Indeed, let $N_{t}:=\sup\{n:U_{1},\dots,U_{n}\notin S(0)\}$. Then, $\mathcal{N}=(N_{t})_{t\geq 0}$ is a death-immigration process with immigration rate $\theta$ and, if it is in state $n$, death rate $n\rho$. For this process, it is well-known that $N_{t}\xRightarrow{t\to\infty}N$; hence, necessarily, $V_{1},\dots,V_{N}$ are independent and uniform draws from $[0,1]$ for a stationary distribution. However, $V_{N+1},V_{N+2},\dots$ are states also found in $S(0)$, so if $S(0)=(1,1,\dots)$, then $S_{N_{t}+k}(t)=1$ for all $t$ and $k$ by construction.

In order to formulate our results, we require the ordered independent loss model not only along a single line, but also along an ultra-metric tree. For this, we introduce some notation.

Remark 2.4** (Ultrametric trees; tree-indexed processes).**

Recall that a tree $\mathbbm{T}$ with leaves $\mathbbm{L}\subseteq\mathbbm{T}$ is called ultra-metric if $d(\ell_{1},\ell_{3})\leq\max(d(\ell_{1},\ell_{2}),d(\ell_{2},\ell_{3}))$ for all $\ell_{1},\ell_{2},\ell_{3}\in\mathbbm{L}$. Here, $d(.,.)$ denotes the graph distance and $\ell\in\mathbbm{T}$ is a leaf if $\mathbbm{T}\setminus\{\ell\}$ only has a single connected component. Note that for such an ultra-metric tree there is a unique $r\in\mathbbm{T}$ (called the root) such that $d(r,\ell)$ does not depend on $\ell\in\mathbbm{L}$. In addition, there is an order on $\mathbbm{T}$ such that $s\leq t$ iff $s\in[r,t]$, where $[r,t]$ is the unique path from $r$ to $t$. Then, we also say that $s$ is ancestor of $t$. Using this order the most recent common ancestor $s\wedge t\in\mathbbm{T}$ of $s$ and $t$ is the largest element in $\mathbbm{T}$ which is ancestor of both, $s$ and $t$. 2. 2.

Usually, a Markov process $X=(X_{t})_{t\geq 0}$ has the property that $(X_{t})_{t\geq s}$ is independent of $(X_{t})_{t\leq s}$ conditional on $X_{s}$. Moreover, we call a Markov-process time-homogeneous if there is a family of transition kernels $(p_{s}(.,.))_{s\geq 0}$ such that

[TABLE]

If a Markov-process is piecewise constant, it is usually called a Markov jump process. If it is time-homogeneous, it follows from (2.1) that the waiting time to the next jump has an exponential distribution. Its parameter depends on the current state and is usually referred to as the rate of the exponential distribution. We call a tree-indexed process $X=(X_{t})_{t\in\mathbbm{T}}$ time-homogeneous Markov if (2.1) holds for all $r,t\in\mathbbm{T}$ with $r\leq t$, where $t-r:=d(r,t)$, and, conditional on the value of the process at a node, the processes in the two descending subtrees are independent. For more work on Markov processes indexed by trees, see e.g. Benjamini and Peres (1994).

Definition 2.5** (The ordered independent loss model along an ultrametric tree).**

Let $\mathbbm{T}$ be an ultra-metric tree with root $r\in\mathbbm{T}$. The $\mathbbm{T}$-indexed time-homogeneous Markov process $\mathcal{S}=(S_{t})_{t\in\mathbbm{T}}$ with $S_{r}\sim\pi$ (recall from Lemma 2.2) and transition kernels given through gain and loss events as in Definition 2.1, is denoted the ordered independent loss model (along $\mathbbm{T}$).

Remark 2.6**.**

Note that it is straightforward to formulate the presented results for non-ultrametric trees. However, since the formulas simplify significantly for ultrametric trees, we assume that $\mathbb{T}$ is ultrametric.

In the sequel, we fix the ultra-metric tree $\mathbbm{T}$, its set of leaves $\mathbbm{L}$ and the process $\mathcal{S}$ from Definition 2.5. Recall that $S_{t}\in[0,1]^{\mathbbm{N}}$ and we will use the shorthand notation

[TABLE]

Moreover, we will identify the vector $S_{t}$ with the set of its entries, i.e. $S_{t}=\{S_{t,i}:i=1,2,\dots\}$. Note that all entries of $S_{t}$ are different almost surely.

Definition 2.7** (Equal spacers).**

For $\mathbbm{K}\subseteq\mathbbm{L}$ and $\ell\in\mathbbm{K}$, define recursively

[TABLE]

Here, $E^{\ell,\mathbbm{K}}_{n}$ is the spacer position in $\ell$ of the $n$th spacer which is also contained in all $S_{k},k\in\mathbbm{K}$, but in none of $S_{k^{\prime}},k^{\prime}\in\mathbbm{L}\setminus\mathbbm{K}$.

An illustration of the process from Definition 2.5 and the corresponding equal spacers are shown in Figure 2.1.

3 Results

3.1 Trees with two leaves $|\mathbbm{L}|=2$

In the special case that $\mathbbm{L}$ consists of only two points, we denote these leaves by $\mathfrak{1}$ and $\mathfrak{2}$. In addition, we set $d(\mathfrak{1},\mathfrak{2})=2T$ for some $T\geq 0$ and define the following random variables:

[TABLE]

i.e. $V_{i}$ is the $i$th element of $S_{\mathfrak{1}}$ which is also contained in $S_{\mathfrak{2}}$ and $W_{i}$ is the $i$th element of $S_{\mathfrak{2}}$ which is also contained in $S_{\mathfrak{1}}$.

Theorem 1** (Distribution of equal spacer sequence in two leaves).**

Let $(A,B),(A_{1},B_{1}),(A_{2},B_{2}),\dots$ be iid pairs of random variables with joint distribution

[TABLE]

In addition, let $C_{\mathfrak{1}},C_{\mathfrak{2}}$ be iid Poisson distributed r.v. with parameter $\frac{\theta}{\rho}(1-e^{-\rho T})$. Then,

[TABLE]

are independent with

[TABLE]

Remark 3.1** **($A$ and $B$ are geometrically

distributed).

In the sequel, we will use the identity

[TABLE]

for $|x|<1$ on several occasions. It can easily be proven by induction. 2. 2.

For the distribution of $(A,B)$ from (3.1), we note that both, $A$ and $B$ are geometrically distributed with parameter $e^{-\rho T}$. (We come back to this observation in Remark 3.3.)

Indeed: For

[TABLE]

we have that

[TABLE]

So, we have shown that $A$ has the desired distribution. Symmetry then gives the same for $B$.

Before we come to the proof of Theorem 1, we state the sampling formula which results from Theorem 1.

Corollary 3.2** **(Sampling formula for equal spacer sequence in two

leaves).

The joint distribution of $(V_{1},V_{2},\dots,V_{m},W_{1},W_{2},\dots,W_{m})$ is given by

[TABLE]

for $m=1,2,\dots$, where $\Delta v_{i}:=v_{i}-v_{i-1},\Delta w_{i}:=w_{i}-w_{i-1},i=2,3,\dots$

Proof.

By the independence from Theorem 1 and the distribution given in (3.1), the only thing which remains to be proven is that

[TABLE]

From Theorem 1, we know that $(V_{1},W_{1})\sim(C_{\mathfrak{1}}+A,C_{\mathfrak{2}}+B)$, where the distribution of $(A,B)$ is as in (3.1) and $C_{\mathfrak{1}}$ and $C_{\mathfrak{2}}$ are independent and Poisson distributed with mean $z:=\frac{\theta}{\rho}(1-e^{-\rho T})$. Hence, with $x$ from (3.3),

[TABLE]

Plugging in $x$ and $z$ gives the result. ∎

Proof of Theorem 1.

For the leaves $\mathfrak{1}$ and $\mathfrak{2}$, recall that $\mathfrak{1}\wedge\mathfrak{2}$ denotes their MRCA (most recent common ancestor) and $d(\mathfrak{1},\mathfrak{1}\wedge\mathfrak{2})=T$. Let $C_{\mathfrak{1}}(C_{\mathfrak{2}})$ be the number of gain-events between $\mathfrak{1}\wedge\mathfrak{2}$ and $\mathfrak{1}$ ($\mathfrak{2}$), which don’t get lost until $\mathfrak{1}$ (until $\mathfrak{2}$). Then, by construction, $C_{\mathfrak{1}}$ and $C_{\mathfrak{2}}$ are independent (since they depend on independent gain events) and Poisson distributed with mean

[TABLE]

For $(S_{\mathfrak{1},C_{\mathfrak{1}}+i})_{i=1,2,\dots}$ and $(S_{\mathfrak{2},C_{\mathfrak{2}}+i})_{i=1,2,\dots}$, i.e. the spacers after the just mentioned gain-events, we note that, by construction,

[TABLE]

Moreover, we have that (by independence of loss-events at all positions) the events

[TABLE]

are independent with

[TABLE]

for all $s\in S_{\mathfrak{1}\wedge\mathfrak{2}}$. Equivalently, we find that independently for all $s\in S_{\mathfrak{1}\wedge\mathfrak{2}}$,

[TABLE]

Moving along $S_{\mathfrak{1}\wedge\mathfrak{2}}$, toss a four-sided die, numbered $1,\dots,4$. If it comes up 1, keep the spacer in $\mathfrak{1}$ but not in $\mathfrak{2}$; if it comes up 2, keep it in $\mathfrak{2}$ but not in $\mathfrak{1}$; if it comes up 3, neither keep it in $\mathfrak{1}$ nor in $\mathfrak{2}$; if it comes up 4, keep it in both, $\mathfrak{1}$ and $\mathfrak{2}$. Then, we need to compute the joint distribution of the number of die rolls $A$ with 1 and $B$ with 2 before the first 4 (which indicates an equal spacer). Letting $K$ be the first die roll with 4, we have that

[TABLE]

since $p_{1}=p_{2}$, where we have used (3.2) in the second to last equality. Before the first equal spacer, we have for $(V_{1},W_{1})$ the sum of new and old spacers, while for $(V_{i},W_{i}),i\geq 1$, we only have old spacers. Since loss events of spacers in $\mathfrak{1}\wedge\mathfrak{2}$ are independent, the result follows. ∎

Remark 3.3** **($\Delta V_{i}$ and $\Delta W_{i}$ are geometrically

distributed).

We have seen in the proof of Theorem 1 that the distribution of $(\Delta V_{i},\Delta W_{i})$ for $i=2,3,\dots$ can be obtained by rolling a die with probabilities $p_{1},\dots,p_{4}$ given by (3.4) and counting the number of occurrences of 1 and 2 before the first 4. For the marginal distribution of $\Delta V_{i}$, we are asking for the number of occurrences of 1 before the first 4. Clearly, this is geometrically distributed with success probability $p_{4}/(p_{1}+p_{4})=e^{-\rho T}$. See also the formal calculation in Remark 3.1.

3.2 Trees with three leaves $|\mathbbm{L}|=3$

In the special case that $\mathbbm{L}$ consists of three points, denoted $\mathfrak{1},\mathfrak{2}$ and $\mathfrak{3}$, we define the following random variables:

[TABLE]

i.e. $X_{i}$ is the $i$th element of $S_{\mathfrak{1}}$ which is also contained in both, $S_{\mathfrak{2}}$ and $S_{\mathfrak{3}}$, $Y_{i}$ is the $i$th element of $S_{\mathfrak{2}}$ which is also contained in $S_{\mathfrak{1}}$ and $S_{\mathfrak{3}}$, and $Z_{i}$ is the $i$th element of $S_{\mathfrak{3}}$ which is also contained in $S_{\mathfrak{1}}$ and $S_{\mathfrak{2}}$. In between $X_{i}$ and $X_{i+1}$, for example, we find spacers of three classes: those, which are only in $S_{\mathfrak{1}}$, i.e. in $S_{\mathfrak{1}}\setminus(S_{\mathfrak{2}}\cup S_{\mathfrak{3}})$, those which are shared with $S_{\mathfrak{2}}$, i.e. in $S_{\mathfrak{1}}\cap S_{\mathfrak{2}}\setminus S_{\mathfrak{3}}$, and those shared with $S_{\mathfrak{3}}$, i.e. in $S_{\mathfrak{1}}\cap S_{\mathfrak{3}}\setminus S_{\mathfrak{2}}$. Recalling the notation from Definition 2.7, we write for $\mathbbm{K}\subseteq\mathbbm{L}=\{\mathfrak{1},\mathfrak{2},\mathfrak{3}\}$

[TABLE]

where the right hand side does not depend on which $k\in\mathbbm{K}$ we choose, because the order is not altered by losses. In words, $F_{i}^{\mathbbm{K}}$ is the number of spacers within $S_{k}$, which are between the $i$th and $(i+1)$st spacer shared among $\mathfrak{1},\mathfrak{2},\mathfrak{3}$ and that appear in all elements of $\mathbbm{K}$, but in no element of $\mathbbm{L}\setminus\mathbbm{K}$. As an example, we can rewrite

[TABLE]

Note that $X_{i+1}-X_{i}=F_{i}^{\mathfrak{1}\setminus(\mathfrak{2},\mathfrak{3})}+F_{i}^{(\mathfrak{1}\cup\mathfrak{2})\setminus\mathfrak{3}}+F_{i}^{(\mathfrak{1}\cup\mathfrak{3})\setminus\mathfrak{2}}+1$ by construction, $i=0,1,2,\dots$

Since there is exactly one tree topology of a tree with three leaves, we will pick an ultra-metric tree, where $\mathfrak{1}$ and $\mathfrak{2}$ are closer relatives than $\mathfrak{1}$ and $\mathfrak{3}$ (and therefore also closer than $\mathfrak{2}$ and $\mathfrak{3}$). So, $T,T^{\prime}$ denote the distances from $\mathfrak{1}$ to $\mathfrak{3}$ and from $\mathfrak{1}$ to $\mathfrak{2}$, respectively, with $T\geq T^{\prime}$. See Figure 3.1 for an illustration of the tree topology.

Theorem 2** (Distribution of equal spacer sequence in three leaves).**

For $\mathbbm{L}=\{\mathfrak{1},\mathfrak{2},\mathfrak{3}\}$, let $(A^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq\mathbbm{L}},(A_{1}^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq\mathbbm{L}},(A_{2}^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq\mathbbm{L}},\dots$ be iid tuples of random variables with joint distribution

[TABLE]

with

[TABLE]

In addition, let

[TABLE]

be independent. Then,

[TABLE]

are independent with

[TABLE]

Remark 3.4** (Old and new spacers).**

In the proof of the Theorem, we will divide the spacers into old spacers, which were already present in $\mathfrak{1}\wedge\mathfrak{2}\wedge\mathfrak{3}$, the most recent common ancestor of $\mathfrak{1},\mathfrak{2}$ and $\mathfrak{3}$, and new spacers, which were gained after $\mathfrak{1}\wedge\mathfrak{2}\wedge\mathfrak{3}$. As in the proof of Theorem 1, we will argue that every old spacer has a chance to be kept until $\mathfrak{1},\mathfrak{2}$ and $\mathfrak{3}$. Here, we have to keep in mind that losses in $\mathfrak{1}$ and $\mathfrak{2}$ are not independent due to the chosen tree topology. 2. 2.

Although we formulate Theorem 2 analogously to Theorem 1 for two leaves, there are some conceptual differences. In particular, in the case of two leaves, $\mathfrak{1}$ and $\mathfrak{2}$, it is clear that only after the first equal spacers we can be sure that spacers are old in the sense that they are also contained in $\mathfrak{1}\wedge\mathfrak{2}$. In the case of three leaves, $\mathfrak{1},\mathfrak{2}$ and $\mathfrak{3}$, and the tree-topology from above, we know that the number of new spacers shared between $\mathfrak{1},\mathfrak{3}$, and between $\mathfrak{2},\mathfrak{3}$, is zero (since $C^{\{\mathfrak{1},\mathfrak{3}\}}=C^{\{\mathfrak{2},\mathfrak{3}\}}=0$). Hence, if the first spacer in $(S_{\mathfrak{1}}\cap S_{\mathfrak{3}})\cup(S_{\mathfrak{2}}\cap S_{\mathfrak{3}})$ arises, we know that subsequent spacers must be old. In particular, $(S_{\mathfrak{1},1},\dots,S_{\mathfrak{1},X_{1}}),(S_{\mathfrak{2},1},\dots,S_{\mathfrak{2},Y_{1}}),(S_{\mathfrak{3},1},\dots,S_{\mathfrak{3},Z_{1}})$ contain some information about which spacers are new and old, which cannot be gathered from $F^{\{\mathfrak{1}\}}_{1},F^{\{\mathfrak{2}\}}_{1},F^{\{\mathfrak{3}\}}_{1}$.

Proof of Theorem 2.

For the leaves $\mathfrak{1}$ and $\mathfrak{2}$, we denote their MRCA by $\mathfrak{1}\wedge\mathfrak{2}$ and $d(\mathfrak{1},\mathfrak{1}\wedge\mathfrak{2})=T^{\prime}$.

Let $C_{\mathfrak{1}}(C_{\mathfrak{2}})$ be the number of gain-events between $\mathfrak{1}\wedge\mathfrak{2}$ and $\mathfrak{1}$ ($\mathfrak{2}$), which don’t get lost until $\mathfrak{1}$ (until $\mathfrak{2}$). Then, by construction, $C_{\mathfrak{1}}$ and $C_{\mathfrak{2}}$ are independent (since they depend on independent gain events) and Poisson distributed with mean

[TABLE]

For $(S_{\mathfrak{1},C_{\mathfrak{1}}+i})_{i=1,2,\dots}$ and $(S_{\mathfrak{2},C_{\mathfrak{2}}+i})_{i=1,2,\dots}$, i.e. the spacers after the just mentioned gain-events, we note that, by construction,

[TABLE]

We take the same route as in the proof of Theorem 1. First, consider the new spacers, i.e. spacers gained after $\mathfrak{1}\wedge\mathfrak{2}\wedge\mathfrak{3}$. Note that there are no new spacers in $F_{0}^{\{\mathfrak{1},\mathfrak{3}\}},F_{0}^{\{\mathfrak{2},\mathfrak{3}\}}$ (hence we set $C^{\{\mathfrak{1},\mathfrak{3}\}}=C^{\{\mathfrak{2},\mathfrak{3}\}}=0$). Since gain-events follow a Poisson process, the number of new spacers in $(F_{0}^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq L}$ are independent and have a Poisson distribution which is uniquely determined by their mean. As an example, consider $F_{0}^{\{\mathfrak{1}\}}$, here, we have to take into account spacers gained between $\mathfrak{1}\wedge\mathfrak{2}\wedge\mathfrak{3}$ and $\mathfrak{1}\wedge\mathfrak{2}$, which are kept in $\mathfrak{1}$ but lost in $\mathfrak{2}$, and spacers gained between $\mathfrak{1}\wedge\mathfrak{2}$ and $\mathfrak{1}$, which are not lost until $\mathfrak{1}$. We compute

[TABLE]

for the mean of $F_{0}^{\{\mathfrak{1}\}}$, which corresponds to $C^{\{\mathfrak{1}\}}$ above. The rates of all other classes of new spacers are obtained accordingly.

Moving along the spacers in $S_{\mathfrak{1}\wedge\mathfrak{2}\wedge\mathfrak{3}}$, we now distinguish the new spacers into eight cases (spacers kept only in $\mathfrak{1}$, only in $\mathfrak{2}$, only in $\mathfrak{3}$, only in $\mathfrak{1}$ and $\mathfrak{2}$, only in $\mathfrak{1}$ and $\mathfrak{3}$, only in $\mathfrak{2}$ and $\mathfrak{3}$, kept in none and kept in all). We count the number of die rolls with a $1,2,\dots,6$ until the first $8$ appears. The probabilities for the die are defined analogous to (3.4) as

[TABLE]

Then we can write (3.5) as (note that $r=e^{\rho T}(1-q_{7})$)

[TABLE]

Again, before the first equal spacer, we split $(F_{0}^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq\mathbbm{L}}$ into the sum of new and old spacers, while for $(F_{i}^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq\mathbbm{L}},i\geq 1$, we only have old spacers. Since loss events of spacers in $\mathfrak{1}\wedge\mathfrak{2}\wedge\mathfrak{3}$ are independent, the result follows. ∎

As in Corollary 3.2, it is straight-forward to translate the last Theorem into a sampling formula, i.e. a formula for the distribution of the family of random variables

[TABLE]

3.3 Trees with $n$ leaves $|\mathbbm{L}|=n$

For the general case with $n$ leaves the distribution of equal spacer sequences is based on a recursion along the tree $\mathbb{T}$, which we introduce in Definition 3.5. Throughout this section, we fix the tree $\mathbb{T}$.

Definition 3.5** (Survival function).**

Let $\mathbb{T}$ be a binary tree with root $r$, a set of leaves $\mathbb{L}$ and internal vertices $\mathbb{V}$, including $r$. Recall that there is a semi-order on $\mathbb{T}$ such that $r$ is the smallest element, and $\mathbb{L}$ is the set of maximal elements. For an internal branch point $t\in\mathbb{V}$, we denote by $t_{+}$ and $t^{+}$ two points infinitesimally close to $t$, pointing in the two directions in $\mathbb{T}$ leading to bigger elements. The survival function $p:\mathbb{T}\rightarrow[0,1]$ is defined by

[TABLE]

using the initial condition $p(t)=1\text{ for }t\in\mathbb{L}.$ In addition, we set $q(t):=1-p(t),t\in\mathbb{T}$.

Roughly speaking, $p(t)$ is the probability that a spacer present at $t\in\mathbb{T}$ is present in at least one of the leaves of the subtree starting in $t$. Analogously, $q(t)$ is the probability that a spacer present in $t$ is lost in the whole subtree starting in $t$. This is formalized in the next Proposition.

Proposition 3.6** **(probability to loose a spacer along all paths from

$t$ to the leaves).

Let $\mathcal{S}$ be defined as in Definition 2.5. Then, for a spacer $s\in S_{t}$ the probability to be absent in all $l\in\mathbb{L}$, i.e. to get lost along each path from $t$ to $\mathbb{L}$, is

[TABLE]

Proof.

The result follows directly from the definition of $p$. Note, however, that if $p(t_{+})$ and $p(t^{+})$ are the probabilities to lose a spacer at least once in the subtree of $t_{+}$ and $t^{+}$, the probability to not lose it at all on the subtree of $t$ is $(1-p(t_{+}))(1-p(t^{+}))$. This explains the form of $p(t)$ for $t\in\mathbb{V}$. See also Proposition 5.2 in Baumdicker et al. (2010). ∎

For $n>3$ the tree topology may differ between trees of the same size. It is hence more involved to compute the distribution of spacer sequences for a given tree $\mathbbm{T}$. To formulate the result for the general distribution of equal spacer sequences in $n$ leaves we need further notation.

Definition 3.7**.**

For a tree $\mathbb{T}$ with root $r$, $n-2$ further internal vertices $v\in\mathbb{V}$, and a set of $n$ leaves, $l\in\mathbb{L}=\{\mathfrak{1},...,\mathfrak{n}\}$, we define for any $\mathbb{K}\subseteq\mathbb{L}$ (we use the semi-order $\leq$ also set-valued, meaning that $t\leq\mathbb{K}$ means $t\leq l$ for all $l\in\mathbb{K}$)

[TABLE]

Moreover, let $\mathbb{S}\subseteq\mathbb{T}$ be a subtree. We define its length by $\lambda(\mathbb{S})$. Note that $\mathbb{T}\setminus\mathbb{S}$ is a forest, which again consists of trees. So, we identify $\mathbb{T}\setminus\mathbb{S}$ with this set of trees. For such a tree $\mathbb{S}^{\prime}\in\mathbb{T}\setminus\mathbb{S}$, we set

[TABLE]

which – from Proposition 3.6 – is the probability that a spacer present at the root of $\mathbb{S}^{\prime}$ is absent in all leaves of $\mathbb{S}^{\prime}$.

Theorem 3** (The distribution of equal spacer sequences in $n$ leaves).**

For $\mathbbm{L}=\{\mathfrak{1},\mathfrak{2},\dots,\mathfrak{n}\}$, let $\mathcal{M}:=\mathcal{P}(\mathbb{L})\setminus\{\mathbb{L}\}$ be the set of true subsets of $\mathbb{L}$. Let $(A^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq\mathbbm{L}},(A_{1}^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq\mathbbm{L}},(A_{2}^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq\mathbbm{L}},\dots$ be iid tuples of random variables with joint distribution

[TABLE]

with

[TABLE]

the probability that a spacer gained at $v\in\mathbb{V}$ is present in $\mathbb{K}$ but absent in $\mathbb{L}\setminus\mathbb{K}$. In addition, let

[TABLE]

be independent. Then,

[TABLE]

are independent with

[TABLE]

Proof.

We use the same reasoning as in the proof of Theorem 2. Moving along the old spacers in $S_{r}$ that are present in the root of $\mathbbm{T}$, we have to distinguish $2^{n}-1$ cases, one for each $\mathbbm{K}\in\mathcal{M}$. The number $a^{\mathbbm{K}}$ corresponds to the number of spacers that have been lost in all $l\in\mathbbm{L}\setminus\mathbbm{K}$ and kept otherwise. The probability for a single spacer $s\in S_{r}$ to end up in the correct subset of $\mathbbm{L}$ is given by

[TABLE]

Just as in Theorem 2 we are not interested in the spacers $s\notin\bigcup\limits_{l\in\mathbbm{L}}S_{l}$ that have been lost in all leaves and

[TABLE]

Finally we wait until the corresponding $2^{n}-1$-sided die falls on the side for an equal spacer with probability

[TABLE]

To set everything together we count all possible combinations that end up in $(A^{\mathbbm{K}}=a^{\mathbbm{K}})_{\mathbbm{K}\in\mathcal{M}}$. This gives the multinomial coefficent in equation (3.7) and we are done.

For new spacers that have been gained along $\mathbbm{T}$ note that the number of new spacers found in $\mathbbm{K}$ and $\mathbbm{K}^{\prime}\neq\mathbbm{K}$ are independent Poisson random variables. Hence, we can first calculate the expected number of new spacers in $S_{v}$ that are gained between the two ancestral nodes $v_{\downarrow}$ and $v$ of $\mathbbm{K}$. The probability that a spacer arises between $v_{\downarrow}$ and $v$, and is still present in $v$, is

[TABLE]

To complete the result sum over all $v\in(r,v_{\mathbbm{K}}]\cap\mathbb{V}$ and multiply the summand with the probabilty

[TABLE]

∎

4 Estimating parameters

Our main results, Theorems 1 and 2, are directly applicable for maximum likelihood estimation of gain and loss parameters. However, some care must be taken in order to produce reliable results. In practice, several difficulties arise for a direct application:

There is only a finite number of spacers:

In our model, we assume that the observed equal spacers $i=1,...,m$ are a finite subset of an infinite number of joint spacers. Using (simulated or real) data, the number of spacers is limited. In particular, under the ordered independent loss model we can compute the likelihood restricted to the first $m$ equal spacers common to everyone. In contrast, in reality the number of equal spacers is a poisson distributed random variable $M$. The model could account for this fact by assuming a finite number of spacers in the ancestral species, but in this case we would lose the regenerative property. 2. 2.

The phylogeny must be given in advance:

While Kupczok and Bollback (2013) estimate the phylogeny (precisely, the time to the MRCA in a sample of size two) together with the spacer insertion/deletion ratio, we assume that the true phylogeny of the spacer arrays is known. In practice this phylogeny needs to be inferred from data, e.g. from a multiple alignment of the cas-genes or of conserved genomic regions. 3. 3.

There are more than three leaves:

Theorems 1 and 2 give likelihoods only for trees with two and three leaves, respectively. If one wants to study larger samples, the corresponding full likelihood is computationally expensive. In order to directly apply our results, we will only use sample sizes of at most three in the sequel.

Often, insertion and loss rates can only be estimated relative to each other, e.g. by estimating the ratio of both rates, see e.g. Baumdicker et al. (2012). In contrast, for the ordered CRISPR spacer arrays, we can estimate the spacer loss rate $\rho$ without any information about the insertion rate $\theta$, if we restrict ourselves to spacers older than the first common spacer; see Theorem 1 and Theorem 2. Having estimated $\rho$, the gene gain rate can e.g. be estimated by observing that the average number of spacers is $\theta/\rho$ in an equilibrium situation; see Remark 2.3.

In the following, we focus on the first step of estimating $\rho$ from the spacer arrays $S_{\mathfrak{1}},S_{\mathfrak{2}}$ in the case of a sample of size 2 and from $S_{\mathfrak{1}},S_{\mathfrak{2}},S_{\mathfrak{3}}$ in the case of a sample of size 3.

4.1 Estimating $\rho$ based on samples of size two

Turning to the statistical problem of estimating $\rho$, we will use the notion of an exponential family. Recall that a family of distributions $\mathbbm{P}_{\rho}$ is called a $k$-parameter exponential family, if, for some $\mathbb{N}^{m}$-valued $X\sim\mathbbm{P}_{\rho}$, it has the form

[TABLE]

for some functions $h$, $c=(c_{1},...,c_{k}),t=(t_{1},...,t_{k})$ and $d$. Notice that $t$ is called the (minimal) sufficient statistics for $(\mathbbm{P}_{\rho})_{\rho}$, i.e. maximum likelihood inference about $\rho$ depends on $x$ only through $t(x)$.

Assume we have two spacer arrays. To estimate the spacer deletion rate using a maximum likelihood approach rate we define:

Definition 4.1** **(number of (un)equal spacers in samples of size

$n=2$).

Let $S_{\mathfrak{1}},S_{\mathfrak{2}}$ be given. Calculate

[TABLE]

the total number of equal spacers in $\mathfrak{1}$ and $\mathfrak{2}$ and

[TABLE]

with $V_{i}:=E_{i}^{\mathfrak{1},\{\mathfrak{1},\mathfrak{2}\}}\text{and }W_{i}:=E_{i}^{\mathfrak{2},\{\mathfrak{1},\mathfrak{2}\}}.$

Theorem 4** (Sufficient statistics and maximum likelihood estimator for $\rho$).**

a) Let $m\geq 2$. The joint distribution of $\Delta V_{i}:=V_{i}-V_{i-1},\Delta W_{i}:=W_{i}-W_{i-1}$, $i=2,\dots,m$ only depends on $\rho$ and is a single-parameter exponential family with

[TABLE]

b) Let $D$ and $M$ be fixed to $d,m\geq 0$ respectively. Then, the maximum likelihood estimator of $\rho$ is given by

[TABLE]

and can be computed explicitly:

[TABLE]

Proof.

a) Using Theorem 1, we write $\underline{\Delta v}:=(\Delta v_{2},\dots,\Delta v_{m})$ and $\underline{\Delta w}:=(\Delta w_{2},\dots,\Delta w_{m})$, such that

[TABLE]

with

[TABLE]

and $t$ from above.

b) Formula (4.1) follows directly from Corollary 3.2. Further, we compute for $p=e^{-\rho T}$

[TABLE]

Hence, the maximum likelihood estimator is

[TABLE]

∎

Remark 4.2** **(Assuming independence leads to the same estimator for

$n=2$).

Recall from Remark 3.1 that $V_{i}-V_{i-1}$ (as well as $W_{i}-W_{i-1}$) is geometrically distributed with success parameter $p/(2-p)$ for $p=e^{-\rho T}$. Moreover, $V_{M}-V_{1}$ (as well as $W_{M}-W_{1}$) has a negative binomial distribution with parameters $M-1$ and $p$. Assuming that $V_{M}-V_{1}$ and $W_{M}-W_{1}$ are independent (which is clearly false in our model), $D$ has a negative binomial distribution with parameters $2(M-1)$ and $p$. A simple calculation shows that the maximum likelihood estimator of $p$ is again given by $p^{\ast}$ from (4.2). Therefore, the maximum likelihood estimator for $\rho$ in the case of a sample of size 2 can also be obtained by assuming independence of $V_{M}-V_{1}$ and $W_{M}-W_{1}$.

4.2 Estimating $\rho$ based on samples of size three

For $n=3$, we again obtain an exponential family, since $((F_{i}^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq\mathbbm{L}})_{i=1,2,...}$ are independent and identically distributed according to (3.5).

Theorem 5** (Sufficient statistics for $\rho$).**

Let $m\geq 1$. The joint distribution of $((F_{i}^{\mathbbm{K}})_{\emptyset\subsetneq\mathbbm{K}\subsetneq\mathbbm{L}})_{i=1,...,m}$ is a 4-parameter exponential family with

[TABLE]

[TABLE]

[TABLE]

Analogous to the case of two spacer arrays we can use the results of Theorem 2 and 5 to build an estimator based on three spacer arrays.

Definition 4.3** **(loss rate estimation in samples of size

$n=3$).

Let $S_{\mathfrak{1}},S_{\mathfrak{2}},S_{\mathfrak{3}}$ be given, and denote the time to the MRCA from $\mathfrak{1}$ and $\mathfrak{3}$ ($\mathfrak{1}$ and $\mathfrak{2}$) by $T$ ($T^{\prime}$) with $T\geq T^{\prime}$. Calculate

[TABLE]

the total number of equal spacers in $\mathfrak{1},\mathfrak{2}$ and $\mathfrak{3}$ and

[TABLE]

then, set (compare with (3.5))

[TABLE]

Maximizing the right hand side of the last expression is hardly done explicitly. We therefore rely on numerical optimization in our applications.

4.3 Comparison of estimation methods

We used simulated data in order to compare our estimates based on

the exact likelihoods based on equal spacers from Definition 4.1 for $n=2$ and Definition 4.3 for $n=3$, (darkgray boxes in Figure 4.1) and 2. 2.

the software panicmage, which is based on Baumdicker et al. (2012) and uses the spacer frequency spectrum in order to simultaneously estimate $\theta$ and $\rho$ (lightgray boxes in Figure 4.1).

In the simulations, we consider gains and losses in equilibrium, i.e. the root of the tree has a Poisson number of spacers with mean $\theta/\rho$. For the trees relating the leaves, we take the coalescent as a null model, i.e. pairs of lines all coalesce independent from each other at constant rate $1$ independently of $\rho$; see e.g. Wakeley (2008). This has the advantage that we can use results from Baumdicker et al. (2012) in order to interpret our results. For example, we know that in a sample of size $n=2$ (and $n=3$), the number of spacers common to both arrays (previously called $M$) has expectation $\frac{\theta}{\rho(1+2\rho)}$ (and $\frac{\theta}{\rho(1+2\rho)(2+2\rho)}$). We simulated $1000$ replicates. For each replicate a new coalescent tree was drawn. Throughout all simulations we set $\theta=100\rho$, such that the number of spacers in an array is Poisson distributed with mean $100$.

As Figure 4.1 shows, both estimation procedures give accurate results (in the sense that the estimator $\hat{\rho}$ is close to $\rho$, or $\hat{\rho}/\rho\approx 1$ in many cases), but there are some differences. Since the equal spacer likelihood methods are based on extensions of geometric distributions, the resulting estimators are (slightly) biased. In contrast, the likelihood based on the spacer frequency spectrum can be calculated from independent Poisson random variables (Baumdicker et al., 2010) and the corresponding estimator is thus unbiased. Concerning the variance of the estimators, the equal spacer likelihood methods give a slightly larger variance of the estimators for small $\rho$ (for $n=2,\rho\leq 0.25$, Fig. 4.1(a); and for $n=3,\rho\leq 0.5$, Fig. 4.1(b)). This might be due to the fact that this estimation procedure only takes spacers after the first common spacer into account. However, after the first equal spacer, most spacers are equal for low $\rho$, such that there is less signal in the data, making estimation of $\rho$ more difficult. Note that for larger loss rates (for $n=2,\rho\geq 0.5$, Fig. 4.1(a); and for $n=3,\rho\geq 1$, Fig. 4.1(b)), the estimate from panicmage has a much higher variance. This might be due to the fact that panicmage cannot estimate $\theta$ and $\rho$ separately and the likelihood surface is flat in one direction in the $\theta/\rho$-plane.

5 Discussion

The ordered independent loss model for CRISPR spacer evolution paves the way for various statistical applications based on equal spacers. However, when comparing the estimates for the spacer deletion rate with the unordered model in small sample sizes and limited deletion rates, the extra information of the order of spacers does not seem to produce more reliable results; see Figure 4.1. Nevertheless, this extra information will prove to be useful in further applications. In particular, since spacers are ordered by time of their insertion in each array, we can see which spacers in the arrays were present in the nodes of the phylogeny. If a spacer is present in a specific node in the phylogeny, all subsequent spacers were inserted before and hence must also be present within this node. This information is only available in the ordered model and will help to understand phenomena such as non-constant insertion rates or horizontal transfer.

The ordered independent loss model, as presented above, assumes constant insertion and deletion rates. The separate estimation of spacer insertion and deletion rates allows to investigate each mechanism independently. Besides selective effects, the insertion of spacers is more strongly dependent on the environment of the population than the random deletion of spacers is. The insertion might hence be time-dependent (Hynes et al., 2016b). In contrast, it is more reasonable to assume that the deletion of a neutral spacer occurs at constant rate. To test this hypothesis in future statistical work, we therefore aim at methods based on equal spacer distances. This will allow us to infer deletion rate variation between different branches in the given phylogeny, different CRISPR-Cas systems or positions in the spacer array.

To apply the presented approach to genmic data, the CRISPR spacer arrays and a reconstructed phylogeny are necessary. In addition the equal spacers among all the CRISPR arrays need to be identified. There are several CRISPR prediction tools that can be used to scan for spacer arrays in genomic data (Bland et al., 2007; Edgar, 2007; Grissa et al., 2007a). For archaeal and bacterial reference genomes predicted spacers are available online (Grissa et al., 2007b). While a typical array has less than 20 spacers, there are also arrays with several hundred spacers.

The ordered independent loss model will be useful to study the evolutionary ecology of CRISPR-Cas systems. For an extensive review of the evolution and ecology of CRISPR we refer the reader to Westra et al. (2016). Here we will highlight the role of spacer deletions in CRISPR-enhanced gene drives and metagenomic studies concerning CRISPR. Gene drive is the dependency of the offspring distribution of an individual on its genetic background. As such, this concept in concert with CRISPR can be used for the artificial alterations of natural populations (Burt, 2003; Hynes et al., 2016b). However, Unckless et al. (2017) argue that in natural populations resistance to a gene drive will likely occur. As a counter measure, different improvements of the CRISPR gene drive have been proposed (Esvelt et al., 2014). One suggestion is to use a CRISPR gene drive with more than one spacer, such that resistance can not evolve so easily. The effect of losing spacers in such a gene drive should be considered carefully before using any designed CRISPR system in natural populations. The same is true for the idea to ”vaccinate” bacteria possessing CRISPR via defective phage and plasmid sequences (Hynes et al., 2014, 2016a). Since any artificial CRISPR spacer insertion is transmitted to the offspring we should assess the fixation probabilities and extinction times of vaccinated bacterial strains before releasing them into natural environments.

To better understand the dynamics of CRISPR and the corresponding phages in natural populations metagenomic studies have been used. Using the sequenced CRISPR spacer arrays, the ancestral environment of the population can be inferred (Sun et al., 2015). In addition, different lineages can be distinguished, even for highly related bacterial strains (Kunin et al., 2008). Since the ordered independent loss model is based on a genealogical tree connecting the individuals, it cannot directly be used for metagenomic data. Nonetheless, the model provides a theoretical framework to have a more detailed look at the evolutionary pattern in the spacer array. In particular, the often observed conservation of old spacers at the trailer-end of the array is of interest (Weinberger et al., 2012). In the ordered independent loss model conservation of spacers correspond to small equal spacer distances, i.e. a low spacer deletion rate $\rho$. For large loss rates the number of spacers between equal spacers rises until finally all observed spacers differ. In a non-neutral setting with selection simulations of CRISPR arrays and the viral population show that the distances between equal spacers still increase for high spacer deletion rates. To which amount the CRISPR trailer-end is conserved due to selective forces compared to the effect of a low deletion rate remains a future challenge.

Acknowledgments

We thank Omer Alkhnbashi and Rolf Backofen for fruitful discussion and two anonymous referees for their constructive comments and hints. We thank an anonymous referee to encourage us to give our most general result, Theorem 3. This research was supported by the DFG through the priority program SPP1590. In particular, FB is funded in parts through Pf672/9-1.