Azzam Alwan

Rémi Cogranne

Pierre Beauseroy

Edith Grall-Maës

Laurent Debelle

Nicolas Belloy

Stéphanie Baud

Manuel Dauchez

Sébastien Almagro

An efficient statistical methodology for peptide 3D structure clustering

Keywords: Clustering, flexible peptide structure, Hierarchical classification, Outliers detection, Kernel PCA

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

According to the World Health Association [START_REF] Organization | Global health risks: mortality and burden of disease attributable to selected major risks[END_REF], cardiovascular diseases are the first cause of death: more people die every year due to cardiovascular problems than to any other health issue. One of the main actors of the vascular system integrity is Elastin, a polymer of a protein called tropoelastin. It confers resilience to the artery wall and its degradation, generated by aging and/or diseases (such as diabetes or atherosclerosis), is one of the symptoms/causes of such diseases. This degradation produces small fragments of protein called peptides. These peptides could be considered as circulating "molecular signals" that potentially activate response from organs and cells. The activity of a peptide is related to its three-dimensional structure and consequently to its physico-chemical properties. Therefore, for a better understanding of cardiovascular aging, it is crucial to identify among elastin peptides, the repetitive conformations, that may be related with their biological functions. These repetitive sequences, and their three-dimensional (3D) structures, are frequently observed in trajectories of molecular dynamics simulations. From examination of trajectories, it appears that these repetitive shapes may lead to detect a specific signature of a "key element" in molecular signaling. The aim of this first approach is to automatically identify, in this very large set of potential molecules, redundant molecular structures that will allow us to discriminate some specific 3D structures which could be related to the "active peptides" produced by the degradation of elastin. Based on peptides' characteristics, many methods have been proposed to identify repetitive conformations of peptides with the goal to understand their functionality. These studies could be divided into two main categories. The first one relies on the peptide sequence [START_REF] Bandyopadhyay | An efficient technique for superfamily classification of amino acid sequences: feature extraction, fuzzy clustering and prototype selection[END_REF]- [START_REF] Vipsita | An efficient technique for protein classification using feature extraction by artificial neural networks[END_REF]. Most of the methods from this category extract biological properties from peptide sequences and then use them in machine learning to determine some key functions. Actually, methods from this category are not sufficient to deduce peptide main functions because the functionality is related to both structures and thermodynamics. This especially explains why two peptides with different sequences may, or may not, have similar function depending on whether they may present the same conformation during a dynamic behavior. On the other hand, the second category, more closely related with the present paper, relies on the threedimensional geometric structure of peptides [START_REF] Wu | An efficient multiple protein structure comparison method and its application to structure clustering and outlier detection[END_REF]- [START_REF] Chazal | Persistencebased clustering in riemannian manifolds[END_REF]. The majority of methods in this category aim at solving the problem of comparing several protein structures by improving the methods of structures alignment e.g DALI and SSAP [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF], [START_REF] Orengo | Ssap: sequential structure alignment program for protein structure comparison[END_REF]. These two methods use the distance matrix to represent each protein structure. They have the major advantage of comparing different proteins with different sizes. Despite its accuracy and sensitivity, this method is still very computationally expensive because it explores all the possibilities of comparisons. Thus, as we will treat structures having the same size, these methods are not convenient for us as we cannot profit from their main advantages. We also have to keep in mind that elastin peptides are very specific (82% the tropoelastin sequences are composed of 5 amino-acids, with numerous repetitive patterns) and lead to the elasticity of the tissues, faraway from globular domains or transmembrane structures found in proteins with regular secondary structures (i.e. helices, strands or sheets and turns). Recently, several state-of-the-art works proposed to model the protein structure as a continuous parameterized curves rather than a sequence of discrete points [START_REF] Wu | An efficient multiple protein structure comparison method and its application to structure clustering and outlier detection[END_REF], [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF]. Although the efficiency of these methods, there is a risk of losing information related to the amino acid side chain due to the function used to interpolate between the discrete points to obtain the curves. Beside that, another approach was suggested and consists in using Self-Organizing Maps (SOM) method that tends to replace protein structures by neurons vectors [START_REF] Fraccalvieri | Conformational and functional analysis of molecular dynamics trajectories by selforganising maps[END_REF], then to submit it to a hierarchical clustering to extract the representative conformations. The only inconvenient of this method is that the SOM is very sensitive to the intuitive of its parameters which can significantly affect the results. Furthermore, a new approach is proposed in [START_REF] Chazal | Persistencebased clustering in riemannian manifolds[END_REF]. It is based on the data density function. The authors use the root-meansquared deviation (RMSD) as a distance metric for clustering. Additionally, the suggested method was applied on a protein with relatively well-understood dynamics, but even so, did not return the expected result for the number of clusters. Thus, as our objective is to observe the relation between the different clusters in the peptide dynamics and to avoid the alignment between structures as it is computationally expansive, we decided not to take this method into account. Indeed, all the clustering methods, that are mentioned previously, are only applied on the C α of the protein structures. Implicitly it leads to the conclusion that the side chains have no effect on the conformation variation, and that was proven wrong in the scope of our work. For those reasons, prior state-of-the-art methods are not adapted to tackle our problem, that is, to find the main patterns of peptides which are so flexible and are not regular as the protein treated in the previews methods.Thus the main contributions of the present paper are the following:

1. The proposed method includes a highly efficient algorithm to detect outliers conformations (that is, relies on a first step that aims at detecting outliers conformation) in order not to take into account the non-repetitive structures which may heavily impact the results of clustering. 2. An accurate description of the similarity between the main conformations can be easily obtained and illustrated after the automatic classification. 3. A novel approach combining two statistical methods is proposed to enhance the performance of the clustering method and to make it simpler and easier. The rest of the paper is organized as follows. The problem formulation is presented in Section 2. In Section 3, the proposed Methodology is described. Experimental results evaluating the efficiency of the proposed method are presented in Section 4. Finally, conclusions are drawn in Section 5.

II. PROBLEM FORMULATION

This section formally states the problem of identifying patterns that frequently occur among peptide structures before presenting the proposed method to address this problem. A peptide is essentially a set of related atoms in space. The three-dimensional position of the N atoms that compose a peptide will be referred to as the structure of the peptide and will be noted as S. Of course, such structure is a set of N three positions, or vectors : S = (a 1 , a 2 , . . . , a N )

T with a n ∈ R 3 the position of n-th atom. This paper studies sequence over time of peptide structures since, as for any molecule, the atom move relatively to each other under thermodynamical influence. Let us denote S the inspected sequence of T structures, that is S = (S 1 , . . . , S T ).

As described in the introduction, the goal of the present paper is to propose a method for automatically finding, without manual intervention, conformations that occur the most frequently within S a sequence of peptide structures. Formally speaking, two structures S t and S u have exactly the same conformation if and only if there exist a rotation, characterized by matrix R, and translation, characterized by vector ζ ∈ R 3 , such that

S u = S t R + 1 N ζ T (1)
with 1 N a vector of N row containing only ones. And consequently the conformation F can be considered as the set of all structures that have the same conformation:

F = S|S u = S t R + 1 N ζ T (2)
with S u the structure from which the conformation is defined, or any instance of structure in the set F. The main conformations of a peptide sequence S are the ones in which falls the higher number of structures from sequence S. However, it is obvious that because the moves of each atom is somewhat random, finding two structures with exactly the same conformation, as defined in Eq. ( 2), is very unlikely. Because we are interested in finding structures that are associated with the same function as the peptides, one can define a conformation in a wider sense as:

F = S|distance(S u , S t R + 1 N ζ T ) < ε (3)
where ε defines the upper bound on distance from structure S u that is allowed to consider a structure S as having the same conformation; ε is thus the maximal distance that affects the function of peptides. Based on the "practical" definition (3) of the conformation, it is obvious that any clustering algorithm should assign all the observations, or structures, from the inspected sequence S, with the conformation they fall within, if any. Structures that do not fall within any main conformation can be considered as outliers, those can typically occur when changing from one main conformation to another, also referred to as "transitional conformations". Let us denote L(t) the label of structure S t and {F 1 , F 2 , . . . , F K } the set of all K main conformations; the described assignment rule can be defined formally as follows:

     L(t) = k if S t ∈ F k L(t) = 0 if ∀k ∈ 1, . . . , K , S t / ∈ F k L(t) = 0 if S t ∈ F k and S t ∈ F l (4)
The Equations ( 1) -( 4) state the core problems of this paper, that is defining the K main conformations from a sequence of peptide structures S, and also allow to formalize the main difficulties addressed in this work. The first one is due to the alignment between the peptide structures, see Equations ( 1) -( 3), which is required to defined conformations independently from the rotation and translation of peptide structures. However, this is computationally expensive and is not always reliable in our purpose. Hence, the proposed method must avoid using structure alignments but must carefully take into account rotation and translation properly. The second difficulty is the unknown number of outliers within the sequence of structures S. Moreover, the definition of outliers (3) is itself fuzzy, this is due to the fact that the distance that defines a conformation is not known and may changes for various classes. The proposed approach must be able to avoid the effects of the structures that do not resemble to all conformations without knowing how to set the tolerance of similarity. The third difficulty is the unknown number of the main conformations in the sequence S. Indeed, several methods exist to perform clustering without any prior knowledge of the number of clusters [START_REF] Birant | St-dbscan: An algorithm for clustering spatialtemporal data[END_REF], [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]. However, it is not enough for us as we are interested in studying the similarities between the different main conformations. Thus, the proposed method must be flexible allowing thus the user to visualize the different clusters found, without additional computational costs, and to adjust the results manually when user needs to divide a conformation into two or more sub-classes III. METHODOLOGY This section presents the proposed methodology for peptides clustering that achieves the goals laid out in the previous section. In the first place, regarding the difficulty that is related to the rotation and translation effects of structures in space, it is proposed to represent each peptide structure as a distance matrix M, which characterizes the distance between atoms of structure [START_REF] Torda | Algorithms for clustering molecular dynamics configurations[END_REF]. Consequently, the rotation and translation effects are automatically overcome. Taking the same set of peptide structures defined in the previous part which is composed of N atoms and has T structures. The "distance matrix" M Si of a structure S i is defined as follows:

