Yaqun Liu

Xia Wu

Yves D'aubenton-Carafa

Claude Thermes

Chun-Long Chen
email: chunlong.chen@curie.fr

OKseqHMM: a genome-wide replication fork directionality analysis toolkit

Keywords: DNA replication, Replication fork direction, Replication initiation, Transcription-Replication conflicts, Hidden Markov Model Genome analysis

Motivation: During each cell division, tens of thousands of DNA replication origins are coordinately activated to ensure the complete duplication of the entire human genome.

However, the progression of replication forks can be challenged by numerous factors. One such factor is transcription-replication conflicts (TRC), which can either be co-directional or headon with the latter being revealed as more dangerous for genome integrity.

Results:

In order to study the direction of replication fork movement and TRC, we developed a bioinformatics tool, called OKseqHMM, to directly measure the genome-wide replication fork directionality (RFD) as well as replication initiation and termination from data obtained by Okazaki fragment sequencing (OK-Seq) and related techniques.

Availability and Implementation:

We have gathered and analyzed OK-seq data from a large number of organisms including yeast, mouse and human, to generate high-quality RFD profiles and determine initiation zones and termination zones by using Hidden Markov Model (HMM) algorithm (all tools and data are available at https://github.com/CL-CHEN-Lab/OK-Seq). In addition, we have extended our analysis to data obtained by related techniques, such as eSPAN and TrAEL-seq, which also contain RFD information. Our works, therefore, provide an important tool and resource for the community to further study TRC and genome instability, in a wide range of cell line models and growth conditions, which is of prime importance for human health.

Introduction

The faithful transmission of genetic information to daughter cells is crucial for maintaining genome stability. In humans, at each cell division, tens of thousands of replication origins need to be coordinately activated to ensure the complete duplication of >6 billion base pairs (bp) of the human genome. However, the DNA replication program is routinely exposed to endogenous and exogenous stresses, which play an important role in many human diseases. In . CC-BY-NC-ND 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

The copyright holder for this this version posted January 13, 2022. ; https://doi.org/10.1101/2022.01.12.476022 doi: bioRxiv preprint 3 particular, replication stress-induced genome alterations can represent an important early cause of cancer [START_REF] Gnan | The impact of transcription-mediated replication stress on genome instability and human disease[END_REF].

The progression of replication forks can be challenged by numerous factors. One such factor is transcription-replication conflicts (TRC) since the replication and transcription machineries share the same DNA template. TRC can either be co-directional (CD) or head-on (HO), and the latter has been revealed as more detrimental for genome integrity [START_REF] Hamperl | Transcription-Replication Conflict Orientation Modulates R-Loop Levels and Activates Distinct DNA Damage Responses[END_REF].

Previous bioinformatics analyses have revealed that, in large numbers of genomes from bacteria [START_REF] Merrikh | Spatial and Temporal Control of Evolution through Replication-Transcription Conflicts[END_REF] to human [START_REF] Chen | Replication-associated mutational asymmetry in the human genome[END_REF][START_REF] Huvet | Human gene organization driven by the coordination of replication and transcription[END_REF], most genes are codirectionally oriented with replication forks to avoid the more deleterious HO TRC. Recently, a new method to directly measure the genome-wide replication fork directionality (RFD) along the human genome by sequencing of Okazaki fragments (OK-Seq) [START_REF] Petryk | Replication landscape of the human genome[END_REF], which are present only on the lagging replicating strand, allows quantitatively analyzing and accurately detecting replication initiation and termination. The analysis of OK-seq data of human cells has also demonstrated a significant co-direction of replication fork progression with gene transcription within active genes [START_REF] Petryk | Replication landscape of the human genome[END_REF].

More and more techniques are now being developed, for instance, Pu-seq [START_REF] Daigaku | A global profile of replicative polymerase usage[END_REF], eSPAN [START_REF] Li | DNA polymerase α interacts with H3-H4 and facilitates the transfer of parental histones to lagging strands[END_REF], SCAR-seq [START_REF] Petryk | MCM2 promotes symmetric inheritance of modified histones during DNA replication[END_REF], GLOE-seq [START_REF] Sriramachandran | Genome-wide Nucleotide-Resolution Mapping of DNA Replication Patterns, Single-Strand Breaks, and Lesions by GLOE-Seq[END_REF] and TrAEL-seq [START_REF] Kara | Genome-wide analysis of DNA replication and DNA double-strand breaks using TrAEL-seq[END_REF], which also provide genome-wide RFD information.

Moreover, in recent years, strong evidence shows that replication-and transcription-related mutational asymmetries are widespread across cancer development [START_REF] Haradhvala | Mutational Strand Asymmetries in Cancer Genomes Reveal Mechanisms of DNA Damage and Repair[END_REF].

Especially APOBEC-associated mutations (also called APOBEC mutation signatures) in humans are represented in up to 15% of all sequenced tumors and contribute to 50% of all mutations in many tumors. APOBEC-associated mutations preferentially occur on the laggingstrand template during DNA replication, and are also highly associated with mismatch repair and transcription-coupled damage repair in cancer [START_REF] Cortez | APOBEC3A is a prominent cytidine deaminase in breast cancer[END_REF][START_REF] Hoopes | APOBEC3A and APOBEC3B Preferentially Deaminate the Lagging Strand Template during DNA Replication[END_REF][START_REF] Shi | APOBEC-mediated Mutagenesis as a Likely Cause of FGFR3 S249C Mutation Over-representation in Bladder Cancer[END_REF]M. J. Shi et al., 2020;[START_REF] Mas-Ponte | DNA mismatch repair promotes APOBEC3-mediated diffuse hypermutation in human cancers[END_REF]. Furthermore, N6methyladenosine (m6A) modifications have been considered as one of the most prevalent internal modifications in mammalian mRNAs and the abnormal m6A modification caused by m6A modulators, e.g., methyltransferase-like 3 (METTL3), is a common feature of various tumors (Y. Shi et al., 2020;[START_REF] Wang | N6-methyladenosine METTL3 promotes the breast cancer progression[END_REF][START_REF] Huang | m6A Modification in Coding and Non-coding RNAs: Roles and Therapeutic Implications in Cancer[END_REF]. Evidence has shown that METTL3 and m6A could promote homologous recombination-mediated repair of doublestrand breaks (DSBs) by modulating DNA-RNA Hybrid (R-loops) accumulation (Zhang et al., 2020). Importantly, R-loops have been recently shown, by others and us, to be frequently accumulated at transcription termination sites of actively transcribed genes displaying high HO . CC-BY-NC-ND 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted January 13, 2022. ; https://doi.org/10.1101/2022.01.12.476022 doi: bioRxiv preprint TRCs [START_REF] Promonet | Topoisomerase 1 prevents replication stress at R-loop-enriched transcription termination sites[END_REF][START_REF] Liu | Topoisomerase I prevents transcription-replication conflicts at transcription termination sites[END_REF]. Therefore, systematically unveiling the genome-wide DNA replication panorama is essential for human health.

Despite its importance, to date, there is no published available tool to analyze RFD data and to determine the replication initiation and termination positions genome-wide, although several methods have been previously described for OK-seq data analysis, such as using the Hidden Markov Model (HMM) to analyze human OK-seq data [START_REF] Petryk | Replication landscape of the human genome[END_REF] or the origin efficiency metric (OEM) to analyze yeast OK-seq data [START_REF] Mcguffee | Quantitative, genome-wide analysis of eukaryotic replication initiation and termination[END_REF]. It is, therefore, important to have a uniform framework of OK-seq data (and related data) analyses. Here, we developed a bioinformatics toolkit, called OKseqHMM, to directly obtain the high-resolution RFD profile genome-wide. Besides the fork direction, the toolkit also deciphers the information of replication initiation/termination zones using an algorithm based on HMM, calculates the OEM to visualize the transition of RFD profile at multiple scales, and finally generates the average metagene profiles and heatmaps to provide RFD/OEM distributions along the regions of interest (Fig. 1). We have gathered a large number of published available OK-seq data (13 in total) from S. cerevisiae, mouse and human cells, and successfully obtained the highresolutive (~1 kb for mouse and human cells and ~50 bp for yeast) RFD profiles and the accurate calling of corresponding replication initiation and termination zones genome-wide.

Implementation and availability

OKseqHMM toolkit is an R package for profiling OK-Seq data to study the genome-wide replication program. This R package contains multi-functions and is served for analyzing OK-Seq data from the original mapping bam file(s) to count matrices, RFD calculation, initiation/termination zone calling and average metagene profiles/heatmaps. The package is available at https://github.com/CL-CHEN-Lab/OK-Seq.

Function OKseqHMM

This function transforms OK-seq data into RFD profiles for a primary visualization (e.g., with the genomic visualization browsers, such as IGV), then, it can accurately identify replication initiation zones (IZs, upward transitions on RFD profile), termination zones (TZs, downward transitions on RFD profile) and also the intermediate states (flat RFD profile) along the genome by using the HMM.

For each window, RFD was computed as follows:

𝑹𝑭𝑫 = 𝑪 -𝑾 𝑪 + 𝑾
where C and W correspond to the number of reads mapped on the Crick and Watson strands, which reveal, respectively, the proportions of rightward-and leftward-moving forks within each window (e.g., 1 kb window was used for OK-seq data of human cells). Since the total amount of replication on both strands should be constant across the genome, we normalized the difference between the two strands by the total read count to account for the variations in read-depth due to copy number, sequence bias and so on. RFD ranges from -1 (100% leftwardmoving forks) to +1 (100% rightward-moving forks), and 0 means equal proportions of leftward-and rightward-moving forks. Data obtained from biological replicates produced RFD profiles that strongly correlated to each other, for HeLa cells, Pearson R=0.92, p<10 -15 (t-test)

and for GM06990 cells, R=0.93, p<10 -15 . Similar correlations were observed between RFD profiles with EdU or EdC labeling [START_REF] Petryk | Replication landscape of the human genome[END_REF]. As in [START_REF] Petryk | Replication landscape of the human genome[END_REF], a four-state HMM was used in OKseqHMM to detect within the RFD profiles the ascending (AS), descending (DS) and flat (FS) segments representing regions of predominant initiation ('Up' state), predominant termination ('Down' state) and constant RFD ('Flat1' and 'Flat2' states) (Fig. 2A). In the HMM segmentation process, the RFD values were computed within 15 kb (for human OK-seq data) sliding windows (by default, stepped by 1 kb across the autosomes). The HMM used the ∆𝑅𝐹𝐷 values between adjacent windows, in

which ∆𝑅𝐹𝐷 ! = "#$ !"# %"#$ ! &
for the window n. By default, windows with <30 reads on both strands were masked. The ∆𝑅𝐹𝐷 values (also between -1 and 1) were divided into five quantiles and the HMM package of R (http://www.r-project.org/) was used to perform the HMM prediction with probabilities of transition and emission, which are manually defined by the training dataset (Fig. 2B). The same transition and emission probabilities used in our previous study [START_REF] Petryk | Replication landscape of the human genome[END_REF] were set as default values and used in all OK-seq data analyses in the current study. The choice of a 15 kb sliding window is based on a compromise between spatial resolution and reproducibility of AS detection among biological replicates.

Finally, the efficiency of the detected AS was estimated as:

∆𝑹𝑭𝑫 𝒔𝒆𝒈𝒎𝒆𝒏𝒕 = 𝑹𝑭𝑫 𝒆𝒏𝒅 -𝑹𝑭𝑫 𝒔𝒕𝒂𝒓𝒕 𝟐
where RFDstart and RFDend correspond, respectively, to the RFD values computed in 5 kb windows around the left and right extremities of each segment.

Function OKseqOEM

For a further investigation of origin efficiency (i.e., ∆𝑅𝐹𝐷), we provide here a second function to visualize it directly at multiple scales.

As defined in the previous publication for yeast OK-seq data analysis 23 , the density of Okazaki are referred to as origin-firing efficiency; and the degree of termination at each position can be measured from 0 to -1.

Here, we further extend OEM calculation within a fixed window size into multiple-scales to better fit OK-seq data analysis of other organisms, such as human cells.

𝑶𝑬𝑴 𝒊 𝒇𝒐𝒓 𝒍𝒊𝒔𝒕[𝒏] = 4𝑾 𝒊1𝒍𝒊𝒔𝒕[𝒏] -𝑾 𝒊 5 4𝑾 𝒊1𝒍𝒊𝒔𝒕[𝒏] -𝑾 𝒊 5 + (𝑪 𝒊1𝒍𝒊𝒔𝒕[𝒏] -𝑪 𝒊 )
Where list[n] can be defined by users as a list of windows (e.g. [1,10,20,50,100]), i is from 1 to the total length of the data -list [n]. C and W correspond to the number of reads mapped on, respectively, the Crick and Watson strands within corresponding windows.

Using the two bam files of reads within, respectively, Watson and Crick strands generated by the previous OKseqHMM function and the annotation coordinates, the function OKseqOEM can automatically calculate the OEM profiles at a series of defined scales (e.g., from 1 kb to 1

Mb for human cells), which allows us to directly visualize the transition states of replication and also to validate the IZs identified by OKseqHMM then to double-check the size and boundary of IZs.

Average metagene profile/heatmap

To analyze RFD distributions around and/or among the genomic regions of interest, such as the identified IZs or the annotated genes, we developed an additional module for the metadata analysis. With the gene coordinates (or IZs) together with the RFD and/or OEM big wiggle files generated from OKseqHMM and/or OKseqOEM functions, we can easily obtain the corresponding profiles/heatmaps by using the computeMatrix and plotProfile/plotHeatmap functions of deepTools (https://deeptools.readthedocs.io/en/develop/index.html) [START_REF] Ramírez | deepTools2: a next generation web server for deep-sequencing data analysis[END_REF] via defining the genomic distances of interest for the upstream and downstream borders.

Performance

Genome-wide replication fork directionality and origin detection in yeast

To evaluate the performance, we first applied our tool to the available yeast OK-seq data [START_REF] Hennion | FORK-seq: Replication landscape of the Saccharomyces cerevisiae genome by nanopore sequencing[END_REF]. OKseqHMM was successfully applied to the yeast OK-seq data to generate the RFD profile at a fine resolution (50 bp), the OEM profiles at different scales (from 50 bp to 25 kb) and a precise IZ/Origin calling (Fig. 3A). About 350 IZs were identified by OKseqHMM, which range from 0.5 kb to 5.5 kb with an average length of 1.5 kb (Fig. 3B, Table 1). To check the accuracy of IZ calling results, we compared OK-seq IZs with the known yeast origins, i.e. autonomously replicating sequence (ARS) from OriDB 2.1.0 [START_REF] Siow | OriDB, the DNA replication origin database updated and extended[END_REF], and up to 80% of our detected IZs were found at ≤ 1 kb distance from a known ARS (Fig. 3C andD).

Genome-wide replication fork directionality and initiation zone detection for different human cell lines

We then further applied the OKseqHMM to analyze the OK-seq data of human cells. In addition to the published OK-seq data of HeLa MRL2 cells [START_REF] Petryk | Replication landscape of the human genome[END_REF], we also generated new additional OK-seq data from HeLa S3 cells, a widely used Encode Tier 2 cell line. The RFD profiles of the two HeLa cell lines are very similar (R=0.86, p<10 -15 ), suggesting that a similar replication program and IZ positions are used (Fig. 4A). About 10,000 IZs have been identified in each HeLa cell line (Table 1), two-third of which are common between the two cell lines (Fig. 4B). The conservation of IZs is even higher in the early-replicating regions, with 80% of early IZs being shared between the two HeLa cell lines (Fig. 4B). A very striking difference of human RFD data compared to those of yeast is that, instead of a sharp 1 kb upward transition of RFD at fixed yeast origins, the size of upward transition of RFD, therefore the IZ length, of human cells is around 10-50 kb (average ~30 kb, ~20-folds larger than the IZ of yeast) (Fig. 4A, Table 1). The heatmaps of OEM profiles computed around IZs at different scales show the strongest positive signals at the corresponding scales, i.e., 10 kb scale for the small IZs (<10 kb), 20 kb scale for the IZs of mid-size (20-50 kb) and 50 or 100 kb for the large IZs (>50 kb), respectively (Fig. 4C), confirming that RFD transition is associated with the detected IZ length. This further supports the difference between the yeast and human OKseq pattern and the accuracy of IZ detection obtained by OKseqHMM.

Replication initiation has been previously reported to be enriched within intergenic regions between active genes [START_REF] Petryk | Replication landscape of the human genome[END_REF]. To demonstrate how our toolkit can help in the analysis of the association between DNA replication and gene transcription, we analyzed the average profiles and the corresponding heatmap of expression level (RNA-seq and GRO-seq) for all detected IZs sorted by their length and confirm that gene transcription presents immediately surround the IZs while with a much lower level within IZs (Fig. 5A). To further compare the distribution of RFD and gene transcription, we calculated the average RFD profile and the corresponding heatmap around TSSs (transcription start sites) and TTSs (transcription termination site) of 16,336 active genes (RPKM > 0) in HeLa cells with an extension ± 50 kb upstream or downstream (Fig. 5B). This clearly indicates a frequent replication initiation (upward transition of RFD) at both regions upstream of TSS and downstream of TTS, which leads to a co-direction between replication and transcription at TSS while a higher head-on TRC at TTS, in agreement with previous publications [START_REF] Petryk | Replication landscape of the human genome[END_REF][START_REF] Promonet | Topoisomerase 1 prevents replication stress at R-loop-enriched transcription termination sites[END_REF]. In addition to OK-seq data of HeLa cell lines, we have gathered and reanalyzed with OKseqHMM, the OK-seq data from previous publications for a large amount of human cell lines of different cell types (Wu et al., 2018;[START_REF] Petryk | Replication landscape of the human genome[END_REF], such as fibroblast (IMR90), lymphoblastoid (GM06990) and lymphoma (Raji, BL79, IARC385), leiomyosarcoma and leukemia (IB118, TLSE19, K562), and erythroblast (TF1) (Table 1, Fig. 6). OKseqHMM generated high-quality cell-type-specific RFD profiles and robust IZ calling for all data analyzed. The sizes of IZs in different cell types are within the same range (average size between 26 to 36 kb), demonstrating that it is a common feature of human cells. [START_REF] Li | DNA polymerase α interacts with H3-H4 and facilitates the transfer of parental histones to lagging strands[END_REF]. Mean and standard error bands are shown for both data, while the standard error bands of OK-seq data are too narrow to be seen.

Extend OKseqHMM to analyze the RFD profiles from other sequencing data

In addition to OK-seq, the OKseqHMM toolkit can be applied to calculate RFD profiles from the sequencing data obtained with other related techniques. As a demonstration, we further extended our toolkit to analyze the published eSPAN [START_REF] Li | DNA polymerase α interacts with H3-H4 and facilitates the transfer of parental histones to lagging strands[END_REF] and TrAEL-seq [START_REF] Kara | Genome-wide analysis of DNA replication and DNA double-strand breaks using TrAEL-seq[END_REF] data. The RFD data computed from the yeast TrAEL-seq data are very close to those obtained by OK-seq (Fig. 7A, R=0.93, p<10 -15 ), and the RFD profile obtained by TrAEL-seq even shows a higher quality with less local noise than the OK-seq RFD profile, which is probably because the TrAEL-seq data used in the analysis contain about two-fold more reads compared with the available OK-seq data. The comparison between the TrAEL-seq IZs, OKseq IZs and yeast ARSs showed that up to 96% of detected IZs from TrAEL-seq were found within 1 kb distance from a known ARS and around 76% of OK-seq IZs associated with ARSs were also detected by TrAEL-seq (Fig. 7B). We also successfully applied OKseqHMM to the OK-seq and eSPAN data of mouse embryonic stem cells (mESC). However, due to the lower amount of reads for the available dataset of eSPAN data, although we used a larger window size (e.g., 10 kb smoothing window instead of 1 kb window) we still got too noisy RFD profiles to perform a robust IZ calling. Nevertheless, we still obtained a mean RFD profile similar to those of OK-seq around the IZs identified in the mESC OK-seq data (Fig. 7C).

Discussion

Genome-wide replication fork directionality data have become an important key in understanding numerous biological processes, such as transcription-replication conflicts, replication-associated mutagenesis, replication couple epigenetic maintenance, etc. Here, we present OKseqHMM, a comprehensive R package, to analyze OK-seq data from various cell types and species to generate and visualize high-resolutive RFD and OEM profiles along the genomes, as well as generate the average profiles/heatmaps on the regions/genes of interest.

The toolkit also allows accurate detection of replication initiation/termination zones with an HMM algorithm. To our knowledge, this is the first bioinformatics tool available to date to handle and analyze the RFD data obtained from various techniques.

We successfully applied OKseqHMM to a large amount of available OK-seq data from different species, including yeast, mouse cells as well as numerous normal and cancer human cell lines (Table 1). This provides an important resource for large research communities, who are interested in studying DNA replication programs, transcription-replication conflicts, replication-associated chromatin organization, replication-associated mutations, genome instability and cancer genomics, among others. Importantly, in addition to OK-seq, more and more new techniques have been developed to study DNA replication and are also able to provide the replication fork direction information. These include the methods like eSPAN and SCAR-seq performing stranded sequencing of BrdU or EdU labeled nascent replicated DNA associated with specific histone modifications, or like TrAEL-seq and GLOE-seq based on the single-stranded end presented on specific replicative templates. Here, we demonstrated that

OKseqHMM can be applied to analyze data obtained by both kinds of techniques, i.e., eSPAN and TrAEL-seq, and obtained high-quality results (Fig. 7). Notably, techniques like TrAELseq, which do not need to incorporate labels and need fewer cells to generate a high-quality RFD profile compared to OK-seq, will provide a good alternative to study DNA replication and genomic instability in different cell types within various stress conditions.

It should be noted that the initiation parameters, such as the transition and emission probabilities, are defined based on the OK-seq datasets of human cells. Although we have shown in the current study that they are quite robust and can be also applied to OK-seq data of yeast (Fig. 1) and mouse cells (Fig. 7, Table 1) to obtain satisfactory results, they might need to be adjusted based on the sequencing-depth and data quality of other datasets, in order to have an optimal IZ/TZ calling. In the future, with technical improvement, we might be able to further extend the OKseqHMM to study the extrinsic (cell-to-cell) or intrinsic (homolog-to-homolog) variability of DNA replication, if we can further extend the relative techniques to obtain data at the single-cell level and/or in an allele-specific manner as recently achieved for the replication timing study [START_REF] Dileep | Single-cell replication profiling to measure stochastic variation in mammalian replication timing[END_REF][START_REF] Gnan | Kronos scRT: a uniform framework for single-cell replication timing analysis[END_REF].

via targeting Bcl-2. Gene, 722, 144076.

Figure 1 .

 1 Figure 1. Schematic presentation of data analysis pipeline of OKseqHMM toolkit. The raw sequencing data can be pre-processed into aligned files by corresponding bioinformatics tools indicated in blue (left panel). The middle panel shows the major functions of the OKseqHMM toolkit. The first function of OKseqHMM checks the input aligned bam files to determine whether they are single-or paired-end sequencing data, then automatically splits the reads into Watson and Crick strands and calculates the replication fork directionality (RFD). By default, the calculation is performed within 1 kb adjacent windows (recommended for human cells) and then smoothed into 15 kb sliding windows with 1 kb step. These parameters can be easily adjusted based on the nature of the data. Different replication features, i.e., initiation zones (IZ), two intermediate states and termination zones, are predicted based on an HMM algorithm (See Implementation for detail). The second function (OKseqOEM) uses the reads on Watson and Crick strands to generate origin efficiency metrics (OEMs) at multiple scales to visualize the RFD transition. And the last function allows users to generate an average metagene profile and heatmap to analyze distributions of RFD and OEM around the genes/regions of interest. Results can be visualized in the genomic visualization browsers (such as IGV) as shown in the right panel.

Figure 2 .

 2 Figure 2. Schematic presentation of HMM algorithm for initiation and termination zone detection. (A) A 4state HMM model used in the segmentation process: Up, regions of predominant initiation (IZ); Down, regions of predominant termination (TZ); Flat1 and Flat2, two intermediate transition states. (B) Transition and emission probabilities. (C) and (D) Examples of RFD profile of HeLa cells from chromosome 1 with the corresponding IZs, TZs and 2 Flat states identified by OKseqHMM. Each point on the RFD profile gives the RFD value calculated within each 1 kb adjacent window, and the windows with positive and negative RFD values are shown in red and blue, respectively.

(

  fragments on the Watson and Crick strands are compared within 4 fixed-size sliding bins, which are strand-specific 10 kb quadrant values to calculate an Origin Efficiency Metric (OEM), computed as 𝑂𝐸𝑀 = WL and WR measure, respectively, the read density in the left and right quadrants on the Watson strand, while CL and CR refer to the density on the Crick strand), ranging from -1 to 1 for each base in the genome. Maximal values in the OEM scores represent replication origins, while the minimal ones are considered as regions of replication termination. In addition, the different amplitudes of positive OEMs (from 0 to 1)

Figure 3 .

 3 Figure 3. Analysis of OK-seq data by OKseqHMM in Yeast. (A) Yeast RFD profile was calculated at 50 bp resolution with the corresponding IZs identified by OKseqHMM, which are highly correlated with the confirmed origins from OriDB (Siow et al., 2012). RFD profile as in Fig. 2C, but with 50 bp resolution. Below, the OEM profiles calculated from 50 bp to 25 kb scales, and the windows with positive and negative OEM values are shown in red and blue, respectively. (B) Length distribution of detected OK-seq IZs. (C) The density profile shows the distribution (in red) of distances between an IZ detected by OkseqHMM and the closest confirmed origin from OriDB, which is much closer compared with the random simulation control (in black). (D) Venn diagram showing the overlap between OK-seq IZs and published origins from OriDB, in which overlap means that the closest distance between each other is less than 1 kb. In case of overlap, the numbers of OK-seq IZ are shown together with the corresponding numbers of OriDB origins indicated in brackets.

Figure 4 .

 4 Figure 4. Analysis of OK-seq data of HeLa cells by OKseqHMM. (A) Replication timing profile obtained by Repli-seq (Chen et al., 2010), RFD profiles and corresponding detected IZs for published HeLa MRL2 OK-seq data (Petryk et al., 2016) and OK-seq data of HeLa S3 cells generated in the current study, the OEM profiles of HeLa S3 cells from 1 kb to 1 Mb scales, and the transcription data provided by GRO-seq and RNA-seq along a ~4 Mb region on chromosome 1. (B) Venn diagrams showing that two-third of Ok-seq IZs matched between the two HeLa cell lines and the overlap goes up to 80 % for the early IZs (with replication timing S50 < 0.4). (C) Mean OEM profiles and heatmaps of OEM (heatmap color scale is indicated on the right) around the HeLa S3 IZ centers at indicated scales (i.e., 10, 20, 50 and 100 kb) sorted by the length of detected IZs.

Figure 5 .

 5 Figure 5. OKseqHMM reveals the coordination between DNA replication and gene transcription. (A) Mean profiles and heatmaps of RNA-seq and GRO-seq around the HeLa S3 OK-seq IZ centers. (B) Mean profile and heatmap of HeLa S3 RFD between TSS (transcription start site) and TTS (transcription termination site) of active genes with an extension of +/-50 kb. The heatmap color scales are indicated in each panel.

Figure 6 .

 6 Figure 6. Genome-wide RFD profiles of different human cell lines show the cell-type-specific replication program. Cell-type-specific RFD profiles and the corresponding detected IZs for indicated human cell lines, IMR90, TF1, K562, TLSE19, GM06990, Raji and BL79, respectively.

Figure 7 .

 7 Figure 7. Genome-wide RFD profiles obtained from TrAEL-seq and eSPAN data. (A) RFD profiles and the corresponding IZs in 50 bp bin size of OK-seq and TrAEL-seq data of yeast (Kara et al., 2021). The known origins (ARSs) are downloaded from OriDB. (B) Venn diagram showing the overlap between OK-seq IZs, TrAEL-seq IZs and published origins (ARSs) from OriDB, in which overlap means that the closest distance between each other is less than 1 kb. (C) Metagene average RFD profiles computed from OK-seq of mouse embryonic stem cells (mESC) and H4K20me2 eSPAN data of MCM2-2A mutant cells[START_REF] Li | DNA polymerase α interacts with H3-H4 and facilitates the transfer of parental histones to lagging strands[END_REF]. Mean and standard error

  Wu,X. et al. (2018) Developmental and cancer-associated plasticity of DNA replication preferentially targets GC-poor, lowly expressed and late-replicating regions. Nucleic Acids Res., 46, 10157-10172. Zhang,C. et al. (2020) METTL3 and N6-Methyladenosine Promote Homologous Recombination-Mediated Repair of DSBs by Modulating DNA-RNA Hybrid Accumulation. Mol. Cell, 79, 425-442.e7.

Table 1 .

 1 All OK-seq data analyzed by OKseqHMM. If data of biological replicates are available, the profiles obtained with the combined data are used in the figures, and only the segments (i.e. IZs and TZs) reproducibly identified in both biological replicates were retained.

				Initiation zones	Termination zones
	Name	Cell type	Replicate		Size (kb)		Size (kb)	Accession number (Reference)
				Number	Mean ± SD	Number	Mean ± SD
	BL79	Burkitt's lymphoma		7798	29±18	7791	211±244	ENA: PRJEB25180 (Wu et al., 2018)
	GM06990	Lymphoblastoid cells	2*	5684	33±19	5715	182±166	SRA: SRP065949 (Petryk et al., 2016)
	HeLa MRL2	Epithelial cell of adenocarcinoma	2*	9836	31±18	9441	141±144	SRA: SRP065949 (Petryk et al., 2016)
	HeLa S3	Epithelial cell of adenocarcinoma		9089	32±19	9084	223±245	(Current study)
	IARC385	B lymphocytes from Burkitt's lymphoma		4465	36±19	4455	125±164	ENA: PRJEB25180 (Wu et al., 2018)
	IB118	Leiomyosarcoma		3645	26±16	3640	428±440	ENA: PRJEB25180 (Wu et al., 2018)
	IMR90	Fibroblast		12482	26±17	12468	151±147	ENA: PRJEB25180 (Wu et al., 2018)
	K562	Late-stage chronic myeloid leukemia		6982	28±15	6967	136±158	ENA: PRJEB25180 (Wu et al., 2018)
	mESC E14	Mouse embryonic stem cells		3370	27±14	3347	483±554	GEO: GSE142996 (Li et al., 2020)
	Raji	Burkitt's lymphoma		8096	29±16	8080	143±135	ENA: PRJEB25180 (Wu et al., 2018)
	TF1	BCR-ABL negative cell line from erythroblast		8377	27±17	8371	196±193	ENA: PRJEB25180 (Wu et al., 2018)
	TLSE19	Leiomyosarcoma		10500	27±17	10492	146±144	ENA: PRJEB25180 (Wu et al., 2018)
	Yeast	S. cerevisiae	2*	348	1.5±0.7	787	14±13	ENA: PRJEB36782 (Hennion et al., 2020)

*

Acknowledgements and funding information. This work was supported by the YPI program of I. Curie; the ATIP-Avenir program from Centre national de la recherche scientifique (CNRS) and Plan Cancer [grant number ATIP/AVENIR: N°18CT014-00]; the Agence Nationale pour la Recherche (ANR) [grant number ReDeFINe -19-CE12-0016-02, TELOCHROM -19-CE12-0020-02]; and Institut National du Cancer (INCa) [grant number PLBIO19-076]. YL thanks ANR for providing her Ph.D. fellowship. The authors would like to thank Joseph Josephides for critical reading, and to acknowledge the high-throughput sequencing facility of Institute for Integrative Biology of the Cell (I2BC) for its sequencing expertise. We declare no

Data Availability. The bioinformatics tool and all processing data underlying this article are available at the GitHub page of the team https://github.com/CL-CHEN-Lab/OK-Seq. And the raw sequencing data of OK-seq are available with the corresponding accession numbers indicated in Table 1. The OK-seq data of HeLa S3 cells were generated as described in [START_REF] Petryk | Replication landscape of the human genome[END_REF] and are available at GEO. The RNA-seq and GRO-seq data of HeLa cells are from [START_REF] Promonet | Topoisomerase 1 prevents replication stress at R-loop-enriched transcription termination sites[END_REF] and [START_REF] Andersson | Nuclear stability and transcriptional directionality separate functionally distinct RNA species[END_REF], respectively. The known yeast origins (ARSs) are downloaded from OriDB (http://cerevisiae.oridb.org/) [START_REF] Siow | OriDB, the DNA replication origin database updated and extended[END_REF].

competing interests and no potential conflicts of interest.


Mehdi Jouhari

Francis Chopin

Mohamed El Houicha

Jean-François Ghienne

Karel Schulmann

Jitka Míková

Michel Corsini

Mohamed El Houicha

U-Pb detrital zircon geochronology and source provenance in the Moroccan Meseta (Variscan belt): A perspective from the Rehamna massif

Keywords: Detrital zircon geochronology, Rehamna, Moroccan Variscan belt, Paleozoic

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

This paper presents a detrital zircon U-Pb geochronology study of samples from the Rehamna massif in the Moroccan Variscan belt [START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF], which aims at better constraining its lithostratigraphy and the evolution of source provenance during the Lower Paleozoic development. The hardness and resistance to chemical weathering make zircon one of the most durable minerals. Therefore, it can survive wide range of geological process (e.g. metamorphism, erosion and sedimentary transport) and preserve primary source signatures [START_REF] Hanchar | Zircon zonation patterns as revealed by cathodoluminescence and backscattered electron images: Implications for interpretation of complex crustal histories[END_REF]. In particular, detrital zircon grains in sedimentary rocks can survive multiple sedimentary cycles involving the recycling of old sedimentary, igneous and metamorphic rocks and younger primary crystalline sources [START_REF] Thomas | Detrital-zircon geochronology and sedimentary provenance[END_REF]. Age distribution of detrital zircon in sedimentary rocks can be used for correlation of similar stratigraphic units, provenance analysis via comparison with known sources and determination of maximum depositional ages in the absence of datable contemporaneous volcanic rocks [START_REF] Reiners | Geochronology and thermochronology[END_REF].

The paleogeographic evolution of the Moroccan Variscan belt (Fig. 1) is mostly based on geological mapping, tectonic studies and lithological correlations [START_REF] Hoepffner | The Moroccan Hercynides[END_REF][START_REF] Simancas | The tectonic frame of the Variscan-Alleghanian orogen in Southern Europe and Northern Africa[END_REF][START_REF] Simancas | Tectonic relationships of Southwest Iberia with the allochthons of Northwest Iberia and the Moroccan Variscides[END_REF][START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF][START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF]. These ideas were [3] recently tested by a new generation of detrital zircon studies mostly used to precise paleogeography and source areas in northern Gondwana (e.g. [START_REF] Pereira | New U-Pb zircon dating of late Neoproterozoic magmatism in Western Meseta (Morocco)[END_REF]Perez-Caseres et al., 2017;Letsch et al., 2017;[START_REF] El Houicha | Recycling of the Proterozoic crystalline basement in the Coastal Block (Moroccan Meseta): New insights for understanding the geodynamic evolution of the northern peri-Gondwanan realm[END_REF][START_REF] Ghienne | The impact of the end-Ordovician glaciation on sediment routing systems: A case study from the Meseta (northern Morocco)[END_REF][START_REF] Perez | Reconstructing Extensional Basin Architecture and Provenance in the Marrakech High Atlas of Morocco: Implications for Rift Basins and Inversion Tectonics[END_REF][START_REF] Accotto | Mixed and recycled detrital zircons in the Paleozoic rocks of the Eastern Moroccan Meseta: Paleogeographic inferences[END_REF][START_REF] Accotto | Tectonic Evolution of the Eastern Moroccan Meseta: From Late Devonian Forearc Sedimentation to Early Carboniferous Collision of an Avalonian Promontory[END_REF]Accotto et al., , 2021aAccotto et al., ,b, 2022 and references therein) and references therein). So far examined detrital zircon populations concern sedimentary rocks deposited in various late Neoproterozoic to Paleozoic basins preserved in the Moroccan Meseta. Their origin generally refers to recycling of the West African Craton (Archean domains, Eburnean Orogeny), the Pan-African Trans-Saharan belt, the Saharan Metacraton to the east, and their overlying Paleozoic strata (e.g., [START_REF] Accotto | Mixed and recycled detrital zircons in the Paleozoic rocks of the Eastern Moroccan Meseta: Paleogeographic inferences[END_REF]. Some studies focused on complete stratigraphic sections allow evaluating evolution of source areas through time or possible impact of remote geodynamic processes (e.g. [START_REF] Gärtner | The zircon evidence of temporally changing sediment transport-the NW Gondwana margin during Cambrian to Devonian time (Aoucert and Smara areas, Moroccan Sahara)[END_REF]Accotto et al., 2021). In this context, the Rehamna massif still remains understudied.

To fill the gap, a detrital zircon U-Pb geochronology study of a deformed Paleozoic sequence was performed in the southern part of the Rehamna massif of the Moroccan Meseta (Fig. 1). The present study aims at determining maximum depositional age of the strata and to give new insights into sedimentary provenance in an area where structural and stratigraphic uncertainties persist to date. Finally, this study allows to contribute to better understand the tectonic evolution of the northern Gondwana during late Paleozoic times.

Geological setting

The Paleozoic evolution of the Moroccan Variscan belt (Fig. 1) is well integrated in the evolution of the northern margin of Gondwana from its early Paleozoic extensional basins until its amalgamation with Laurussia during the late Carboniferous-early Permian [START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF][START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF][START_REF] Martínez Catalán | The Mid-Variscan Allochthon: Keys from correlation, partial retro deformation and plate-tectonic reconstruction to unlock the geometry of a non-cylindrical belt[END_REF]. Three main Paleozoic events can be distinguished: [START_REF] Vermeesch | How many grains are needed for a provenance study?[END_REF] Since the late Cambrian formation of the Rheic ocean, the northern edge of the West African Craton (WAC) evolved in a passive margin environment (Stampfli and Borel,[4] 2002; [START_REF] Ribeiro | Geodynamic evolution of the SW Europe Variscides[END_REF][START_REF] Oriolo | Early Paleozoic accretionary orogens along the Western Gondwana margin[END_REF]. [START_REF] Avigad | Coupled U-Pb-Hf of detrital zircons of Cambrian sandstones from Morocco and Sardinia: Implications for provenance and Precambrian crustal evolution of North Africa[END_REF] showed that until the late Devonian, the sedimentation reflects the decreasing contribution of detrital material derived from local sources. (2) Subsequent formation of the Paleotethys ocean was associated to late Devonian-early Carboniferous intracontinental rifting [START_REF] Stampfli | A plate tectonic model for the Paleozoic and Mesozoic constrained by dynamic plate boundaries and restored synthetic oceanic isochrons[END_REF][START_REF] Frizon De Lamotte | Evidence for Late Devonian vertical movements and extensional deformation in northern Africa and Arabia: Integration in the geodynamics of the Devonian world[END_REF]. (3) This event was followed by the late Carboniferous-early Permian intracontinental Variscan Orogeny triggered by the collision between Gondwana with the European Variscan belt and Laurussia, which ultimately lead to the formation of Pangea [START_REF] Hoepffner | The Moroccan Hercynides[END_REF][START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF][START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF].

The Moroccan Meseta shows rare Precambrian basement units (e.g. [START_REF] Pereira | Evidence of a Paleoproterozoic basement in the Moroccan Variscan belt (Rehamna Massif, Western Meseta)[END_REF][START_REF] El Houicha | Recycling of the Proterozoic crystalline basement in the Coastal Block (Moroccan Meseta): New insights for understanding the geodynamic evolution of the northern peri-Gondwanan realm[END_REF] unconformably overlain by thick Paleozoic successions deformed during the late Paleozoic [START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF][START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF]. The Variscan belt in Morocco can be divided into two main tectonic domains (Fig. 1; e.g., Piqué et al., 1993a;[START_REF] Hoepffner | The Moroccan Hercynides[END_REF][START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF]: ( 1 [START_REF] Hoepffner | La tectonique hercynienne dans l'Est du Maroc[END_REF][START_REF] Accotto | Mixed and recycled detrital zircons in the Paleozoic rocks of the Eastern Moroccan Meseta: Paleogeographic inferences[END_REF][START_REF] Accotto | Tectonic Evolution of the Eastern Moroccan Meseta: From Late Devonian Forearc Sedimentation to Early Carboniferous Collision of an Avalonian Promontory[END_REF]. The Western Meseta [START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF] consists of almost undeformed Cambrian-Devonian strata of its western "Coastal Block" and a central highly deformed and metamorphosed zone composed of crystalline basement covered by Cambrian-Carboniferous strata [START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF][START_REF] Delchini | The Geological Evolution of the Variscan Jebilet Massif, Morocco, Inferred From New Structural and Geochronological Analyses[END_REF] outcropping in the Western High Atlas, Jebilet, Rehamna and Central massifs. The exotic Sehoul Block is located to the northwest and is either assumed to have an Avalonia [START_REF] El Hassani | La Zone de Rabat-Tiflet: bordure nord de la Chaîne Calédono-Hercynienne du Maroc[END_REF]Piqué et al., 1993b;[START_REF] Tahiri | Geochronological data on the Rabat-Tiflet granitoids: Their bearing on the tectonics of the Moroccan Variscides[END_REF] or Meguma affinity (Michard et al.,[5] 2010) and was accreted to the Moroccan Meseta during the early Carboniferous [START_REF] Tahiri | Geochronological data on the Rabat-Tiflet granitoids: Their bearing on the tectonics of the Moroccan Variscides[END_REF].

Geology of the Rehamna massif

The Paleozoic Rehamna massif is situated midway south of Casablanca and north of Marrakech in the Western Meseta (Fig. 1). Its southern part is sometimes called the Rehamna metamorphic dome in contrast to the northern part (not shown in Fig. 2), which reveals solely the presence of comprehensive Cambrian to Permian sedimentary sequences. In this work we deal exclusively with the metamorphic southern part, which is also characterized by Precambrian inliers. The Rehamna massif is traditionally divided into three Western, Central and Eastern Rehamna zones separated by the Western Meseta Shear Zone (WMSZ, Fig. 1) or Median Fault and the Ouled Zednes Fault (Fig. 2; [START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF][START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF]. Based on the intensity of deformation and the degree of metamorphism, the southern part is further subdivided into upper and lower metamorphic units (Fig. 3) according to the intensity of deformation and the degree of metamorphism [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF][START_REF] Razin | Carte géologique du Maroc au 1/50 000, feuille de Jebel Kharrou, Mémoire explicatif[END_REF][START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF][START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF].

Lithostratigraphy

In contrast to the northern part of the Rehamna massif, where the age of the temporal framework of the succession is well understood (e.g. [START_REF] Hollard | Le massif Paléozoïque des Rehamna (Maroc): Stratigraphie, Tectonique et Petrogenese d'un segment de la Chaine Varisque[END_REF][START_REF] El Kamel | Sédimentologie, magmatisme pré-orogénique et structuration du Paléozoïque des Rehamna et d'Ouled Abbou (Meseta occidentale, Maroc)[END_REF], the lithostratigraphy of its southern part is still controversial (Fig. 3; [START_REF] Destombes | Le massif Paléozoïque des Rehamna (Maroc): Stratigraphie, Tectonique et Petrogenese d'un segment de la Chaine Varisque[END_REF][START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF][START_REF] Razin | Carte géologique du Maroc au 1/50 000, feuille de Jebel Kharrou, Mémoire explicatif[END_REF][START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF]. In the lower metamorphic units of the Central and Eastern Rehamna, the amphibolite facies metamorphism obliterated fossiliferous records, and the tentative stratigraphy is solely based on facies analogies with fossiliferous sedimentary sequences cropping out to the north (Hoepffner, 1974;[START_REF] Jenny | Contribution à la géologie structurale des Rehamna (Meseta marocaine méridionale). Le matériel paléozoïque et son évolution hercynienne dans le centre du massif[END_REF].

The Paleoproterozoic basement of the Western Rehamna is characterized by metamorphosed porphyroids and rhyolites [START_REF] Pereira | Evidence of a Paleoproterozoic basement in the Moroccan Variscan belt (Rehamna Massif, Western Meseta)[END_REF] covered by low-grade early- [6] middle Cambrian sediments (Fig. 3; [START_REF] Corsini | Découverte de la série basale du Cambrien et de son substratum dans les Rehamna centraux, haut-fond au Cambrien (Meseta marocaine). Prémices de l'orogenèse hercynienne[END_REF]. The sequence starts with early Cambrian arkoses and limestones showing a gradually increasing detrital component (e.g. Lalla Mouchaa Formation;[START_REF] El Houicha | Recycling of the Proterozoic crystalline basement in the Coastal Block (Moroccan Meseta): New insights for understanding the geodynamic evolution of the northern peri-Gondwanan realm[END_REF]. This formation is conformably overlain by the "Paradoxides Shale Formation" composed of siltstone, greywacke and sandstone preserving a middle Cambrian fauna [START_REF] Guézou | Note sur la structure du môle côtier mésétien dans l'ouest des Rehamna (Maroc hercynien)[END_REF][START_REF] Corsini | Découverte de la série basale du Cambrien et de son substratum dans les Rehamna centraux, haut-fond au Cambrien (Meseta marocaine). Prémices de l'orogenèse hercynienne[END_REF].

To the east of the WMSZ, a similar succession is present in the lower metamorphic unit of the Central Rehamna. Here, an early Ediacaran (593±8 Ma; [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF] rhyolitic basement crops out in the core of the Sidi Ali dome (Fig. 2; [START_REF] Corsini | Découverte de la série basale du Cambrien et de son substratum dans les Rehamna centraux, haut-fond au Cambrien (Meseta marocaine). Prémices de l'orogenèse hercynienne[END_REF] and is conformably overlain by early Cambrian terrigenous clastic sediments (conglomerate, displaying granite and quartzite pebbles, interbedded with arkosic sandstone and mudstone; Fig. 3). A limestone sequence (including calcschists interbeds) marks the transition to early Cambrian siliciclastic sedimentation [START_REF] Corsini | Découverte de la série basale du Cambrien et de son substratum dans les Rehamna centraux, haut-fond au Cambrien (Meseta marocaine). Prémices de l'orogenèse hercynienne[END_REF]. The succession continues with earlymiddle Cambrian metapelitic strata (Draa El Kebir Formation) unconformably overlain by the Kef El Mouneb metaconglomerate assumingly Devonian (e.g. Hoepffner et al., 1972). The age of the quartzites from the Skhour Formation of the upper metamorphic allochthonous unit is still controversial and was considered to be Ordovician by [START_REF] Destombes | Le massif Paléozoïque des Rehamna (Maroc): Stratigraphie, Tectonique et Petrogenese d'un segment de la Chaine Varisque[END_REF]. However, based on structural and sedimentological arguments, [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF] and [START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF] alternatively proposed a possible Upper Devonian age for this formation (Fig. 3). The structurally lower metamorphic unit of the Eastern Rehamna is divided into the Devonian (Ouled Hassine-Ouled Zednes, OH-OZ) and Carboniferous (Lalla Tittaf) formations (Fig. 3).

The former formation corresponds to a succession of pelitic beds interbedded with poorlypreserved Crinoidea-bearing limestones [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF] and limy/sandy conglomerates or quartzites, generally attributed to the Devonian (Hoepffner, 1974). Preserved Eifelian-Frasnian corals were recently reported from this formation (marble of the Ouled Zednes area; El Houicha et al., 2019). The Lala Tittaf Formation consists of pelites intercalated with andesites, gabbros [7] and rhyolites (Hoepffner, 1974;[START_REF] Cornée | Étude lithostratigraphique et tectonométamorphique des Rehamna sudorientales. Plissements et nappes. Contribution à la connaissance de la chaîne hercynienne en Meseta marocaine[END_REF]. Considered first as Carboniferous, [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF] proposed a Paleoproterozoic age based on the dating of zircon grains from a gabbroic intrusion (2136 ± 17 Ma). However, this was challenged by a new U-Pb zircon data which provide early Carboniferous ages (339.5 ± 0.89 Ma) for the intraformational rhyolite sills [START_REF] Aït Lahna | The Lalla Tittaf Formation (Rehamna, Morocco): Paleoproterozoic or Paleozoic age?[END_REF].

The upper unit to the north is formed by the allochthonous and very low grade Jbel Kharrou and Koudiat El Adam Formations (Fig. 2) which are made of Middle Ordovician to Silurian shale and sandstone [START_REF] Destombes | Le massif Paléozoïque des Rehamna (Maroc): Stratigraphie, Tectonique et Petrogenese d'un segment de la Chaine Varisque[END_REF]. To the southeast and south of the Ouled Ouggad Fault, the upper unit is represented by various virtually unmetamorphosed rocks (Fig. 2). The latter are composed of Ordovician strata overlain by Devonian (-Tournaisian?) sandstone and pelite and by the middle (?) late Visean to Serpukhovian sandstone and limestone including a sandy flysch facies typifying in the Dalaat unit (Fig. 3;[START_REF] Cornée | Étude lithostratigraphique et tectonométamorphique des Rehamna sudorientales. Plissements et nappes. Contribution à la connaissance de la chaîne hercynienne en Meseta marocaine[END_REF][START_REF] Razin | Carte géologique du Maroc au 1/50 000, feuille de Jebel Kharrou, Mémoire explicatif[END_REF].

Magmatism and metamorphism

Two main metamorphic events were described in the southern part of the Rehamna massif. The first one corresponds to a prograde Barrovian metamorphism (D1-M1) reaching the staurolite (kyanite) stability field with estimated peak P-T at conditions 6-9 kbar and 410-560 °C. This metamorphism affects the lower metamorphic unit during the late Carboniferous-early Permian and was interpreted as a result of tectonic burial [START_REF] Hoepffner | Le massif Paléozoïque des Rehamna (Maroc): Stratigraphie, Tectonique et Petrogenese d'un segment de la Chaine Varisque[END_REF][START_REF] Aghzer | Détachements et tectonique extensive dans le massif hercynien des Rehamna (Maroc)[END_REF][START_REF] El Mahi | L'évolution structurale et métamorphique de la zone hercynienne des Rehamna centraux (Maroc)[END_REF][START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF][START_REF] Wernert | Tectonometamorphic evolution of an intracontinental orogeny inferred from P-T-t-d paths of the metapelites from the Rehamna massif (Morocco)[END_REF]. The D1-M1 event was followed by a greenschist facies retrograde metamorphism and D2 upright E-W folding of the entire sequence [START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF]. This late tectonic event was probably contemporaneous with the intrusion of the Ras El Abiod granitoid [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF]. The M1 Barrovian and M2 retrograde metamorphic events were dated at 310-295 Ma and 295-285 Ma, respectively using 40 Ar/ 39 Ar (micas, amphiboles) and U-Pb (monazite) techniques [START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF] [8] [START_REF] Wernert | Tectonometamorphic evolution of an intracontinental orogeny inferred from P-T-t-d paths of the metapelites from the Rehamna massif (Morocco)[END_REF]. Late to Post orogenic felsic magmatism (Piqué, 1979;[START_REF] Lagarde | Les plutons granitiques tardi carbonifères marqueurs de la déformation crustale. L'exemple des granitoïdes de la méseta marocaine[END_REF] is represented by the intrusions of the Sebt Brikiyine granite (Fig. 2) dated at ca. 275 Ma (e.g. Chopin et al., 2014 and reference therein).

Methods -zircon separation, imaging and measurements

Two to four kilograms of rock per sample were processed at the laboratories of Czech Geological Survey in Prague. The zircon grains were separated by sieving (0.2 and 0.1 mm), heavy liquid (acetylene tetrabromide and methylene iodide) and magnetic separation. Then, they were selected randomly from the zircon concentrates by hand-picking under binocular microscope. The grains were mounted in epoxy discs, ground and polished.

Cathodoluminescence (CL) and secondary electrons imaging were performed at the SEM Laboratory of the Czech Geological Survey (FEG-SEM TESCAN MIRA3 GMU working at 15 kV and ~22 mm working distance) in order to visualize internal zircon texture and zoning patterns of individual grains. We tried to get around 150 analyzes per sample which allowed to obtain representative and statistically significant age populations after discordant analyses are removed [START_REF] Vermeesch | How many grains are needed for a provenance study?[END_REF].

The U-Pb dating of zircon grains was performed at the laboratories of the Czech Geological Survey. Measurements were performed on an Analyte Excite 193 nm excimer laserablation system (LA; Photon Machines), equipped with a two-volume HelEx ablation cell, in tandem with an Agilent 7900x ICPMS (Agilent Technologies Inc., Santa Clara, USA). Samples were ablated in He atmosphere (0.8 l min -1 ) at a pulse repetition rate of 10 Hz using a spot size of 25 μm and laser fluence of 4.71 Jcm -2 . Each measurement consisted of 20 seconds of blank acquisition followed by ablation of the sample for a further 40 seconds of signal collection at masses 202,204,206,207,208,232 and 238 using the SEM detector, with one point per mass peak and the respective dwell times of 10, 10, 15, 30, 20, 10 and 15 milliseconds per mass (total sweep time of 0.134 seconds). Instrumental drift was monitored by repeated measurements of [9] the 91500 (ca. 1063 Ma, [START_REF] Wiedenbeck | Three natural zircon standards for U-Th-Pb, Lu-Hf, trace element and REE analyses[END_REF] reference zircon. Data deconvolution using Iolite software followed the method described by Paton et al. (2010), including an 'on peak' gas blank subtraction followed by correction for laser-induced elemental fractionation (LIEF) by comparison with the behaviours of the 91500 reference zircon [START_REF] Wiedenbeck | Three natural zircon standards for U-Th-Pb, Lu-Hf, trace element and REE analyses[END_REF] which yielded a concordia age at 1062.6 ± 1.6 Ma. No common Pb correction was applied in this study. In addition, zircon reference samples GJ-1 (~609 Ma; [START_REF] Jackson | The application of laser ablation-inductively coupled plasma-mass spectrometry to in situ U-Pb zircon geochronology[END_REF] and Plešovice (~337 Ma; [START_REF] Sláma | Plesovice zircon -A new natural reference material for U-Pb and Hf isotopic microanalysis[END_REF] were analysed periodically during this study and yielded concordia ages at 609.3 ±1.5 Ma and 338.5 ±1.0 Ma (2σ), respectively.

The analytical data were plotted using the IsoplotR toolbox [START_REF] Vermeesch | IsoplotR: A free and open toolbox for geochronology[END_REF]. Data reduction followed the recommendation of [START_REF] Spencer | Strategies towards statistically robust interpretations of in situ U-Pb zircon geochronology[END_REF] whereby data is treated as concordant if the 2σ ellipse overlaps with the concordia, and the preferred ages are single zircon concordia ages following the recommendations of [START_REF] Nemchin | Discordance of the U-Pb system in detrital zircons: Implication for provenance studies of sedimentary rocks[END_REF] and [START_REF] Zimmermann | Maximising data and precision from detrital zircon U-Pb analysis by LA-ICPMS: The use of core-rim ages and the singleanalysis concordia age[END_REF]. For all concordant ages younger than 1000 Ma the 206 Pb/ 238 U ages were used, whereas older ages are based on the obtained 207 Pb/ 206 Pb ratios [START_REF] Gehrels | Detrital Zircon U-Pb Geochronology: Current Methods and New Opportunities[END_REF]. We kept U-Pb zircon ages with a degree of concordance situated in between 95% and 105% which corresponds to the number of analyzed zircon needed to get sufficient representative data.

4.Results

In this study, we collected early Cambrian to Mississippian (meta)sandstone samples from eight sites for U-Pb detrital zircon geochronology. From each location, we took several pieces of rock from the main sedimentary strata to ensure the representativeness of the sampling. The geographic, geological and stratigraphic characteristics of each sampling sites are shown in figures 2 and 3, and table 1. Note that no zircon grain were found after mineral separation of (meta)clastic rocks samples from the Draa El Kebir, Kef El Mouneb and Lalla Tittaf Formations. Therefore, only five age spectra are presented in the following. subordinate clusters at 1700-2200 Ma (18%) and another one at ca. 1 Ga (1%). A small peak at ca. 580 Ma is also present (Fig. 8a).

Sample R257-OH

R257-OH was collected in quartzite lenses of the eastern extremity of the OH-OZ Formation (near the Souk el Had; Fig. 2, Fig. 3). The formation can be either Lower or Middle to Upper Devonian [START_REF] El Houicha | Discovery of Eifelian-Frasnian corals in metamorphic rocks from the Rehamna massif (Moroccan Meseta): biostratigraphic and paleogeographic implications. 11th Colloquium 3MA « Magmatism, Metamorphism and associated Mineralization[END_REF] depending on the interpretation of the stratigraphic polarity (e.g. Hoepffner, 1974;[START_REF] Cornée | Étude lithostratigraphique et tectonométamorphique des Rehamna sudorientales. Plissements et nappes. Contribution à la connaissance de la chaîne hercynienne en Meseta marocaine[END_REF][START_REF] Razin | Carte géologique du Maroc au 1/50 000, feuille de Jebel Kharrou, Mémoire explicatif[END_REF]. These weakly metamorphosed and deformed impure quartzite (Fig. 4d) is composed of recrystallized, elongate grains of quartz with lobate boundaries and irregular undulatory extinction and minor recrystallized muscovite (Fig. 5c). Large detrital zircon grains (100-150 µm) have a long prismatic shape whereas smaller zircon population (50-75 µm in size) are sub-rounded. Most of the large zircons show a sector zoning but some are characterized by prismatic zonation patterns (Fig. 6c). Several rounded grains are homogeneous and weakly luminescent. A set of 138 concordant ages has been obtained from 161 analyses (Fig. 8a). They show a main Cryogenian-Ediacaran population (540-850 Ma; 79%) with a peak at ca. 632 Ma and a smaller Paleoproterozoic cluster (1700-2200 Ma; 17%) with two peaks at ca. 1764 and 2071 Ma. Distributed Mesoproterozoic and Archean ages are also present (Fig. 8b).

Sample R255-OZ

The sample R255-OZ was collected next to the fossiliferous Eifelian-Frasnian limestones [START_REF] El Houicha | Discovery of Eifelian-Frasnian corals in metamorphic rocks from the Rehamna massif (Moroccan Meseta): biostratigraphic and paleogeographic implications. 11th Colloquium 3MA « Magmatism, Metamorphism and associated Mineralization[END_REF] of the OH-OZ Formation along the Ouled-Zednes Fault zone (Fig. 2, Fig. 3). The site shows a narrow belt of steeply SE-dipping quartzite, shale and limestone of the OH-OZ Formation surrounded by Skhour Formation rocks in the northwest and Koudiat El Adam rocks in the SE (Fig. 4g,e). The sampled rock is a light-grey low-grade quartzite (Fig. 4f) showing incipient bulging type recrystallization at margins of detrital quartz grains (100 to 200 µm in size) and presence of small sericite grains (Fig. 5d). The age of this rock is probably [12] Upper Devonian (post-Frasnian) to Dinantian in stratigraphic continuity with the Middle Devonian Ouled Zednes Formation of the lower metamorphic unit [START_REF] Jenny | Contribution à la géologie structurale des Rehamna (Meseta marocaine méridionale). Le matériel paléozoïque et son évolution hercynienne dans le centre du massif[END_REF]Razin, 2003;[START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF]. However, [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF] did not exclude that it might represent a tectonic slice from the upper metamorphic unit, either from the western Cambro-Ordovician or Upper Devonian Skhour Formation, or from the eastern Ordovician Jbel Kharrou-Koudiat El Adam nappe (Fig. 2). The zircon grains reveal a complex texture in their CL images with variable luminescence and patterns. Most of the extracted detrital zircon from this sample is about 60-100 µm in size, showing oscillatory or sector zoned with irregular or straight grain boundaries. Less frequently, low-to high luminescent homogenous elongated grains 120-140 µm in size are present (Fig. 6d). One hundred twenty-two concordant ages from 178 analyses have been obtained from this sample (Fig. 8c). They show multimodal age spectra ranging between 310-540 Ma (13%), 540-850 Ma (52%), 850-1200 Ma (15%) and 1500-2200 Ma (15%). A younger population shows subordinate peaks at 477 and 559 Ma. Nevertheless, the youngest detrital zircon ages are close to 350 Ma (Fig. 6d, Fig. 8c).

Sample R256-DT

The sample R256-DT was collected from the Dalaat unit (member 2 in the Dalaat Formation s.s.; [START_REF] Razin | Carte géologique du Maroc au 1/50 000, feuille de Jebel Kharrou, Mémoire explicatif[END_REF]Fig. 2;Fig. 3). This unit corresponds to middle (?)-late Visean basal fine to medium grained arkosic sandstones alternating with limestone (Fig. 4h,i). The sandstone is composed of quartz and feldspar grains, 150-400 µm in size, some detrital muscovite, chloritized biotite and interstitial calcite (Fig. 5e). The zircon grains show either subhedral shapes and oscillatory and sector zoning or sub-rounded shapes and homogeneous weakly or highly luminescent patterns (Fig. 6e). The studied zircon grain size falls primarily into the 120-170 µm range and less frequently into the 60-100 µm range. A total of 182 analyses were carried out, of which 136 displayed concordant ages (Fig. 8d). The sample R256-DT is dominated by a Cryogenian-Ediacaran population (68%) ranging between 540 and 850 Ma with a peak at 633.9 ± 0.5 Ma and three subordinate populations at 850-1200 Ma (12%), 1500-2200 Ma (13%) and 2200-2500 Ma (13%) suggesting potential peaks around 991 Ma, 1943 Ma and2446 Ma, respectively (Fig. 8b).

Most of the analyses from all the samples are characterized by Th/U ratios > 0.1 (ranging between 0.1 and 1.7) suggesting a primary magmatic origin [START_REF] Rubatto | Zircon trace element geochemistry: Partitioning with garnet and the link between U-Pb ages and metamorphism[END_REF]. However, some zircon rims show Th/U ratios lower than 0.1, in particular those in the 600-850 Ma range, which is compatible with a metamorphic remobilization.

Discussion

Evolution of detrital zircon sources in time

South of the Sebt Brikiyine granite, [START_REF] Pereira | Evidence of a Paleoproterozoic basement in the Moroccan Variscan belt (Rehamna Massif, Western Meseta)[END_REF] described an Eburnean basement dated at 2050.6 ± 3.0 Ma (SHRIMP U-Th-Pb, rhyolitic porphyry, sample SBR-21 in Fig. 2 and Fig. 3). The basement is unconformably overlain by clastic (pelites, arkoses) and carbonate (limestones and calcschists) strata similar to the early Cambrian succession outcropping in the northern Lalla Mouchaa Formation and in the Sidi Ali dome in the lower metamorphic unit of the Central Rehamna [START_REF] Corsini | Découverte de la série basale du Cambrien et de son substratum dans les Rehamna centraux, haut-fond au Cambrien (Meseta marocaine). Prémices de l'orogenèse hercynienne[END_REF]. The sample R16-1A (Fig. 2), an arkosic sandstone directly overlying the Paleoproterozoic basement, only show zircon grains with Paleoproterozoic ages (mean age 2048.3 ± 1.4 Ma). This cluster is similar to the one obtained from early Cambrian microbreccia from the Lalla Mouchaa Formation (Sample NB1; Fig. 2, Fig. 3) with detrital zircons at ca. 2047 ± 12 Ma [START_REF] El Houicha | Recycling of the Proterozoic crystalline basement in the Coastal Block (Moroccan Meseta): New insights for understanding the geodynamic evolution of the northern peri-Gondwanan realm[END_REF]. Both samples (R16-1A, this study, and NB1, El Houicha et al., 2018) are characterized by angular feldspar and euhedral zircons, the latter with concentric zoning suggesting a magmatic origin and a low roundness indicative for a very short transport distance and thus, low-order recycling of sediments. Correlation with published ages shows that the detrital zircon likely derived from first-cycle erosion of the local Eburnean basement [START_REF] Pereira | Evidence of a Paleoproterozoic basement in the Moroccan Variscan belt (Rehamna Massif, Western Meseta)[END_REF](Pereira et al., ) indicating proximal al., 2021b) ) show that such an age spectrum does not allow differentiating Ordovician from Devonian strata and that similar zircon populations are identified in both the Ordovician of some of the autochthonous strata (Ezzhelligha area, EZ on the Fig. 1) and allochthonous strata (Azrou area, AZR on the Fig. 1). Nevertheless, because sample SK1-Z has affinities with both Ordovician samples from the Ezzhelligha area (Accotto et al., 2021b) and the late Cambrian (El Hank) from the nearby Coastal Block (sample 15DL12, IMF on the Fig. 1) [START_REF] Letsch | A missing link in the peri-Gondwanan terrane collage: The Precambrian basement of the Moroccan Meseta and its lower Paleozoic cover[END_REF], we favor an Ordovician age for the SK1-Z rock unit, and not favor an Upper Devonian age as suggested by [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF] and [START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF]. Whatever the stratigraphic attribution of the samples might be, SK1-Z and R257-OH and similar ones from the Meseta show that the southern-derived West Gondwana source (i.e. dominant Cadomian/Pan-African and Eburnean source) essentially pertained from the Cambrian-Ordovician up to the Devonian (Fig. 9b; [START_REF] Gärtner | The zircon evidence of temporally changing sediment transport-the NW Gondwana margin during Cambrian to Devonian time (Aoucert and Smara areas, Moroccan Sahara)[END_REF]Accotto et al., 2021b). The minor presence of Meseoproterozoic zircons, also described elsewhere in the Meseta might reflect a minor lateral source from the Sahara Metacraton (see however discussion in Accotto et al., 2021b).

The sample R255-OZ collected in the OH-OZ Formation shows a major change with two grains indicating a maximum depositional age as young as early Carboniferous (350.6 ±10 Ma and 352.5 ±7 Ma, i.e Tournaisian), i.e. significantly younger than the Devonian strata hosting the sample R257-OH. This is in full agreement with a distinctive age spectrum, showing Ordovician to Lower Devonian subordinate peaks and a less prominent Pan-African population.

As in R256-DT (see next paragraph), Mesoproterozoic zircon grains are present but more distributed within the Mesoproterozoic instead of a peak at ~1 Ga. It suggests a different sourcing, confirmed by the morphology of zircons which are mostly euhedral and poorly rounded -in contrast to detrital zircons from the underlying Cambrian to Devonian strata-, indicating a relatively short transportation path. If the Cambrian to Lower Ordovician zircon grains can be derived from nearby sedimentary sources [START_REF] Letsch | A missing link in the peri-Gondwanan terrane collage: The Precambrian basement of the Moroccan Meseta and its lower Paleozoic cover[END_REF]Accotto et al., [16] 2021b), again indicating the recycling of the Meseta early Paleozoic sedimentary cover, occurrences of late Paleozoic and Mesoproterozoic zircon grains suggests in addition great changes in the distribution of the watersheds. In particular, the latter are well known to include renewed volcanic activity at the regional scale [START_REF] Kharbouch | Le volcanisme dévono-dinantien du Massif central et de la Meseta orientale[END_REF][START_REF] Moreno | Age and depositional environment of the Draa Sfar massive sulfide deposit, Morocco[END_REF][START_REF] Aït Lahna | The Lalla Tittaf Formation (Rehamna, Morocco): Paleoproterozoic or Paleozoic age?[END_REF]. The source of Mesoproterozoic zircons might reflect the presence of an allochthonous unit positioned to the north and docked in the course of the Variscan collision as suggested from studies in the Eastern Meseta (e.g. Accotto et al., 2020b). However, the limited number of "Variscan" late Paleozoic zircons grains in the Rehamna massif, which are largely more represented in the Eastern Meseta, makes this sourcing from the northern watersheds unlikely. On the other hand, the sediment transport scenario from the WAC cannot be excluded.

The samples R255-OZ and R256-DT show some increased Mesoproterozoic component which is well known from the sedimentary cover of the western edge of the WAC (Taoudeni basin; [START_REF] Bradley | Synthesis of geological, structural, and geochronologic data (Phase V, Deliverable 53[END_REF][START_REF] Gärtner | The zircon evidence of temporally changing sediment transport-the NW Gondwana margin during Cambrian to Devonian time (Aoucert and Smara areas, Moroccan Sahara)[END_REF]. Moreover, detrital zircon data indicates a dispersal of Mesoproterozoic detritus since Cambrian times (e.g. [START_REF] Gärtner | The zircon evidence of temporally changing sediment transport-the NW Gondwana margin during Cambrian to Devonian time (Aoucert and Smara areas, Moroccan Sahara)[END_REF]Fig 10b).

Derivations from either the southwest, running parallel to the continental margin and the coevally forming orogen, or from the northwest (Avalonia or European Variscan Orogen; i.e., orthogonal to the Variscan collision) might be postulated.

The middle (?)-late Visean sample R256-DT collected from the basal siliciclastic Dalaat [START_REF] Ghienne | The impact of the end-Ordovician glaciation on sediment routing systems: A case study from the Meseta (northern Morocco)[END_REF]Accotto et al., 2020a;2021). A combination of the two schemes, a western Gondwana source combined with the recycling of older strata, is in fact likely.

Geodynamic implications

This study indicates a modification through time of detrital sources feeding the continental shelf and the sedimentary basins in the current Rehamna area, which is compatible with the geodynamic model recently proposed for the Variscan belt [START_REF] Martínez Catalán | The Mid-Variscan Allochthon: Keys from correlation, partial retro deformation and plate-tectonic reconstruction to unlock the geometry of a non-cylindrical belt[END_REF]. The reconstructions based on existing paleomagnetic and geological data and models of [START_REF] Domeier | Plate tectonics in the late Paleozoic[END_REF] and [START_REF] Matthews | Global plate boundary evolution and kinematics since the late Paleozoic[END_REF] allow to visualize evolution of zircon sources as suggested in this work, in complementing and agreeing with the study of Accotto et al. (2021a). Since the Ordovician, the northern Gondwana has evolved in a passive margin setting, suggesting that the detrital zircon populations deposited in Ordovician to Lower Devonian basins were entirely governed by recycling Gondwanan sources from the south and SE, represented either by remote exposed basements, and/or their sedimentary covers (Fig. 10a). However, the situation radically changed owing to two separate late Devonian to early Carboniferous events. The first one is related to the propagation of the Paleotethys ocean [START_REF] Stampfli | Remnants of the paleotethys oceanic suture-zone in the western tethyan area[END_REF][START_REF] Stampfli | The formation of Pangea[END_REF] into the interior of Gondwana thereby forming the network of basins typical for Eastern and Western Meseta [START_REF] Cózar | Lower-middle Viséan transgressive carbonates in Morocco: Palaeobiogeographic insights[END_REF] and oceanic realm further east [START_REF] Edel | Late Palaeozoic palaeomagnetic and tectonic constraints for amalgamation of Pangea supercontinent in the European Variscan belt[END_REF]Fig. 10b), which is typified in the Rehamna massif by the formation of the Devonian to early Carboniferous basins and associated magmatism (Hoepffner, 1974;1982;[START_REF] Aït Lahna | The Lalla Tittaf Formation (Rehamna, Morocco): Paleoproterozoic or Paleozoic age?[END_REF]. At the same time, a ribbon fragment of Gondwana called the Mid-Variscan Autochthon was obducted by giant allochthonous nappe complexes [START_REF] Martínez Catalán | The Mid-Variscan Allochthon: Keys from correlation, partial retro deformation and plate-tectonic reconstruction to unlock the geometry of a non-cylindrical belt[END_REF] involving both oceanic and continental rocks (Fig. 10b). This event occurred simultaneously with possible westward indentation of Avalonian promontory with the Variscan autochthon as discussed by Accotto et al. (2021a). It was recently shown that the Mid-Variscan allochthon also contains sedimentary sequences with abundant [18] Mesoproterozoic zircon signature [START_REF] Lindner | The Dobra Gneiss and the Drosendorf Unit in the southeastern Bohemian Massif, Austria: West Amazonian crust in the heart of Europe[END_REF] previously considered as typical for Avalonia [START_REF] Mazur | Single zircon U-Pb ages and geochemistry of granitoid gneisses from SW Poland: Evidence for an Avalonian affinity of the Brunian microcontinent[END_REF].

Therefore, the opening of Lower Devonian to Carboniferous sedimentary basins within the Moroccan Meseta was coeval with a major convergent event in the north. In this context, Mesoproterozoic zircons grains of the late Paleozoic Rehamna record are thought to have principally originated from the convergent deformation zone, derived from either an Avalonian (e.g. [START_REF] Gärtner | Provenance of detrital zircon from siliciclastic rocks of the Sebkha Gezmayet unit of the Adrar Souttouf Massif (Moroccan Sahara) -Palaeogeographic implications[END_REF] or mid-Allochthonous (e.g. [START_REF] Lindner | The Dobra Gneiss and the Drosendorf Unit in the southeastern Bohemian Massif, Austria: West Amazonian crust in the heart of Europe[END_REF] source, even though the massive peak at 1.3 Ga from the latter could be "diluted" in our samples. Indeed, a contribution through recycling of the sedimentary cover of the WAC cannot be totally excluded.

Anyway, the source of Mesoproterozoic zircons from the Variscan core and/or Avalonian indenter (Accotto et al., 2021a) remain questionable as the early Carboniferous paleogeographic position of both possible source areas remains poorly constrained. To the contrary, early Carboniferous zircon grains (sample R255-OZ) are argued to principally derived from early Carboniferous magmatism reported in the Moroccan Meseta [START_REF] Kharbouch | Le volcanisme dévono-dinantien du Massif central et de la Meseta orientale[END_REF], although

we note that the Mid-Variscan Allochthon was intruded by series of magmatic arcs from 380 to 350 Ma (e.g. [START_REF] Lardeaux | The moldanubian zone in the French Massif Central, Vosges/Schwarzwald and Bohemian Massif revisited: differences and similarities[END_REF][START_REF] Deiller | Eclogite subduction wedge intruded by arc-type magma: The earliest record of Variscan arc in the Bohemian Massif[END_REF] which could have also contributed to the detrital zircon record in the Rehamna basin.

Figure Caption

Fig. 1. Simplified geological map of the Moroccan Variscides (modified after [START_REF] Hoepffner | The Moroccan Hercynides[END_REF][START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF][START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF]. Middle Meseta Fault Zone (MMFZ), Western Meseta Shear Zone (WMSZ), South Meseta Fault (SMF), South Meseta Front (SMFr), South Atlas Fault (SAF), Rabat-Tiflet Fault Zone (RTFZ), Sidi Bettache Bbasin (SBB), and Azrou-Khenifra Basin (AKB). Fig. 2. Simplified Geological map of the Rehamna massif (adapted from [START_REF] Corsini | Découverte de la série basale du Cambrien et de son substratum dans les Rehamna centraux, haut-fond au Cambrien (Meseta marocaine). Prémices de l'orogenèse hercynienne[END_REF][START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF][START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF][START_REF] Pereira | Evidence of a Paleoproterozoic basement in the Moroccan Variscan belt (Rehamna Massif, Western Meseta)[END_REF][START_REF] El Houicha | Recycling of the Proterozoic crystalline basement in the Coastal Block (Moroccan Meseta): New insights for understanding the geodynamic evolution of the northern peri-Gondwanan realm[END_REF]. [START_REF] Michard | Le massif Paléozoïque des Rehamna (Maroc) : Stratigraphie, Tectonique et Petrogenese d'un segment de la Chaine Varisque[END_REF]Hoepffner, 1974;[START_REF] Corsini | Découverte de la série basale du Cambrien et de son substratum dans les Rehamna centraux, haut-fond au Cambrien (Meseta marocaine). Prémices de l'orogenèse hercynienne[END_REF][START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF][START_REF] Chopin | Tectonic evolution of the Rehamna metamorphic dome (Morocco) in the context of the [20] Alleghanian-Variscan orogeny[END_REF][START_REF] Pereira | Evidence of a Paleoproterozoic basement in the Moroccan Variscan belt (Rehamna Massif, Western Meseta)[END_REF][START_REF] El Houicha | Recycling of the Proterozoic crystalline basement in the Coastal Block (Moroccan Meseta): New insights for understanding the geodynamic evolution of the northern peri-Gondwanan realm[END_REF]. Table 1. List of samples with details of the analyses carried out on the zircon grains.

[28] [29] [30]

  ) the southern thick-skinned Anti-Atlas belt showing a fold and thrust belt with Paleozoic cover variably detached from its Precambrian basement, and (2) the Moroccan Meseta domain to the north of the Southern Meseta Fault (SMF) subdivided into two main zones, the Eastern and Western Meseta. The Eastern Meseta (e.g. Middelt, Tazekka, Debdou-Mekkam massifs) is characterized by deformed pre-Carboniferous sediments unconformably overlain by middle Visean to Serpukhovian-Bashkirian deposits

Formation

  is also dominated by Ediacaran to Cryogenian ages but with other subordinate, yet significant signals: distributed Mesoproterozoic zircon grains, a small but definite Stenian-Tonian peak, and a early Cryogenian subpopulation. It may be interpreted as a continuum of the Cambrian-Ordovician source, with addition of subordinate new sources, or, alternatively, as the recycling of Meseta-derived Cambrian-Devonian strata (Fig. 9d) because all of the additional signals are also present, admittedly in various proportion, in Lower Paleozoic strata outcropping in the Meseta, and more specifically in northern and north-eastern Meseta [17]

Fig. 3 .

 3 Fig. 3. Synthetic lithostratigraphic logs of the Rehamna massif subdivided into the Western

Fig. 4 .

 4 Fig. 4. Field photographs of the sample location. (a) A quarry at the southwest of the Sebt

Fig. 5 .

 5 Fig. 5. (a-e) Microphotographs of the samples R16-1A, SK1-Z, R255-OZ, R257-OH and R256-

Fig. 6 .

 6 Fig. 6. Cathodoluminescence images of selected zircon grains.

Fig. 7 .

 7 Fig. 7. (a) Concordia diagrams of the sample R16-1A. Green circles (concordant); Gray circles

Fig. 8 .

 8 Fig. 8. (a) Concordia diagrams, zircon U-Pb Kernel Density Estimation (KDE) diagrams and

Fig. 9 .

 9 Fig. 9. Conceptual tectonic and sediment dispersal evolution of the Rehamna massif during the

Fig. 10 .

 10 Fig. 10. Paleogeographic reconstruction (modified after Martinez Catalan et al., 2021) showing

Fig. 1 .

 1 Fig. 1.

Fig

  Fig. 2.

Fig

  Fig. 3.

Fig

  Fig. 4. 790

Fig

  Fig. 9. 801

Acknowledgments

The authors want to gratefully acknowledge Nikol Novotná (SEM laboratory) and the technician from the Czech Geological Survey involved in the preparation of the samples. This work was supported by the TelluS Program of CNRS/INSU and institutional grants from the Czech Geological Survey (grant number 310270 and DKRVO) and from the CNRS IPGS UMR7516 (now ITES UMR 7063). The University Chouaïb Doukkali and Martin Simon are thanked for their support during the field work. Fabien Humbert (Southern University of Science and Technology, Shenzhen) is thanked for the processing of photographs.

Sample R16-1A

The sample R16-1A was collected at the southern border of the Sebt Brikiyine granitic intrusion (Fig. 2). It is an early Cambrian arkosic metasandstone overlying the Paleoproterozoic basement [START_REF] Pereira | Evidence of a Paleoproterozoic basement in the Moroccan Variscan belt (Rehamna Massif, Western Meseta)[END_REF] with an essentially preserved unconformable stratigraphic contact, which was, however, possibly tectonically reactivated (Fig. 3, Fig. 4a). This deformed arkosic metasandstone (Fig. 4b) is composed of equigranular fine-grained quartz grains with well equilibrated straight boundaries ranging from 60 to 500 µm in size. Feldspar and muscovite grains show incipient recrystallization (Fig. 5a). Polygonal zircon grains with angular shape are 120-200 µm grains size while sub-rounded grains of irregular shape are smaller (<50 µm).

Most of zircon grains show homogeneous cores surrounded by a brighter oscillatory domain (CL-light) and darker rim (CL-dark; Fig. 6a). The grains are commonly fractured. One hundred seventeen concordant ages from 230 analyses were obtained (Fig. 7a), showing a unimodal Paleoproterozoic age spectrum (Fig. 7b) peaking at ca. 2000-2100 Ma (mean age 2048.3 ± 1.4 Ma).

Sample SK1-Z

The sample SK1-Z was collected from a quartzitic level surrounded by low grade schist of the Ordovician-Devonian (?) Skhour Formation (Fig. 2, Fig. 3, Fig. 4c). This highly deformed quartzite is characterized by well equilibrated quartz foam texture with subequant quartz grains 150-200 µm in size with straight boundaries and an excellent triple point network. (Fig. 5b).

The zircon population consists of elongated, sub-euhedral to euhedral grains 140-180 µm in size and few more rounded smaller grains 60-100 µm in size. The analyzed zircon grains are characterized by a prismatic and sector zoning with some fractures and inclusions (Fig. 6b).

Ninety-four concordant ages from 123 analyses were obtained (Fig. 8a). Cryogenian-Ediacaran ages (540-850 Ma) constitutes 73% of all analyses with a peak at ca. 626 Ma. There are very watersheds collecting the signal from the nearby rocks (Fig. 9a). At regional scale, the detrital zircons from the basal early Cambrian sediments overlying the Meseta basement reflect exclusive, local recycling of the heterogeneous Precambrian basement. Indeed, detrital zircons ages from the early Cambrian microbreccia at Lalla Mouchaa and (meta)arkose from south of Sebt Brikiyine granite are perfectly similar to the age of the underlying rhyolitic basement (2050.6 ± 3.0 Ma, sample SBR-21 in [START_REF] Pereira | Evidence of a Paleoproterozoic basement in the Moroccan Variscan belt (Rehamna Massif, Western Meseta)[END_REF]. Similarly, a single Cadomian/Pan-African cluster is found in detrital zircon from the early Cambrian microbreccia (ca. 583, sample JD2) covering the El Jadida rhyolitic basement dated at 584.2 ± 4.8 Ma, sample JD-1 in El Houicha et al., 2018). Finally, the Neoproterozoic metavolcanic rocks from the Sidi Ali Dome (593 ± 8 Ma, [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF] is covered by early Cambrian metaconglomerates containing granitic pebbles showing zircon ages in the range 600-640 Ma [START_REF] Pereira | New U-Pb zircon dating of late Neoproterozoic magmatism in Western Meseta (Morocco)[END_REF], whereas detrital zircons ages from the matrix are in the range 548-624 Ma (Sample 15DL11 in [START_REF] Letsch | A missing link in the peri-Gondwanan terrane collage: The Precambrian basement of the Moroccan Meseta and its lower Paleozoic cover[END_REF]. This discrimination between Eburnean and Cadomian/Pan-African basements and sources in the Western Meseta, might reflect an inherited suture zone originating the WMSZ since the early Paleozoic (e.g. [START_REF] Piqué | The Western Meseta Shear Zone, a major and permanent feature of the Hercynian belt of Morocco[END_REF][START_REF] Corsini | Influence de la paléogéographie du Paléozoïque sur les déformations hercyniennes de la Meseta nord-occidentale du Maroc[END_REF]. Samples SK1-Z (upper metamorphic unit) and R257-OH (lower metamorphic unit) show very similar density probability curves, largely dominated by a Cryogenian-Ediacaran population, including a late Ediacaran subpopulation. Additionally, there is a Paleoproterozoic population ranging between 1700-2200 Ma. The main Pan-African peak with moderate Paleoproterozoic ages is typical of West African Cambrian to Devonian age spectra [START_REF] Gärtner | The zircon evidence of temporally changing sediment transport-the NW Gondwana margin during Cambrian to Devonian time (Aoucert and Smara areas, Moroccan Sahara)[END_REF][START_REF] Letsch | A missing link in the peri-Gondwanan terrane collage: The Precambrian basement of the Moroccan Meseta and its lower Paleozoic cover[END_REF][START_REF] Ghienne | The impact of the end-Ordovician glaciation on sediment routing systems: A case study from the Meseta (northern Morocco)[END_REF]Accotto et al., 2021b). Similar age populations might indicate that those two samples originally belong to the same stratigraphic unit: the Devonian stratigraphic age of the R257-OH sample (Fig. 2, Fig. 3) might therefore suggest that the SK1-Z sample is also from a Devonian unit, as proposed by [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF] and [START_REF] Michard | The South-Western Branch of the Variscan Belt: Evidence from Morocco[END_REF]. However, detrital zircon data from the Central massif ( [START_REF] Jenny | Contribution à la géologie structurale des Rehamna (Meseta marocaine méridionale). Le matériel paléozoïque et son évolution hercynienne dans le centre du massif[END_REF]Baudin et al,, 2002;Razin 2003) Cambro-Ordovician or Upper Devonian [START_REF] Baudin | Carte géologique du Maroc au 1/50 000, feuille de Skhour des Rehamna, Mémoire explicatif[END_REF]Razin, 2003) Tournaisian