M Si (k, l) = a k -a l 2 , (5)
where v 2 stands for the Euclidean norm of vector v.

The matrix M Si is considered as an observation in the socalled "conformational space". The dimension of this space changes according to the number of atoms in structure. Thus, for the current peptide, the dimension of its space is equal to N (N -1) /2, due to the symmetry of M. In the second place, as it is noted previously, it exists structures that can be considered as outliers because of their unique conformations and can have effects on the clustering results. Hence an automatic outlier detection algorithm is proposed to detect these outliers. Moreover, this algorithm gives the possibility to the user to adjust the detection accuracy in order to ensure his clustering results. Finally, the next step is to use a clustering method on the remaining data to detect the principal conformations and to study the relations between them. Thus, an algorithm that does not need a prior knowledge of the clusters' number has been proposed. It is able to show the relations between the different structures and it is intelligent enough to adjust easily the number of clusters according to the data and the requirements.

In the following part, the method that is used to ignore the transitional conformations (outliers) will be described in Section III-A. Then in Section III-B, the clustering method that is used to classify the data will be explained.

A. Outliers detection

To detect properly the main conformations, the transitional conformations which can be found between them must be eliminated from the data. Recently, Kernel PCA has been investigated for outliers detection, and it has shown a high performance. It will be briefly described in the following subsection.

1) Kernel PCA: Kernel PCA is an extension of PCA to deal with the data non-linearity. Recently this method has been used for outliers detection. It demonstrates a competitive performance as compared to other outliers detection methods [START_REF] Schölkopf | Support vector method for novelty detection[END_REF]- [START_REF] Hoffmann | Kernel pca for novelty detection[END_REF]. Its concept is to map the data from the original space to a higher dimensional feature space F through a nonlinear transformation Φ, then the PCA will be applied in this space and a reconstruction error Recc will be calculated for each observation as an index of outliers detection. In the linear case, for a point x i ∈ R z , Recc is performed as follows :

Recc = x i -V q V T q x i 2 (6)
where V q is the matrix containing the q principal components, V q V T q x i is the estimated point after reconstruction, and q represents the number of the principal components. In F, Equation ( 6) is achieved as following

Recc = (Φ(x i ) • Φ(x i )) -(WΦ(x i ) • WΦ(x i )) (7)
where W is a matrix contains the q eigenvectors that are obtained from Kernel PCA. For more details about this treatment see [START_REF] Hoffmann | Kernel pca for novelty detection[END_REF]. Indeed, the principal components are affected directly by the main clusters that constitute the majority of the data. In other words, the observations that belong to these clusters are closer to the principal components than those which are outside of them. By consequence, the main observations have a low reconstruction error values and the outliers have the greatest reconstruction error values. To employ this method, two parameters have to be optimized; the number of eigenvectors q used for PCA and the kernel parameter (σ in the case of Gaussian kernel), K (x i , x j ) = exp -x i -x j 2 /2σ 2 . The choice of these parameters will be discussed in Section IV-B

B. Clustering Analysis

Hierarchical clustering is one of the popular clustering algorithms [START_REF] Gordon | A review of hierarchical classification[END_REF]. This method has been selected due to its flexibility and its accuracy. In this method it is not necessary to define the number of cluster beforehand. Since our goal is to study the data dispersion as well as to find out the main conformation, hierarchical clustering provides an accurate description of the data structures by using a dendrogram which represents similarity within the data on different levels. Moreover, by cutting the dendrogram at different levels, different clustering could be found with different numbers of clusters. To understand how it works, and for more details, one can refer to [START_REF] Day | Efficient algorithms for agglomerative hierarchical clustering methods[END_REF].

IV. EXPERIMENTS AND RESULTS

This section describes the data that are used in the experiments, discusses the selection of the proposed method's parameters (q, σ), evaluates the proposed method and finally shows the results obtained with the real data set of peptides.

A. Data sets

In this paper, two data sets have been used. The first, DB1, is a real data set. It consists of sequences of 3D related points. They represent the spatial configuration of the atoms composing the various peptides from elastin. Each peptide is composed of N atoms and 40 000 samples have been considered. Each sample correspond to a specific 3D position of all the atoms at sampling time with fixed to 200ns in molecular dynamics trajectory simulations. The second data set, DB2, has the same features of peptides in order to evaluate our approach. It is a simulated dataset that is formed of 9 atoms with each 1300 positions for 4 original conformations, and 350 observations for the transitional conformations. These 350 samples should be removed by the outlier detection step. The dynamics of atoms is assumed to be Gaussian noise as illustrated in Fig. 1.

B. Common core for all experiments

In this section, the selection of the proposed method's parameters is discussed. As mentioned previously, two main parameters have to be optimized, the number of eigenvectors, q, used in kernel PCA, and the Gaussian Kernel standard deviation σ. In fact, we are in unsupervised case. Hence, in order to evaluate the proposed approach and to optimize the parameters, one has to use methods that must be independent of the number of outliers and clusters. A list of these methods has been presented in [START_REF] Pimentel | A review of novelty detection[END_REF]. Since there is no information regarding the real data set, it is proposed to follow the method that consists to apply our algorithm on a simulated data (DB2), where the number of outliers and clusters is known, and study its capacity to detect the outliers. To evaluate detection performance the ROC (Receiver operating characteristic) curve has been used; this ROC curve shows the pobability of correct outlier detection against the probability of misclassifying one sample from main conformations as an outlier. However, to avoid having a curve for each setting, we reduce the evaluation to the Area Under the Curve (AUC) [START_REF] Bradley | The use of the area under the roc curve in the evaluation of machine learning algorithms[END_REF] which summarize the performance of a binary-detector. The effects of σ and q on the outliers detection are shown in Fig. 2a and Fig. 2b respectively. It can be seen in Fig. 2a that for a small value of σ , Kernel PCA has a good performance for almost all the q values used. When the value of σ increases, the AUC value is decreased. That is why a large value of σ is not recommended. Despite that, the optimal value of σ is still unknown. In addition, by observing the yellow and red curves in Fig. 2b, that correspond to small σ values, it can be seen that these curves are rather stable and q are very similar in terms of AUC values. Otherwise, when σ is far from its reasonable value which is estimated to be small, the results are not conclusive. Accordingly, it is not needed to assign a large value for q to have a good performance as kernel PCA could achieve it with a small value. Furthermore, the optimal value of q also depends on the data inertia rate, which represents the rate of information in the data after performing Kernel PCA. Fig. 3 presents q versus σ for four different inertia rates. In order to observe the behavior of Kernel PCA for different inertia rates, the coordinates of points of each curve in Fig. 3 are considered as parameters for Kernel PCA in order to calculate the detection and false alarm probabilities. Looking at the illustrations in Fig. 4, it is clear that the circled zone (where β has its highest values and υ has its lowest values) has a better performance than the other regions. By comparing that to the curves in Fig. 3, we see that this zone corresponds to the inflection points. Consequently, these illustrations help us to estimate the two optimal parameters. σ has a range of optimal values illustrated in Fig. 3. If it has a small value, we discretize the data, in the sense that each point becomes an eigenvector, resulting in a poor representation of the data variance. Furthermore, if the value of σ is large, the kernel PCA becomes meaningless, as explained in [START_REF] Hoffmann | Kernel pca for novelty detection[END_REF]. Meanwhile, increasing the value of q has no influence on the outliers detection. Hence, it is preferable to choose a low value for q, because with large values of q, the Kernel PCA treatment becomes time consuming.

AUC q =2 q=5 q =25 q =70 q =150 (

C. Real data

After performing the proposed algorithm on the simulated data (DB2) and finding the procedure to determine the optimal Fig. 4: Probability of detection and false alarm for various inertia rates. Taking the point coordinates of each curve in Fig. 3 as parameters for the kernel PCA algorithm. Each step in the abscissa corresponds to a σ value and its corresponding q value on the curves of Fig. 3.

parameters, we will apply our proposed approach on the target data (DB1) to see if it fulfills our expectations. Since the data are not labeled, it is not possible to recognize outliers from clusters. However, by looking at the available data of one of the peptides, we find 5 000 successive observations that can be classified into 3 clusters (We will use this amount of peptide structures only for clarity, otherwise we can apply our method on all the 320 000 structures). That can be concluded from the representation of the similarities between them. Knowing that these 5 000 observations can be classified into 3 clusters, it is suggested to apply the proposed classification method on them and to find out the optimal parameters. However, in this case, we cannot use the ROC curve to evaluate the results because the number of outliers is unknown. Hence, it is proposed to use a method that consists in exploiting the knowledge of the class number in these 5 000 observations, and to use the stability of the class center positions, after clustering, as an evaluation criterion. Then, to determine the optimal parameters (σ, q), we apply the following steps: first, we illustrate the curve that represents the number of eigenvectors q versus σ to have a range of reasonable values for σ, as shown in Fig. 3. Next, we take a value for σ from this range and then we use our method with it and with various values of q. In order to ignore the transitional forms, it is proposed to consider 20% of the data as outliers. Indeed, our main objective is to find the main forms and not to know how many outliers there are (it will be studied in a extension of this paper). So, even if we consider such a number of observations as outliers, it is not a problem for us in this approach, and it can be modified according to the requirements and the domain expert's demands. In our case, this percentage can be justified from a biological point of view. Due to space limitation, we will not illustrate the variation of the three centers positions of the three clusters in function of q. Briefly, after q = 12 the positions of the three centers become stable which means that we do not need more than 12 eigenvectors to perform our method. We applied this simulation for different values of σ from the values' range that is noted above and the same results are obtained. That is why, we adopt σ equals 5, and q equals 12 for our algorithm.

D. Results

In this paper, only a restricted number of cases have been presented. To show the performance of the proposed method, the Kernel PCA is applied on the 5 000 observations (From DB1) that are used to determine the two parameters, with σ = 5 and q = 12. Then the 1 000 observations that have the greatest values of Recc are considered as outliers, and they are eliminated from the data. Next, hierarchical clustering is employed on the 4 000 remaining observations to find out their main conformations. The dendrogram in Fig. 5 shows clearly that these latter observations can be divided into 3 clusters. We used CH (Calinski and Harabasz) as index for clustering validity [START_REF] Charrad | Nbclust: An r package for determining the relevant number of clusters in a data set[END_REF]. Fig. 6 shows the 3 main conformations present in these observations. These three structures are chosen to be the closest 3D structures to the three resulting clusters' centers. To ensure that this dendrogram reflects accurately the similarity between structures, we divided each main cluster to two subclusters, then we compared the difference between them in 3D dimension and in theoretical measure. Fig. 7 shows the

V. CONCLUSIONS

This paper studies the clustering of 3D structures of highly flexible elastin-derived peptides computed from numerical molecular dynamics simulations. It introduces a new strategy that combines different statistical methodologies to detect the main conformations of a peptide during a period of time. It proposes to use the distance matrix as a representation for the peptide's structure to avoid the effects of both translations and rotations over time. Moreover, it employs for the first time Kernel PCA to detect outliers conformations that do not resemble to any others. The removal of such outliers has been shown to be very efficient as it ensures a classification that is not impacted by those points that essentially lie between most frequently found conformations. Using Hierarchical Agglomerative Clustering method on the remaining data is eventually proposed to ensure a large flexibility in the classification of different conformations. Indeed by setting the threshold, one is free to choose between a coarse and a fine grain classification and is also able to visualize the difference between cluster centers to help tune the method according to the data and the requirements. Experimental results, performed using an existing database and simulated data, have shown that the proposed method is valid and efficient.

Fig. 1 :

 1 Fig. 1: Illustration of the 4 main conformations in blue, red, green and yellow. The transitions are shown in purple.

Fig. 2 :

 2 Fig. 2: AUC criterion in respect to : σ in (a) and q in (b)with various values of q and σ respectively.

Fig. 3 :

 3 Fig. 3: Number of eigenvectors q versus σ for four different inertia rates.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Dendrogram obtained after applying the hierarchical clustering on the 4 000 observations without Outliers.

(a) subcenter 4 (b) subcenter 5 Fig. 7 :

 457 Fig. 7: The structures which correspond to the centers of sub-clusters 4 and 5 in Fig. 5.

Fig. 8 :

 8 Fig. 8: Difference between the distance matrices M of the centers of th sub-clusters 4 and 5. Blue color means low values and red color high values


Francesca B Lopez

Antoine Fort

Luca Tadini

Aline V Probst

Marcus Mchale

James Friel

Peter Ryder

Frédéric Pontvianne

Paolo Pesaresi

Ronan Sulpice

Peter Mckeown

Galina Brychkova

Charles Spillane
email: charles.spillane@nuigalway.ie

Gene dosage compensation of rRNA transcript levels in Arabidopsis thaliana lines with reduced ribosomal gene copy number

Ribosomal gene reduction does not affect rRNA homeostasis but is correlated with loss of enrichment of heterochromatic marks, signs of genome instability

INTRODUCTION

Ribosomal RNA (rRNA) genes are one of the most abundant repetitive elements in eukaryotic genomes and are typically clustered as tandem repeats containing hundreds or thousands of rRNA gene copies, known as rDNA and usually associated with the nucleolus.

The 45S rDNA loci consist of sections encoding the 18S, 5.8S and 25S/28S rRNAs, which comprise the catalytic components of the mature ribosome, two Internal Transcribed Spacers (ITS), as well as two External Transcribed Spacers (ETS) (Sáez-Vásquez and Delseny, 2019).

The 45S arrays are associated with large chromosomal regions known as Nucleolar Organiser Regions (NORs): when stained with silver nitrate the NORs appear as dark areas on metaphase chromosomes and can be observed forming the nucleolus during interphase (McClintock, 1934).

To meet cellular demand for production of rRNA during cell proliferation, eukaryotes have a dedicated transcription system, RNA Polymerase I (Pol I) and highly duplicated 45S rDNA regions which provides hundreds or thousands of 45S rDNA copies to transcribe. The transcription of 45S rRNA can comprise in excess of 80% of total cellular transcription during proliferation. The epigenetic mechanisms which orchestrate the activation or silencing of rRNA genes have been widely documented in a range of model organisms (McStay and Grummt, 2008;Grummt and Pikaard, 2003). For instance, the Saccharomyces cerevisiae genome features ~150 rRNA copies on chromosome XII which are regulated by a mechanism which has likely evolved to stably maintain the number of rDNA copies in the genome (Kobayashi, 2006). Notably, silenced 45S rDNA copies are involved in maintenance of genome stability and cell senescence (Kobayashi, 2014), highlighting their functional role (s) in the maintenance of cellular homeostasis. It has been demonstrated that yeast cells with 80% reduction in rDNA copies suffer from increased DNA damage due to a reduced ability to repair Double Strand Breaks (DBS) (Ide et al., 2010). Similarly, in mammals a role for silent rDNA copies in the maintenance of genomic stability is now well established [START_REF] Stochaj | Nucleolar Organization and Functions in Health and Disease[END_REF], where the reduction of 45S rDNA Copy Number (CN) (often coupled with Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 amplification of 5S rDNA loci) has been associated with the onset of a range of cancers (Udugama et al., 2018;Xu et al., 2017).

Plants can exhibit extensive Copy Number Variation (CNV) at their 45S rDNA loci.

For instance, in inbred wild accessions of Arabidopsis thaliana sourced from across Sweden, the 45S rDNA CN ranged from 500 to 2500, with rDNA CN strongly correlating with genome size (Long et al., 2013). In Arabidopsis thaliana, the rDNA arrays are located on the acrocentric chromosomes 2 and 4 (Copenhaver and Pikaard, 1996), adjacent to the telomeric repeats. Remarkably, even closely-related accessions of Arabidopsis thaliana can display considerable 45S rDNA CNV, which has been predominantly attributed to variation in 45S rDNA repeats on chromosome 2 (NOR2) (Rabanal et al., 2017). The same study reported that 45S rDNA CNV can also be found within isogenic lines of the same accession, as well as in in Recombinant Inbred Lines (RILs). The underlying molecular basis for such variation appears to be due to crossing-over events, rather than the activity of homologous recombination (HR). In support of this, Sims et al., found that rDNA loci are shielded from HR components during meiosis, likely as a mechanism to avoid deleterious non-allelic interactions (Sims et al., 2019).

In multicellular eukaryotes, the transcriptional silencing of rRNA occurs during development, a phenomenon requiring the reorganization of chromatin and DNA methylation patterns across megabase tracts of chromosomes. The mechanisms to achieve multimegabase silencing which have evolved at cellular and organismal levels are diverse among different model organisms (Bersaglieri and Santoro, 2019). For instance, as mammalian cells exit a pluripotent state and begin differentiating, the silencing of rDNA copies is established de novo by the Nucleolar Remodelling Complex, NoRC (Santoro et al., 2002) through cytosine methylation of the rDNA promoter, accompanied by tight nucleosomal packaging of the rDNA coding sequence (associated with enrichments in silencing histone marks such as H3K9me2).

Nucleolar dominance was initially observed by Navashin as secondary constrictions in chromosomes of F1 progeny of interspecific hybrids, where the constriction only occurred on the chromosomes inherited from one of the parents (Navashin, 1934). Nucleolar dominance was later shown in Xenopus hybrids to be due to the uniparental expression of rDNA loci (Honjo and Reeder, 1973). In plants, many epigenetic mechanisms that regulate rRNA transcription in the context of nucleolar dominance have been studied (Chen and Pikaard, 1997;Lawrence et al., 2004), where the epigenetic mechanisms which maintain rDNA silencing in nucleolar dominance appear to also control silencing of rDNA during Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 development. Unlike in mammals, such transcriptional silencing of rDNA loci is not dependent on the NoRC, but instead achieved by the concerted action of a number of chromatin remodelling factors (particularly HISTONE DEACETYLASE6), both cytosine and histone methyltransferases, as well as methylcytosine binding domain proteins (MBD6 and MBD10) which together mediate the large-scale silencing of rRNA genes (Preuss et al., 2008;Probst et al., 2004;Tucker et al., 2010;Pontvianne et al 2012).

In the Arabidopsis thaliana accession Col-0 (which harbours ~375 45S rDNA copies per Nucleolar Organiser Region) NOR silencing has been characterized as a chromosomal position-dependent phenomenon, where NOR2 is developmentally silenced 10-20 days post germination, while NOR4 remains available for transcription throughout vegetative development, thus leaving about 50% of 45S rDNA copies competent for transcription (Pontvianne et al., 2010;Mohannath et al., 2016). As a result, in plants, as in animals, rRNA genes can be classified as active, inactive or silent depending on their chromatin organization (McKeown and Shaw, 2009). Ascribing any functional role of rDNA copy number variation in plants has remained elusive due to the lack of molecular tools to elicit a targeted reduction in rDNA copy number.

However, in Arabidopsis thaliana reduction of rDNA copy number has been reported in loss-of-function mutants of two of the large subunits of CHROMATIN ASSEMBLY FACTOR1 (CAF-1) (Mozgova et al., 2010). The CAF-1 complex is necessary for H3-H4 deposition and chromatin assembly following DNA replication: it is formed of three protein subunits, FASCIATA1 (FAS1), FASCIATA2 (FAS2) and MULTICOPY SUPPRESSOR OF IRA1 (MSI1). The fas1 fas2 double mutant displays progressive transgenerational shortening of telomeres on all chromosome arms, which is further associated with loss of 45S rRNA copy number on NOR2 and NOR4 (Pontvianne et al., 2013). The phenotypes of the fas1 fas2 double mutants are lost when WT alleles are reintroduced, whether at early or late generations (Pavlištová et al., 2016). Genetic complementation approaches of the fas1 fas2 mutant have generated a FAS1 FAS2 complemented line which features a loss of up to 80% of 45S rDNA copies, but reported to resemble the WT phenotype (Pavlištová et al., 2016).

Here, we use CRISPR-Cas9 induced Double Strand Breaks (DSBs) at 45S rDNA loci that can be used to generate plant lines with altered 45S rDNA copy number. To identify the lower limits of 45S rDNA copy number which still allow for viable plants, over six generations we generated lines with up to 93% reduction in 45S rDNA copy number. We have also investigated how rRNA transcript rates can be maintained despite such drastic reductions in 45S rDNA copy number and demonstrate that dosage compensation of rRNA Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 production appears to be achieved through chromatin re-organization at the rDNA loci, rather than through alterations of the transcription steady state. Using Nanopore genome sequencing, we further screened for genomic alterations in two independent 45s rDNA Low Copy Number (LCN) lines to determine whether genome stability could be compromised by loss of 45S rDNA copy number. CRISPR-Cas9 induced deletion of copies from tandem repeat regions, such as the 45s rDNA in plants, provides a new tool to understand the roles of such loc.

RESULTS

Cas9-induced double strand breaks at the 18S loci cause reduction of 45S rDNA Copy

Number

To determine the impact of reducing rDNA levels to their functional minimum in a model plant, we reduced the number of 45S rDNA copies in Arabidopsis thaliana using transgenerational Cas9 targeting of the 45S rDNA repeats (Figure 1). To achieve such copy number reductions, we designed a single guide RNA specific to the 18S locus within the 45S rDNA repeats (with no predicted off-target site) using the CRISPR-P online tool (http://cbi.hzau.edu.cn/crispr/). Using a previously described vector (Ryder et al., 2017;Wang et al., 2015), we developed a transgene cassette (pHEE-18S) containing the 18S guide RNA (gRNA). This transgene cassette allows expression of Cas9 exclusively in the egg cell (EC) of the haploid female gametophyte, where we hypothesised that Cas9 activity across the 45S rDNA repeats would generate either large deletions or insertions of the repeats, through the subsequent activity of the error-prone Non-Homologous End Joining (NHEJ) DNA repair pathway (Figure 1C) (Cubbon et al., 2018). Spatiotemporally localising Cas9 expression to the EC of the female gametophyte also allowed us to investigate the effects of copy number mutagenesis in the absence of Cas9 activity during other crucial stages of the life cycle such as meiosis, fertilisation and seed development. The T1 transformant seedlings were sown on hygromycin selective media and genotyped for 45S rDNA copy number by qPCR.

We recovered a population of T1 plants displaying large CNV in the 45S rDNA (Figure 2A), ranging from ~20% to 160% CN compared with wild type (WT). Whilst selection was initially performed to identify lines with CN loss (e.g. ~20% of WT copies, line #236 and #289, see Figure 2A) and CN gain, we determined that Cas9 activity predominantly causes transgenerational reduction of 45S CN. Hence a fixed increase in CN of 45S repeats could not be maintained over successive generations.

The Col-0 accession harbours four allelic variants of 45S rDNA which are associated with either NOR2 (VAR1 and 3) or NOR4 (VAR2, VAR3 and VAR4) (see Figure 2C and Pontvianne et al., 2010;Chandrasekhara et al., 2016). Investigation of genomic abundance of the 45S rDNA variants (Figure 2B) revealed that our mutagenesis approach caused a range of gene dosage variation of the 45S rDNA repeats across each independent line. Further, we investigated via RT-PCR whether the expression levels of the different 45S rDNA variants were altered and found qualitative changes in variant expression in the latest generation analysed (T7). For instance, we observed a strong expression signal of VAR4, the least abundant variant, in seedlings of line #236, while VAR1 appears more actively transcribed in rosettes of both LCN lines.

From the T1 generation, we selected two lines with particularly low copy number, lines #236 and #289 (Figure 2A, henceforth termed as LCN lines), and allowed these lines to self-fertilise for 6 generations, after which we recovered plants with CN variation ranging from 7-17% (line #289) and 11-21% (line #236) of WT (Figure 2B). In the LCN lines, (#236 and #289 T7 generations), effects on plant development were characterized from germination onwards (Supplemental Figure 1), where we found that the majority of seeds established a healthy seedling (~80%) and a similar percentage of seedlings successfully underwent their complete life cycle. In this study we have focussed on the WT-like plants to improve understanding of how a multicellular organism, with numerous cell-fate transitions throughout its lifecycle, can successfully complete its lifecycle with only ~10% of its 45S ribosomal RNA genes, under the growth conditions used.

To further confirm the reduction of 45S rDNA copies in line #236, we performed fluorescence in situ hybridisation (FISH) probing the 45S rDNA repeats in line #236 at the T7 generation. As previously shown, we found that both NOR2 and NOR4 are actively transcribed in young seedlings of A. thaliana Col-0, where one or two large signals localise at the nucleolus, indicating that both NORs are available for transcription (n=10 plants).

However, during the course of wild-type development, NOR2 becomes progressively silenced (Chandrasekhara et al., 2016) and moves further away from the nucleolus, eventually organising in form of chromocentres at the nuclear periphery. NOR4 remains in proximity to the nucleolus, typically localised at two large chromocentres on the sides of the nucleolus (n=12). (Pontvianne et al., 2013). While some chromatin condensation is still visible (i.e. the rounder signals within the nucleolus), this localisation of NOR2 within the nucleolus suggests that the 45S rRNA gene copies of NOR2 may remain available for transcription in line #236 (n=30) as a potential mechanism of gene dosage compensation.

While chromatin organization is strongly altered, rRNA homeostasis remains unchanged

We investigated whether transcription of rRNA or its accumulation is altered in the LCN plants, particularly during seedling development, when more copies of rRNA genes are actively transcribed. Hence, we quantified rRNA levels in line #236 by two methods. We performed an RNA gel blot using probes targeting either the functional rRNA (18S) or the precursor 45S (Figure 3B). Our t analysis shows that apart from minor differences in the 5' region of the transcript, there is no major alteration in rRNA level accumulation between the LCN lines and the WT control. We also performed an absolute quantification of rRNA molecules (using RNA spikes, see materials and methods, Figure 3C) on pooled seedlings of line #236 in a previous generation (T4), and found no difference in rRNA accumulation with WT. Cell size was also measured to exclude any underlying difference in global rRNA accumulation: leaf protoplasts were isolated and cell size measured. A 6% decrease in cell size was found in line #236 (Supplemental Figure 1) which is likely negligible when comparing global RNA accumulation per biomass.

Remarkably, this suggests that drastic reduction of 45S CN does not affect levels of rRNA accumulation. Hence, we hypothesised the presence of a gene dosage compensation mechanism affecting transcription which allows rRNA accumulation in the LCN lines to be similar to WT, despite a loss of >90% of 45S rRNA genes. We considered whether such a gene dosage compensation mechanism might be acting through higher rates of transcription at active repeats, and/or by chromatin reorganization allowing additional rDNA repeats to be activated.

To investigate the possibility of a gene dosage compensation mechanism involving increased transcription rates, we performed nuclear run-on transcription assays on line #236 (T7) and WT: we probed the 5' region of the transcript, which appears more abundant in #236 than in WT from our RNA gel blot analysis, but found no evidence of an increase in the global rate of transcription of rRNA in the CN depleted line (Figure 3D). Taken in combination with the steady state rRNA levels remaining the same, this indicates that rRNA stability remains similar in the CN depleted lines. However, it cannot be excluded that rate of transcription per activated 45S rRNA gene copy could potentially be increased, i.e. through increased Pol I loading or occupancy [START_REF] French | In Exponentially Growing Saccharomyces cerevisiae Cells, rRNA Synthesis Is Determined by the Summed RNA Polymerase I Loading Rate Rather than by the Number of Active Genes[END_REF]Lawrence et al., 2004).

We further investigated whether chromatin organization changes were occurring in this same line #236, which could allow for rRNA accumulation to be as high as WT. We employed chromatin immunoprecipitation (ChIP) against active (H3K9ac) and silencing (H3K9me2) marks which are known to be present on the rDNA repeats, as well as global H3 as a proxy for nucleosomal occupancy on the repeats (Benoit et al., 2019). qPCR was then used to estimate the enrichment of both marks across the 45S rRNA gene, relative to a euchromatic control (HXK1) and a heterochromatic control (Ta3 -LTR retrotransposon). We demonstrate that the enrichment of H3K9me2 was significantly reduced at the 45S rDNA loci in the LCN lines, while a trend of higher enrichment of H3K9Ac was also observed at the 3' end of the repeats (Figure 3E). Furthermore, the dramatic loss of H3 occupancy at the 5' region of the 45S rRNA gene suggests a more relaxed chromatin state of the regulatory regions which is consistent with active rDNA repeats devoid of nucleosomes. Taken together, our data indicate that gene dosage compensation of rRNA levels occurs despite rDNA copy number depletion. We suggest that this is most likely due to increased frequency of transcriptional activation of rDNA copies in LCN lines, occurring predominantly through a relaxation of chromatin organization and loss of silencing marks.

Long-read sequencing of two independent LCN lines indicates that genome architecture is largely preserved

To investigate the chromosomal integrity of the genome in our LCN lines, we performed long-read Nanopore sequencing (Table 1) to detect possible deletions, insertions or duplications in the copy number depleted lines #236 (T6) and #289 (T5) (Figure 4). Analysis of the coverage of the nanopore reads, of the CN depleted lines versus WT, identified five loci with a coverage fold change <0.5 or >1.5 (Figure 4A). The two loci with loss of coverage Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 in the CN depleted lines relative to WT mapped to the two 45S rDNA repeats loci annotated in TAIR10 (NOR4 is not annotated in TAIR10, with the 45S genes located in Chr2 and Chr3 in TAIR10 reference genome), confirming the loss of 45S rRNA genes in those two lines (~20% coverage in both lines versus WT, Figure 4D), consistent with our qPCR and FISH assays (Figure 2).

Of the three other loci displaying differences, two mapped to repetitive elements on Chromosomes 1 and 3 (92 bp AT repeat in Chromosome 1 starting at position 13,085,467 and a 70 bp CT repeat in Chromosome 3, starting at position 2,258,906). We consider that these likely represent background differences with the Col-0 WT line used as control for coverage analysis, which was not the original line used for pHEE401 transformation. A single, large duplication event (~220 kb, at positions 208,000 to 430,000 and containing 51 genes) was found in line #289 on chromosome 4 (Figure 4C), indicating that, at least for one of the two LCN lines, genome integrity has been affected. To ascertain that this duplication is correlated with loss of 45S rDNA, we designed primers spanning the duplication segment, which allowed us to confirm its appearance at the T5 generation (Supplemental Figure 2). This suggests that genome integrity has possibly been affected by loss of 45s rDNA, as described in another line which features loss of 45S rDNA (Picart-Picolo et al., 2020)

Genome-wide gene expression changes following 45S rDNA copy number reduction

While the impact of a major reduction of 45S copy number on genome architecture appears limited, we further analysed the transcriptome of 7-day-old seedlings (with three biological replicates) of the two LCN lines. This revealed the presence of transcriptome differences between WT seedlings and the two LCN lines (Figure 4A and Figure 4B). Indeed, the seedlings of the LCN lines #236 and #289 had 581 (233 down-and 348 up-regulated) and 565 (315 down-and 250 up-regulated) differentially expressed genes (DEGs, q-value < 0.05 and fold change > 1.5), respectively, when compared with WT. 50% of the DEGs were detected in both LCN lines (Figure 4B, Supplemental Data Set 1). Similarly, MapMan4 bin enrichment analysis identified some consistent responses between the two LCN lines (Usadel et al., 2009;Schwacke et al., 2019). We identified 12 bins enriched for DEGs in both LCN lines (Figure 4E; Supplemental Data Set 1), representing 54% of enriched bins in both lines #236 (22 bins) and #289 (22 bins), respectively. Among those enriched bins were functional processes and enzyme pathways such as cell wall organization, oxidoreductases (including the Cytochrome P450 family), biosynthesis Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 and degradation of glucosinolates and biosynthesis of pyruvate family amino acids. The enrichment for DEGs in these bins suggests a possible dysregulation of key pathways, such as glucosinolate metabolism or cell wall formation in the LCN lines. Genes involved in cell wall organization displayed opposite expression profiles between the two independent LCN lines, with an up-regulation and down-regulation of genes involved in this process in line #236 and #289, respectively. Finally, 35% of the genes located within the novel tandem duplication detected on Chromosome 4 of line #289 were found to be up-regulated (18 out of 51 genes, mean fold change of the 51 genes within locus = 1.5 Figure 4C) displaying a positive gene-dosage response. Conversely, only two genes in this locus were dysregulated in line #236 (2 downregulated genes out of 51). The mean fold change was 0.93 in line #236, which did not contain the tandem duplication (Figure 4C).

DISCUSSION

Severely depleted rDNA copy number is compatible with plant viability

We hypothesised that targeting the tandemly arrayed 45S rDNA repeats with Cas9 endonuclease in the female gametophyte would cause DSBs across the NORs, causing a transgenerational reduction of 45S rRNA genes. Supporting this we found a high degree of CNV in the transformant T1 population, indicating that Cas9 activity had caused both loss and gain of 45S rDNA copies as hypothesised. The Col-0 accession harbours ~375 copies per NOR. In our study, we aimed to reduce 45S CN by allowing Cas9 activity in the EC by maintaining these generations of these transformed lines until a reduction in copy number was no longer detected. This was achieved at the T7 generation in lines #236 and #289, where rDNA copy number plateaued at ~10%. By using transgenerational Cas9 activity, we have demonstrated that, by the 8 th generation, about 80% of plants in both #236 and #289 can fully undergo their life cycle despite having only ~10% of rDNA copies compared to WT (i.e. ~25/30 copies per NOR). This finding allows us to definitively conclude that as little as 50 rDNA copies (10 % of WT) per haploid genomes are sufficient for viability in A. thaliana accession Col-0, which was previously proposed to be 15-20% (Pontvianne et al., 2013) but impossible until now to causatively establish.

It remains to be determined whether the lowest functional CN is absolute or relative to the WT CN. Since A. thaliana accessions display large CNV, spanning from ~500 copies to Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 ~2500 (Rabanal et al., 2017), generating LCN lines in accessions with higher or lower 45S CN will help to address this question. In Arabidopsis thaliana, such a drastic loss of rDNA copies has only been observed in the fas1 fas2 mutant (Mozgova et al., 2010), which showed a 45S rDNA CN reduction of ~80% of WT. However the question of minimum 45S CN compatible with viability was confounded by the pleiotropic phenotypes exhibited by the fas1 fas2 mutant, such as loss of telomeres on all chromosomes (Mozgova et al., 2010). Our analysis of the 45S rDNA variants -which are associated with either NOR2 (VAR1 and VAR3) or NOR4 (VAR2, VAR3 and VAR4) (Chandrasekhara et al., 2016;Mohannath et al., 2016) suggests that loss of CN from either NOR occurred randomly in each independent line, as one could expect due to the sequence identity and repetitive nature of the 45S rDNA repeats.

While we demonstrate here that plants with as little as 50 rDNA copies per haploid genome are able to complete their lifecycle, the presence of ~20% of unviable seedlings within the progeny of the two independent low CN lines remains to be mechanistically characterized (Supplemental Figure 1). Further, we located the insertion site(s) of the transgenic Cas9 cassette in both lines using the nanopore data. We found homozygous insertions (2 insertions in #289, mapping both to chromosome 2, position 9.73 Mb and 8.58

Mb and 1 insertion in #236, mapping to chromosome 5, position 26.04 Mb (Figure 4A). In line #289, the first transgene insertion is located within the promoter of AT2G22860, and the second insertion in the 5' UTR of AT2G19880 but neither are associated with expression changes. In line #236, the insertion is located in an intergenic region between genes AT5G65165 and AT5G65170, none of which are differentially expressed. Hence, we consider that the insertion sites of the transgenic cassette are unlikely to affect the two low CN plant phenotypes and/or viability, and that the 20% of unviable seedlings found in both LCN lines is likely due to other mechanisms.

Altered chromatin organization of 45S rDNA ensures gene dosage compensation of rRNA transcription levels

One of the key questions we investigated was whether reduction in 45S rDNA copy number causes alteration to rRNA transcription levels, particularly during critical stages of development (e.g. seedling establishment). However, no difference in rRNA levels was found in either of the two LCN lines both for mature rRNA accumulation, as well as the precursor rRNA transcript. While major cell size differences could lead to erroneous absolute quantifications based on biomass, we only found a small (~6%) difference in cell size in line #236 vs WT (Supplemental Figure 1). This indicates that biomass-based quantification of rRNA transcript abundance is appropriate in this system, however this finding also presents interesting new avenues, to be explored in the future detailed phenotypic characterization of LCN lines.

Following the analysis of steady state levels, nuclear run-on assays were used to measure in vitro transcription rates of rRNA. This also indicated little difference in rDNA transcription rates between WT and LCN lines. To identify the underlying mechanism for gene dosage compensation of rRNA levels in the 45S rDNA copy number depleted lines we probed the epigenetic landscape of the rDNA locus and found statistically significant changes in the chromatin organization of the repeats. Our results suggest that the loss of enrichment in silencing marks such as H3K9me2 may be sufficient to allow for gene dosage compensation of rRNA transcription to occur in the LCN lines resulting in similar rRNA levels as in WT plants, despite the dramatic loss of 45S rRNA copy number. It is plausible that removal of silencing marks at rRNA allows for WT-like transcription levels of rRNA, without the apparent necessity of altering transcriptional rates. Despite the large reduction in rDNA copies, such increased levels of rRNA levels may arise through the action of enhanced histone demethylase activity at rDNA loci (Pontvianne et al., 2012). This hypothesis remains to be tested by investigating the levels of activity of different histone lysine methyltransferases and histone demethylases which could be involved, as well as probing global enrichments of heterochromatin in comparison to WT.

Taken together, our findings support a model for gene dosage compensation of rRNA transcript levels in A. thaliana through chromatin organization, where the transcriptional activity of Pol I itself is controlled largely through epigenetic modification, as previously suggested (Grummt and Pikaard, 2003). It cannot be excluded however that increased loading of Pol I on the rDNA promoter could also aid in maintenance of transcription, as reported in yeast [START_REF] French | In Exponentially Growing Saccharomyces cerevisiae Cells, rRNA Synthesis Is Determined by the Summed RNA Polymerase I Loading Rate Rather than by the Number of Active Genes[END_REF].

The relative proportion of active/inactive and silent rDNA genes is still unknown: the FISH images we collected suggest that a degree of chromatin condensation is still present (i.e. from the larger round signals inside the nucleolus). However chromocentres appear hollower and less condensed than in WT nuclei, similarly to what was observed in the atxr5 atxr6 mutants, which suffer a conspicuous loss of heterochromatin (Feng et al., 2017). We considered that Cas9 activity acting on the rDNA repeats in the female gametophyte, coupled with transgenerational loss of rDNA copies, could induce and aggravate genome instability through extensive DNA damage and loss of large amounts of genomic DNA. We used long-read Nanopore sequencing to determine any major chromosomal rearrangements, deletion or insertions. Other than a ~220 Kb segmental duplication on chromosome 4 of line #289, we found no major signs of genomic changes. However, one of the most up-regulated genes in both lines is AT3G61010 (log2 fold change > 7). This gene is part of a family of ferritin / ribonucleotide reductases genes that have been linked to genotoxic stress in Arabidopsis (Roa et al., 2009). The up-regulation of that gene could be a direct result of potential genotoxic stress induced by either the expression of the Cas9 transgene in the egg cell or by the loss of rDNA copies itself. Picart-Picolo et al, 2020 have reported that Line 6 from a fas1 fas2 population (Pavlištová et al., 2016), which bears 20% of rDNA genes, harbours a chromosome segment duplication observed at the heterochromatic knob region of chromosome 4 (Picart-Picolo et al., 2020b). Curiously, the duplication we retrieved in line #289 appears in a neighbouring region of the heterochromatic knob on the same chromosome. On the other hand, the duplication in Line 6 is accompanied by a large number of smaller tandem duplications across the genome, which were not found in our LCN lines.

This leaves open the question of further determination of the level of genomic instability in lines #236 and #289. However, these segmental duplications occurred in neighbouring regions of two independent lines (#289 and Line 6), obtained in two completely different approaches. This suggests the possible presence of a fragile region (e.g. with rDNA copy number dependent fragile sites) which might be crucial for genomic structural integrity and which becomes duplicated following large loss of repetitive DNA (or rDNA) content.

The heterochromatic knob on chromosome 4 hk4s is the only region of the A. thaliana genome which features a Topologically Associated Domain (TAD)-like structure (Grob et al., 2013;Pontvianne and Grob, 2020), suggesting that this region could play an important role in the structural dynamics of the genome of A. thaliana. Furthermore, our results suggest that the genome appears well-buffered against loss of rDNA copies, however DNA Damage & Repair (DDR) pathways in these lines remain to be investigated, together with global Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 heterochromatin enrichment, to provide a full view of the landscape of genome stability in these lines.

We consider that these novel rDNA depleted lines provide exciting new tools to i) investigate the extent of genomic instability singularly caused by loss of 45S rDNA copy numbers, ii) to characterize the effects of extensive DNA damage on rDNA loci in the female gametophyte, iii) to elucidate the role of 45S rDNA in plant cellular senescence, and iv) to investigate the dynamics of rDNA CNV when backcrossing to WT and/or with other accessions with reduced 45S CN.

The transcriptome of LCN lines shows significant effects of 45S rDNA copy number depletion

The transcriptome analysis of both independent LCN lines revealed ~570 dysregulated protein-coding genes compared with WT with a fold change > 1.5 (Figure 4D). Such transcriptome difference could be caused by the loss of rDNA copies itself or by chromatin remodelling due to nucleolus reorganization, possibly to facilitate the activation of all remaining 45S rDNA copies. In this context, we hypothesised that the epigenetic machinery may be responsible for maintaining ribosomal RNA homeostasis (Figure 3), through aberrant deposition/maintenance of epigenetic marks throughout the rest of the genome. Genome-wide chromatin immunoprecipitation on histones and histone marks can be used to better characterize the transcriptome effects observed due to rDNA depletion. Genome-wide aberrant epigenetic marks could lead to changes in gene expression between independent lines. Furthermore, as shown in Figure 4B and4E, 311 of the differentially expressed genes are shared between the two independent LCN lines, leading to twelve enriched common pathways (Figure 4E). Those results indicate that 45S rDNA contraction likely leads to transcriptome changes, with a common set of rDNA CN-sensitive genes and pathways. Of note also is the dysregulation of glucosinolate metabolism in both lines, suggesting elevated cellular stress which might be connected to pathogen response pathways, as observed in Line 6. However, the precise cause and effect of the dysregulation of glucosinolate biosynthesis and degradation genes remains to be characterized. Recently, glucosinolate re-mobilisation following Carbon starvation was observed in Arabidopsis (Brandt et al., 2018), indicating that the glucosinolate metabolism is likely to play a role in adaptation to abiotic stress, in addition to its well-known role in defense. , 2020). The DEGs appear to be correlated with the occurrence of duplication events. Our study supports this hypothesis, since the single duplication event found in line #289 indeed correlates with higher expression levels of 35% of the genes within the duplicated locus (18/51 genes up-regulated, Figure 4C, mean fold change of genes within the region of 1.5 versus WT). However, this cannot explain the deregulation of the remaining ~550 deregulated genes, or indeed of the 581 DEGs in line #236 which does not contain duplications despite a similarly low 45S copy number. Conversely, only 16 DEGs are common between line 6 and our low CN lines #236 and #289. Hence, our analysis reveals an effect on the transcriptome arising simply from rDNA copy number depletion, similar to effects observed in Drosophila (Paredes et al., 2011). The mechanistic basis for such transcriptome deregulation remains to be further investigated in future studies.

In this study we present an approach and tools to elucidate the precise role of 45S rDNA copy number in plants. Indeed, our approach represents a simple, "clean" method to modify rDNA copy number without off-target modifications. A single transformation event allowed us to generate a population of plants with rDNA copy number ranging from ~20% to 160% of those of WT. Intriguingly, we present evidence that gene dosage compensation of rRNA levels is tightly regulated, most likely by chromatin remodelling, with a similar rRNA accumulation in LCN lines despite the loss of hundreds of 45S copies. The possible impacts of 45S copy number depletion on protein dynamics, genome integrity, plant reproduction, development and fitness are fertile avenues for further investigation using the novel 45S LCN lines (Figure 5).

Gene dosage compensation systems have arisen in different evolutionary contexts to adjust transcript levels in response to changes in gene dosage (e.g. as arises on sex chromosomes) (Veitia et al., 2013). Our results reveal a gene dosage compensation mechanism in Arabidopsis for rRNA transcript levels when rDNA copy number is altered, which differs from classical genetic studies in maize where changing the dosage of NORs was found to have a dosage effect on NOR RNA levels (Lin, 1955). Conversely, Buescher et al., 1984 found that a maize line with two NORs did not affect the levels of rRNA, which is consistent with our findings of gene dosage compensation effects on rRNA levels in Arabidopsis thaliana.

Finally, our study adds to the growing body of knowledge on use of CRISPR-Cas for study of plant functional genomics (Lynagh et al., 2018;Beying et al., 2020). Our novel approach paves the way for research on understanding roles of crop rDNA repeats, including in species where genomic architecture is more complex than Arabidopsis (e.g. with a higher abundance of repetitive elements, allopolyploidy, high functional redundancy).

METHODS

Plant materials and growth conditions

Arabidopsis thaliana accession Col-0 was used for all experiments. Unless otherwise specified, plants were grown on either 1% w/v MS medium (seedling samples) or sown onto a mixture of soil, perlite, vermiculite (5:1:1) (leaves samples). Long day conditions (16 hours light at 21 °C, 8 hours dark at 18 °C) were used with white LED lights with Photosynthetically Active Radiation (PAR) of ~100 µmol.m -2 .s -1 photons.

Generation of T1 Cas9 transgenic plants

The vector used in this study (pHEE401) has been previously described in (Wang et al., 2015). 18S-specific single guide RNA was designed using CRISPR-P (http://cbi.hzau.edu.cn/crispr/) and cloned using the following oligonucleotides: 18S Fwd ATTGATACGCTCCTGGTCTTAAT and 18S Rev AAACATTAAGACCAGGAGCGTAT.

Annealed oligos were directly inserted into pHEE401 using BsaI cut-ligation. Recombinant plasmids were electroporated into Agrobacterium tumefasciens, and Col-0 plants transformed using the floral dip method (Clough and Bent 1998). Transformed seeds were surfacesterilised using chlorine gas and sown onto 0.5% MS medium with 20 µg/mL of Hygromycin B. Transformant seedlings were selected as described in (Harrison et al., 2006).

DNA isolation and Copy Number quantification

For genomic DNA extractions, single seedlings or 10 mm 2 piece of rosette leaves were flash frozen in liquid nitrogen and subsequently ground into a fine powder using a QIAGEN TissueLyser. Pulverized tissue was incubated at 65 °C for 1 hour in 2X CTAB extraction buffer (2% Cetyl-trimethyl-ammonium bromide, 1.4 M NaCl, 100 mM Tris-HCl pH 8.0, 20 mM EDTA, 5 mM β-mercaptoethanol). Samples were incubated on ice for 2 minutes followed by addition of 1 volume of chloroform. Samples were vortexed thoroughly and centrifuged, supernatant was collected, and DNA was precipitated with ice-cold isopropanol. Samples were centrifuged and pellet washed with 70% Ethanol and air dried. DNA was suspended in 50 µL of sterile ddH 2 O. To quantify the number of 45S copies in the genome of WT and LCN lines, we followed the method described in (Rabanal et al., 2017), using quantitative real time PCR (Bio-Rad CFX™ thermocycler) in 5 µL reaction volume with SensiFast Sybr No-ROX kit (Bioline). Copy number per haploid genome was calculated by comparing the threshold cycle from the 18S subunit of the 45S repeats with the threshold cycle from the single copy genes TZF1 and HXK1 using the following formula 𝐶𝑜𝑝𝑦𝑁𝑢𝑚𝑏𝑒𝑟𝑝𝑒𝑟ℎ𝑎𝑝𝑙𝑜𝑖𝑑𝑔𝑒𝑛𝑜𝑚𝑒 = 2 2

Copy Number per haploid genome was normalized on WT values to express the copy number in our various lines as a percentage of WT copy number.

Absolute quantification of ribosomal RNA molecules

For the quantification of rRNA molecules, we used the protocol described in (Piques et al., 2009;Pal et al., 2013;Vaid et al., 2020;Ishihara et al., 2017) with modifications. 7 days old seedlings of WT and #236 lines were harvested and immediately flash-frozen in liquid nitrogen, with three biological replicates. The frozen biomass was ground into a fine powder using a ball mill (QIAGEN TissueLyser II). Aliquots of ~20 mg of ground frozen tissue were resuspended in 1 mL of Bioline ISOLATE II lysis buffer, and the suspension diluted 10 fold in lysis buffer. 10 µL of that final dilution was used for RNA extraction. A RNA spike mix (Ambion ArrayControl spikes 2,5,6,7 and 8) containing a known number of RNA molecules in the dynamic range of 3.9 x 10 7 to 7. amount of biomass in the aliquots used for RNA extraction (Ishihara et al., 2017) . Since rRNA quantity could be affected by a reduction in rDNA copy number. The spikes were added prior to RNA extraction and the results normalized on a biomass basis. Adding the spikes to a fixed amount of RNA rather than biomass could give erroneous results in the case of a differential expression of ribosomal RNAs. Indeed, since rRNAs represent the majority of total RNAs within a cell, normalizing by the amount of total RNA (as is done in RNA gel blots) could potentially erase the possible intrinsic difference in rRNA levels between samples, leading to the risk of false negative results.

3' EST variant assays

gDNA from individual rosette leaves of T4 and T7 generations was extracted by CTAB as previously described, from lines #236 and #289. PCR was carried out as follows using primers P3-P4 (see Primers Table ) and MyTaq Red Mix by Bioline®: 94° for 5', 25 cycles of 94° for 30'', 58° for 30'' and 72° for 60''.

For cDNA analysis of variants, RNA extracted as previously described, and treated with DNAse I by Invitrogen® as per manufacturer's instruction. 1.5 µg RNA were retrotranscribed with SuperScript™ III First-Strand Synthesis System by Invitrogen®, following manufacturers instruction, PCR for variant analysis was carried out as follows using P3-P4 (see Primers Table ) and MyTaq Red Mix by Bioline®: 94° for 5', 29 cycles of 94° for 30'', 58° for 30'' and 72° for 60''.

ACT2 was amplified from cDNA using primers listed in primers table and PCR was carried out as follows: 94° for 5', 26 cycles of 94° for 30'', 58° for 30'', 72°.

RNA Gel Blot and Nuclear Run-on Assay

For RNA gel blot studies, total RNA was isolated from 7-day-old (after germination) pooled

Col-0 and LCN seedlings. RNA gel blot analyses were performed using 4 μg of total RNA for each sample. 32 P-labelled DNA probes were generated using primers listed in Table 1. For Run-on transcription assays, nuclei were extracted from 1.2 g of 7-day-old pooled seedlings and isolated according to (Folta and Kaufman, 2007) 1 µg of DNA probes (A', A and B, Figure 2A) were blotted onto a Hybond-N + membrane (Amersham, Little Chalfont, UK) and hybridized with 32 P -labelled RNA. DNA probes were generated using primers as listed in Supplemental Table 1.

Nanopore sequencing and data analysis

Genomic DNA preparation was performed as previously described. DNA was further purified using Genomic DNA Clean & Concentrator kit (Zymo Research, USA). Qubit (dsDNA High Sensitivity (Thermo Fisher Scientific, USA) quantification was performed before library preparation using the 1D Genomic DNA by ligation kit SQK-LSK109 (Oxford Nanopore Technologies, UK), following manufacturer's instructions. The R9.5 ONT flow-cell FLO-MIN106D (Oxford Nanopore Technologies, UK) was used with MinKNOW version 3.6.5, Guppy 3.2.10 through the software MinION release 19.12.5. Reads were aligned onto the Arabidopsis genome using minimap2 (Li, 2018). Next, the Arabidopsis TAIR 10 genome was split into 100 kb windows using bedtools makewindows and the coverage counts of WT and lines #236 and #289 was calculated against the 100 kb windows using bedtools coverage (Quinlan and Hall, 2010) and normalised by the total number of reads. Normalised counts were analysed using R, and the fold change per window versus WT calculated by dividing normalised counts of the LCN lines against WT. The final normalised Fold Change was calculated by dividing the fold change per window by the geometric mean of the fold change of all 100 kb windows for lines #236 and #289. The normalised fold change per window was visualised using Circos (Krzywinski et al., 2009) and 100 kb windows with a fold change <0.5 or >1.5 highlighted as having differential coverage from TAIR10 genome.

RNA extraction and transcriptome analysis

RNA extraction was performed on pooled seedlings (collected in triplicates at 7 days after germination) of Col-0 WT, #236 and #289 with WT phenotype grown on three Murashige and Skoog plates using Aurum™ Total RNA Mini Kit (Bio-Rad). RNA Integrity and purity were verified by gel electrophoresis and Nanodrop quantification, 10 µg of total RNA was concentrated and purified using RNA Clean & Concentrator-5™ (Zymo Research) and sent for sequencing to Novogene UK. Transcriptome sequencing was performed on mRNAenriched RNA libraries using Illumina technology, 150 bp paired end reads, generating >21 million reads for each library. The reads were trimmed of residual adaptor sequences using Trimmomatic (Bolger et al., 2014) and transcript abundance estimated using Kallisto (Bray et al., 2016), using the latest reannotation of Arabidopsis thaliana reference transcriptome (Araport 11). Differential expression was assessed by Likelihood Ratio Testing with Sleuth in R (Pimentel et al., 2017) using "genotype" (i.e #236, #289 and WT) as factor in the full model, against a reduced model without genotype information. The minimum detection frequency filter was set to > 0.3 to allow for detection of transcripts detected in at least one of the three genotypes. Transcripts were aggregated into genes during the Sleuth analysis. The comparison of the fits between full and reduced models for the abundance of each gene highlights those whose expression is more likely determined by the genotype than by the null hypothesis. Finally, log2 fold change approximations were extracted relative to WT using a Wald test, also in Sleuth. Genes with an absolute log2 fold change (b-value) > 0.58, representing a fold change of > 1.5 were selected and p-value correction performed using Benjamini-Hochberg on the p-values of the genes passing the fold change threshold. The fold change cut-off was selected since individual genotypes were grown on separate plates, potentially creating subtle batch effects. Hence, we only considered targets with a fold change greater than 1.5 as differentially expressed. In summary, genes with i) a likelihood ratio q-value < 0.05, ii) a Wald test q-value < 0.05 and iii) an absolute log2 fold change (bvalue) > 0.58 were deemed significantly dysregulated. For MapMan analysis (Usadel et al., 2009), the lists of up and down-regulated genes, as well as the entire list of expressed genes (background) was used for bin enrichment. Genes were assigned into a MapMan bin structure using Mercator4 (Schwacke et al., 2019). A MapMan bin was classified as enriched if the number of genes belonging to that bin in the up/down lists was statistically higher than expected from the background, using a two-tailed Fisher's exact test and adjusted for multiple testing using Benjamini-Hochberg correction. Only bins enriched in both lines were retained for the analysis, to highlight possible similarities between the two independent LCN lines.

The Log2 Fold Enrichment was calculated by dividing the number of observed versus expected genes in each bin.

Fluorescence in situ hybridization

Nuclei were extracted from either pooled seedlings (5 days after sowing) or rosette leaves (20 days after germination) of WT and line #236. Tissues were fixed in cold 4% paraformaldehyde for 15 minutes at 4 °C and washed twice in Tris Buffer (10 mM Tris-HCl pH 7.5, 10 mM EDTA, 100 mM NaCl). Nuclei were extracted by chopping tissues with a razor maintaining the same exposure time and conditions.

Chromatin immunoprecipitation

Seedlings of WT and line #236 were grown for 15 days after sowing and the areal parts including the meristem were pooled and collected in triplicates in ice-cold ddH 2 O. Chromatin crosslinking was performed in 1% formaldehyde under vacuum for 10'. Quenching was carried out with 0.1 M glycine and samples were placed under vacuum for further 5' and washed twice with ice-cold ddH 2 O. Tissues were flash-frozen and ground in liquid N 2 and chromatin immunoprecipitation was carried out as described in (Benoit et al., 2019).

Chromatin was precipitated using antibodies for H3 (Abcam, ab1791), H3K9me2 (Abcam, ab1220) and H3K9Ac (Actif Motif AB_2561017). qPCR analysis was performed with a Roche Lightcycler 480 in 10 µL reactions of Roche Lightcycler 480 SYBR Green I Master using either primers for controls (Ta3 and HXK1) or for different regions of the 45S locus.

Primers are detailed in Primers Table . To adjust for different IP efficiencies between the independent ChIP experiments, IP levels were normalised to input for H3 and normalised to Ta3 for H3K9me2 and to HXK1 for H3K9Ac.

Cell size measurement

Leaf protoplasts were extracted from 10 leaves of 4 week-old rosette of A. thaliana plants in WT and #236 background (n=5 from each genotype) using the "Tape-Arabidopsis Sandwich method" described in Wu et al., 2009. For cell wall digestion, 1.25% of Cellulase Onozuka R-10 (Serva) and 0.3% of Macerozyme R-10 (Serva) were used. Protoplasts were washed with W5 solution (5 mM Glucose, 154 mM NaCl, 125 mM CaCl 2 , 5mM KCl, 5 mM MES (Sigma-Aldrich) and imaged using an Olympus BX61. 37 and 45 pictures containing ~10,000 protoplasts were analysed using ImageJ (https://imagej.nih.gov/ij/) to calculate size. The images were segmented using the "find edges" process and a colour threshold. Touching cells were separated from the segmented binary images using the watershed process. Cell size was then measured using the analyse particles option.

  Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 We observed a strong reduction in signal intensity (Figure 2D, n=50), confirming the drastic loss of 45S rDNA copies in line #236. NOR signals are localised at the nucleolus, although the fluorescence appears more diffuse within the nucleolus itself, indicating loss of the highly condensed chromatin organization of the NORs. During vegetative development of line #236, the 45S signals still localise inside the nucleolus, indicating that copies from both NORs appear to remain available for transcription, similar to what is observed in the fas1 fas2 mutant

  Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 Chromosome-segment duplication: a symptom of genome instability or a mechanism to preserve genome integrity?

  Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 Picart-Picolo et al, 2020 highlighted that Line 6 displays ~ 350 DEGs (Picart-Picolo et al.

  Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021

Figure 1 :

 1 Figure 1: Overview of 45S rDNA loci, dynamics of transcription and ribosome assembly, and mutagenesis concept. (A) The 45S rRNA gene comprises two External Transcribed Spacers (ETS), two Internal Transcribed Spacers (ITS) and the 18S, 5.8S and 25S (28S in fungi and animals) ribosomal components. A single guide RNA (sgRNA) was designed specifically to guide Cas9 to the 18S locus and generate Double Strand Breaks (DSB) across the 45S repeats. (B) Transcription of 45S rRNA occurs in the nucleolus, the largest structure in the nucleus. Assembly of ribosomes requires import of the 5S transcript from the nucleoplasm, as well as import of ribosomal proteins into the nucleolus. LSU=Large Subunit; SSU= Small Subunit. (C) CRISPR-Cas9 mutagenesis: Random DSB (grey) were generated along the 45S repeats targeting the 18S locus within the 45S rRNA gene. Two outcomes were expected: deletions of large numbers of repeats or insertion of supernumerary copies (blue) by the action of Non-Homologous End Joining (NHEJ) DNA repair machinery. AC = Antipodal cells, V = Vacuole, CCN = Central cell nucleus, SC = Synergid cells.

Figure 2 :

 2 Figure 2: Selection of Low 45S Copy Number lines, 3' ETS variant analysis and fluorescence in situ hybridization at T7 generation. (A) 45S rDNA Relative Copy Number (CN) of T1 transformant population quantified by qPCR against the single copy genes HXK1 and TZF1. Blue bar and green bar indicate lines #289 and #236 respectively. Highlighted are LCN lines on which all experiments were performed. (B) 45S rDNA Relative Copy Number by quantification by qPCR of T4 to T7 generations. 11 individual plants of lines #236 and #289 were selected per generation and 1 WT Col-0 control. Bars represent standard error of technical qPCR replicates. (C) Schematic representation of 3' allelic variants, adapted from Pontvianne, 2010. Analysis of genomic and gene dosage of 45S rDNA variants in WT and

Figure 3 :

 3 Figure 3: Synthesis of 45S rRNA appears largely regulated by chromatin organization. (A) 45S rDNA locus. Letters indicate the probes used for RNA gel blots (Panel B) and transcription run-on (Panel C). Stars represent regions amplified by ChIP-qPCR. (B) RNA gel blot analysis shows accumulation of 45S and other ribosomal RNAs in WT Col-0 and the LCN lines #236 and #289, letters indicate probes used as shown in panel A, Methylene Blue is shown as loading control. Worth noticing that the plastid 16S and 23S rRNAs result equally represented in the LCN mutants and show a WT-like stoichiometric ratio with the 45S-derived rRNAs. (C) Absolute quantification of 18S and 25S rRNA molecules in WT and #236 (T4 generation). No difference in the accumulation of either was detected (Student t test, bars represent standard error among biological replicates, n=3 biological replicates) (D) Nuclear transcription run-on assay shows transcription rate for 45S rRNA in WT Col-0 and the LCN lines #236 and #289. ACT2 transcript was used as control. (E) ChIP-qPCR shows differential enrichment of global H3, H3K9me2 (silencing mark) and H3K9Ac (active mark) across the 45S loci. Ta3 and HXK1 were used as controls for a silent retrotransposon and a transcriptionally active gene respectively. H3 occupancy (left) was determined relative to input. Fold enrichments for H3K9me2 (middle) were normalised against heterochromatic control Ta3; fold enrichments for H3K9Ac (right) were normalised against euchromatic control HXK1. (Student t test, bars indicate standard error among biological replicates, **=p<0.01, *=p<0.05, n=3 biological replicates, no antibody control = average of 45S no antibody amplicons).

Figure 4 :

 4 Figure 4: Genome integrity and transcriptome effects of 45S contraction. (A) DNA coverage and gene expression analysis of two independent LCN lines compared with WT. For both LCN lines, outer layer represents Fold change of the DNA coverage compared with WT, inner layer the identified differentially expressed genes (q-value < 0.05, Fold Change > 1.5). Coverage > 1.5 or < 0.5 versus WT, representing duplications and deletions, respectively, are highlighted in red. (B) Euler diagram representing the number of Differentially Expressed Genes (DEGs) between the two LCN lines and WT. (C) Close-up of the duplicated region identified in line 289. Y axis represent the fold coverage versus WT, X axis represents the position on Chromosome 4. Bars indicate genes, colour coded based on their expression level vs WT. (D) Coverage of 45S rDNA gene loci in LCN lines relative toWT assessed by Nanopore sequencing. E) MapMan enrichment analysis of bins significantly enriched for up or down-regulated genes in both independent LCN lines. X axis represent the fold enrichment of each significant bin, left: down-regulated genes, right: up-regulated genes, the size of the circles corresponds to the -log 10 adjusted p-value.

Figure 5 :

 5 Figure 5: Effects of rDNA Copy Number reduction on rRNA transcription. Analysis of rRNA transcription in our lines revealed that gene dosage may be maintained by rendering both NORs competent for transcription potentially through the removal of silencing histone marks. The rDNA LCN lines generated can be further used to advance understanding of the functional role of rRNA genes in cellular and developmental processes.

Figure 1 :Figure 2 :Figure 3 :C

 123 Figure 1: Overview of 45S rDNA loci, dynamics of transcription and ribosome assembly, and mutagenesis concept. (A) The 45S rRNA gene comprises two External Transcribed Spacers (ETS), two Internal Transcribed Spacers (ITS) and the 18S, 5.8S and 25S (28S in fungi and animals) ribosomal components. A single guide RNA (sgRNA) was designed specifically to guide Cas9 to the 18S locus and generate Double Strand Breaks (DSB) across the 45S repeats. (B) Transcription of 45S rRNA occurs in the nucleolus, the largest structure in the nucleus. Assembly of ribosomes requires import of the 5S transcript from the nucleoplasm, as well as import of ribosomal proteins into the nucleolus. LSU=Large Subunit; SSU= Small Subunit. (C) CRISPR-Cas9 mutagenesis: Random DSB (grey) were generated along the 45S repeats targeting the 18S locus within the 45S rRNA gene. Two outcomes were expected: deletions of large numbers of repeats or insertion of supernumerary copies (blue) by the action of Non-Homologous End Joining (NHEJ) DNA repair machinery. AC = Antipodal cells, V = Vacuole, CCN = Central cell nucleus, SC = Synergid cells. Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021

Figure 4 :Figure 5 :

 45 Figure4: Genome integrity and transcriptome effects of 45S contraction.A) DNA coverage and gene expression analysis of two independent LCN lines compared with WT. For both LCN lines, outer layer represents Fold change of the DNA coverage compared with WT, inner layer the identified differentially expressed genes (q-value < 0.05, Fold Change > 1.5). Coverage > 1.5 or < 0.5 versus WT, representing duplications and deletions, respectively, are highlighted in red. B) Venn diagram representing the number of Differentially Expressed Genes (DEGs) between the two LCN lines and WT. C) Close-up of the duplicated region identified in line 289. Y axis represent the fold coverage versus WT, X axis represents the position on Chromosome 4. Bars indicate genes, colour coded based on their expression level vs WT. D) Coverage of 45S rDNA gene loci annotated in TAIR10 in LCN lines relative to WT assessed by Nanopore sequencing. E) MapMan enrichment analysis of bins significantly enriched for up or down-regulated genes in both independent LCN lines. X axis represent the log 2 fold enrichment of each significant bin, left: down-regulated genes, right: up-regulated genes, the size of the circles corresponds to the -log 10 adjusted p-value.

  Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021

	13 x 10 8 per µL was added to the diluted tissue prior to
	RNA extraction. RNA extraction was performed using Bioline ISOLATE II RNA kit
	following manufacturer's instruction. RNAs were subsequently treated with DNAse1 (Sigma-
	Aldrich), and cDNA synthesis performed on 15 µL RNAs using Bioline SensiFast cDNA
	synthesis kit with random hexamers. The cDNAs were diluted 50-fold with sterile miliQ
	water before RT-qPCRs. qPCRs were performed in 5 µL reactions containing 1 µL of cDNA
	and 4 µL of SensiFast Sybr No-ROX kit (Bioline). Threshold (Ct) values obtained from the
	qPCR reactions against each RNA spike were used to generate a standard curve for each
	sample (R

2 > 0.96 for all reactions). Ct values from the 18S and 25S qPCR reactions were compared with the slope and intercept to obtain the absolute amount of RNA molecules for both rRNA transcripts. Finally, the number of rRNA molecules was normalised by the

  Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 blade in Galbraith Buffer(Galbraith et al., 1983) and filtered through a 30 µM disposable filter (CellTrics, Sysmex, Germany). Filtrate was dripped onto a microscope slide and, after air-drying, slides were post-fixed in 2% formaldehyde in PBS. Slides were then washed in 2x SSC and treated with RNAse (1 hour at 37 °C) after which they were washed again in 2x

	SSC twice and further fixed in 2% formaldehyde. Following fixation, slides were washed
	with distilled H 2 O and dehydrated in ethanol series (70%, 90% and 100%), air dried and
	hybridised with the appropriate probes. 45S rDNA probes were generated by amplification of
	18S and 25S from Col-0 genomic DNA and cloned into pGEM-T Easy (Promega). 1 µg of
	each probe was labelled with Cy3-conjugated dUTPs (GE Healthcare) by Nick Translation
	(Roche) according to manufacturer instructions. Probes were resuspended at a concentration
	of 10 ng/µL in Hybridisation Buffer (50% deionised formamide, 2X SSC, 10 mM Sodium
	phosphate pH 7, 20% dextran sulphate), and further applied onto the slides which were
	denatured at 80 °C for 2 minutes and hybridised overnight at 37 °C. Post-hybridisation
	washes were carried out at 60 °C twice in 2x SSC, slides were mounted in Vectashield

containing DAPI (Vector laboratories). Imaging was performed with a Zeiss Axioimager and captured with ORCA-Flash4.0 V2 Digital CMOS camera C11440 (HAMA-MATSU)

Table 1 : Summary statistics of the nanopore sequencing runs Sample Number of reads

 1

			Mean read		Mean	Mean mapping
			length (bp)	Mapped %	coverage	quality
	Line 236	1,043,643.00	3,868	98.01	30.17	55.95
	Line 289	1,237,074.00	4,097	96.74	37.73	56.19
	WT	569,275.00	6,570	92.21	28.97	54.52

  Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021

Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021

Google Scholar: Author Only Title Only Author and Title Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021

ACKNOWLEDGEMENTS

Downloaded from https://academic.oup.com/plcell/advance-article/doi/10.1093/plcell/koab020/6126391 by guest on 05 February 2021 This work was supported by a Science Foundation Ireland Principal Investigator Grant (13/IA/1820) to CS. FL is supported by the Irish Research Council Postgraduate Fellowship GOIPG/2018/1960. FL and this work were partially funded by Cost Action CA16212 "Impacts of Chromatin Domains on Plant Phenotypes (INDEPTH)" (FL, FP and AVP are members of INDEPTH). FL was partly supported by EMBO Short Term Fellowship 8312. Special thanks to Sara Farrona for insightful discussions, Sophie Desset for advice with FISH and Christel Llauro for Nanopore library preparation and sequencing and Ariadna Picart-Picolo and Marie Mirouze for help with data analysis. Parts of Figures 1 and 5 were created with BioRender.

Data availability

Sequencing files for the nanopore and RNA-sequencing reads are available in Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) under BioProject ID PRJNA640267

Accession numbers

Gene models used in this article can be found in the Arabidopsis Genome Initiative database under the following accession numbers: HXK1 (AT4G29130), 18S (AT2G010101), TZF1 (AT2G25900).

Disclosures

The authors declare no conflict of interest.

SUPPLEMENTAL DATA

Supplemental Figure 1 -Plant viability in LCN lines.