Mathieu Rousseau-Gueutin
email: mathieu.rousseau-gueutin@inra.fr

Jérôme Morice

Olivier Coriton

Virginie Huteau

Gwenn Trotoux

Sylvie Negre

Cyril Falentin

Gwenaëlle Deniot

Marie-Madeleine Gilet

Frederique Eber

Sylvie Nègre

Gwennaëlle Deniot

Frédérique Eber

Alexandre Pelé

Sonia Vautrin

Joëlle Fourment

Maryse Lodé

Hélène Bergès

Anne-Marie Chèvre
email: anne-marie.chevre@inra.fr

The impact of open pollination on the structural evolutionary dynamics, meiotic behavior, and fertility of resynthesized allotetraploid Brassica napus L

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

Brassica napus allopolyploidy genome dynamic meiotic behavior fertility

Polyploidy, or whole genome duplication (WGD), has played a major role in the evolutionary history of eukaryotes, especially in flowering plants [START_REF] Soltis | The polyploidy revolution then...and now: Stebbins revisited[END_REF]. Recent genomic analyses revealed that all angiosperms have been subjected to at least one round of polyploidy in their evolutionary history, and are thus considered paleopolyploids [START_REF] Garsmeur | Two evolutionarily distinct classes of paleopolyploidy[END_REF]. Despite the numerous studies performed on the origin and evolution of polyploids, it is still unknown how these species stabilized their genome and successfully established.

To determine the immediate impact of polyploidy on genome evolution, studies of either resynthesized or natural and recent allopolyploids (neoallopolyploids) have been performed. Contrasting results were observed, with some allopolyploid species presenting a rather unchanged genome (i.e., Spartina and Gossypium), while others had highly shuffled genomes (i.e., Arabidopsis, Brassica napus, or Tragopogon) [START_REF] Liu | Polyploid formation in cotton is not accompanied by rapid genomic changes[END_REF][START_REF] Baumel | Retrotransposons and genomic stability in populations of the young allopolyploid species Spartina anglica C.E. Hubbard (Poaceae)[END_REF][START_REF] Pontes | Chromosomal locus rearrangements are a rapid response to formation of the allotetraploid Arabidopsis suecica genome[END_REF][START_REF] Mestiri | Newly synthesized wheat allohexaploids display progenitordependent meiotic stability and aneuploidy but structural genomic additivity[END_REF][START_REF] Szadkowski | The first meiosis of resynthesized Brassica napus, a genome blender[END_REF][START_REF] Szadkowski | Polyploid formation pathways have an impact on genetic rearrangements in resynthesized Brassica napus[END_REF]Xiong et al. 2011;[START_REF] Chester | An assessment of karyotype restructuring in the neoallotetraploid Tragopogon miscellus (Asteraceae)[END_REF][START_REF] Chester | Patterns of chromosomal variation in natural populations of the neoallotetraploid Tragopogon mirus (Asteraceae)[END_REF]. In plants, allopolyploidy can be at the origin of a genomic and transcriptomic shock, including gene conversion events, activation of transposable elements, chromatin remodeling, DNA methylation changes, and transcriptional or post-transcriptional changes [START_REF] Madlung | Remodeling of DNA methylation and phenotypic and transcriptional changes in synthetic Arabidopsis allotetraploids[END_REF][START_REF] Kashkush | Transcriptional activation of retrotransposons alters the expression of adjacent genes in wheat[END_REF][START_REF] Salmon | Genetic and epigenetic consequences of recent hybridization and polyploidy in Spartina (Poaceae)[END_REF][START_REF] Salmon | Homoeologous nonreciprocal recombination in polyploid cotton[END_REF][START_REF] Gaeta | Genomic changes in resynthesized Brassica napus and their effect on gene expression and phenotype[END_REF][START_REF] Marmagne | Analysis of gene expression in resynthesized Brassica napus allotetraploids: transcriptional changes do not explain differential protein regulation[END_REF][START_REF] Szadkowski | The first meiosis of resynthesized Brassica napus, a genome blender[END_REF][START_REF] Buggs | Transcriptomic shock generates evolutionary novelty in a newly formed, natural allopolyploid plant[END_REF]. Overall, these various structural and functional changes increase the potential of polyploid species toward functional plasticity and evolutionary novelties, contributing to the phenotypic variability that may enable the species to exploit a wider range of environmental conditions [START_REF] Pires | Flowering time divergence and genomic rearrangements in resynthesized Brassica polyploids (Brassicaceae)[END_REF][START_REF] Fawcett | Plants with double genomes might have had a better chance to survive the Cretaceous-Tertiary extinction event[END_REF][START_REF] Van De Peer | The evolutionary significance of ancient genome duplications[END_REF][START_REF] Jiao | Ancestral polyploidy in seed plants and angiosperms[END_REF].

A major challenge facing the successful establishment of a novel polyploid species reproducing sexually is its fertility. For several species, such as oilseed rape (Brassica napus), or wheat (Triticum aestivum), it has been established that the initial hybridization and its genome doubling are rare events. For that reason, the first few allopolyploid plants will produce seeds by self-fertilization to generate a founder population of the new species. The size and structure of this initial population in the first generations strongly depend on karyotype instability, and the level of aneuploid progenies produced (Zhang et al. 2013a). Recent studies examining resynthesized allopolyploids (allotetraploid B. napus and allohexaploid Triticum aestivum) or neoallopolyploids (allotetraploid Tragopogon species) have shown that novel allopolyploids may present aneuploid progenies (Xiong et al. 2011;[START_REF] Chester | An assessment of karyotype restructuring in the neoallotetraploid Tragopogon miscellus (Asteraceae)[END_REF][START_REF] Chester | Patterns of chromosomal variation in natural populations of the neoallotetraploid Tragopogon mirus (Asteraceae)[END_REF]Zhang et al. 2013b). Two different classes of aneuploids were considered: (i) compensated aneuploids (also called hidden aneuploids), with the expected number of chromosomes, but with a loss and gain of homeologous chromosomes; and (ii) numerical aneuploids, with an unexpected chromosome number. Relatively low levels of numerical aneuploids ($10%) were observed in resynthesized allohexaploid wheat (selfed generation 2) (Zhang et al. 2013a), and in natural populations of neoallotetraploid Tragopogon species ($40 generations old) [START_REF] Chester | An assessment of karyotype restructuring in the neoallotetraploid Tragopogon miscellus (Asteraceae)[END_REF][START_REF] Chester | Patterns of chromosomal variation in natural populations of the neoallotetraploid Tragopogon mirus (Asteraceae)[END_REF] compared to previous observations ($60%) in resynthesized B. napus populations (S5:6) obtained by controlled selfing (Xiong et al. 2011). In this latter species, a lower seed yield was observed in the numerical aneuploids compared with the compensated aneuploids (Xiong et al. 2011).

Oilseed rape (Brassica napus L., AACC, 2n = 38) is an excellent model species to study the immediate structural evolutionary dynamics of allopolyploid species, and the effect of natural selection on its evolution. B. napus is a recent, partially allogamous allopolyploid species that formed $7500 years ago [START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF] after the hybridization and genome doubling of two closely related diploid species, B. rapa (AA, 2n = 20) and B. oleracea (CC, 2n = 18) [START_REF] Nagaharu | Genome analysis in Brassica with special reference to the experimental formation of B. napus and its peculiar mode of fertilization[END_REF][START_REF] Inaba | Phylogenetic analysis of Brassiceae based on the nucleotide sequences of the S-locus related gene, SLR1[END_REF]. This species was presumably selected by humans because a wild population of B. napus was not discovered. To date, it remains unclear how selection allowed genome stabilization. In resynthesized B. napus populations, genetic studies of A01/C01 chromosomes using molecular markers have revealed that the first meiosis acts as a "genome blender" because 50% of the gametes presented homeologous exchanges between these two chromosomes [START_REF] Szadkowski | The first meiosis of resynthesized Brassica napus, a genome blender[END_REF]. The polyploid formation pathway (somatic doubling vs. unreduced gametes) impacts the size and number of translocations [START_REF] Szadkowski | Polyploid formation pathways have an impact on genetic rearrangements in resynthesized Brassica napus[END_REF]. In addition, genetic and cytogenetic analyses of single seed descent obtained by self-fertilizing of independent synthetic B. napus lines deriving from an initial hybridization of B. oleracea "TO1000" and B. rapa "IMB218" revealed that these lines presented numerous nonreciprocal translocation events [START_REF] Gaeta | Genomic changes in resynthesized Brassica napus and their effect on gene expression and phenotype[END_REF]Xiong et al. 2011), and a high level of numerical aneuploids ($60%) (Xiong et al. 2011). This latter study was performed using resynthesized B. napus allopolyploids obtained by controlled self-fertilization, regardless of the genome stability, pollen viability, and seed set, thus permitting the study of highly unstable genomes with a very low reproductive fitness. However, it is most likely that individuals from the first B. napus population were produced via open pollination to allow outcrossing. To date, the effects of open pollination on the structural evolutionary dynamics, meiotic behavior, and fertility of resynthesized allotetraploid B. napus are unknown. In addition, to our knowledge, the size (in megabases) of the structural variations (SVs) in resynthesized B. napus, and their impact on gene content and allelic diversity, have been analyzed in only one synthetic [START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF], precluding the establishment of a general trend. To better understand the role of the reproductive mode on karyotype stability following allotetraploidization in B. napus, and before a putative speciation event, we created two different resynthesized B. napus populations obtained either by controlled selfing or open pollination, with the latter revealing a higher fertility. We then conducted a detailed analysis of the populations obtained from open pollination by investigating their aneuploidy level, meiotic behavior, and genome-wide SVs. From these analyses, we observed that the resynthesized B. napus populations presented a very low level of numerical aneuploids. The genetic (using a 60K Illumina array) and cytogenetic studies of 33 G3 individuals obtained by open pollination allowed the detection of unstable meiosis and extensive shuffling of the genome of these plants. BAC-FISH experiments were subsequently performed to validate four of the SVs identified using the array, allowing us to determine that these deletions of genomic regions resulted from translocations that were fixed at the homozygous stage (via homeologous exchanges, hereafter referred to as HE). We were also able to determine the size of these SVs, and their impacts on gene content and allelic diversity. Overall, our results shed light on the impact of open pollination on the genome dynamics, fertility, and putative successful establishment of novel allopolyploids.

MATERIALS AND METHODS

Plant material

Two different resynthesized B. napus, named "RCC" and "EMZ," were created following the experimental design detailed below and presented in Figure 1. First, two crosses between a B. oleracea and a B. rapa individual were performed. The first cross was performed between a B. oleracea var. alboglabra, a doubled haploid line "RC34" (mother plant) with a B. rapa plant "C1.3" belonging to a fodder variety named "chicon" var. rapifera. Similarly, a cross between two fully homozygous doubled haploid lines was performed: B. oleracea var. botrytis italica "HDEM" (2n = 2x = 18, mother plant) with B. rapa var. trilocularis "Z1" (2n = 2x = 20). All these genotypes were provided by the Biological Resource Center (CRB BrACySol, Rennes-Ploudaniel, France) except Z1 (provided by K. C. Falk, Agriculture and Agri-Food Canada, Ottawa, ON). The resulting amphiploid hybrids (AC, 2n = 2x = 19) were somatically doubled using colchicine [START_REF] Chevre | Cytological studies of tetraploid kale (Brassica oleracea L. ssp acephala) obtained from diploid lines after colchicine treatment[END_REF]), leading to one RCC-S0 (2n = 4x = 38) and EMZ-S0 (2n = 4x = 38)resynthesized B. napus individual from the first and second cross, respectively. The meiotic behavior and chromosome number of the F1 and S0 plants (2n = 38) were evaluated [START_REF] Szadkowski | The first meiosis of resynthesized Brassica napus, a genome blender[END_REF][START_REF] Szadkowski | Polyploid formation pathways have an impact on genetic rearrangements in resynthesized Brassica napus[END_REF]. By selfing (the hand pollination of floral buds before anthesis) one RCC-S0 and one EMZ-S0 plant, we produced RCC-S1 and EMZ-S1 progenies. Thereafter, two strategies (for each genetic background) were used to produce the subsequent generations. In both cases, one seed per plant was used to produce the progeny (Single Seed Descent method). In the first strategy, manual self-fertilization of 11-18 S1 plants was performed to produce S2 and S3 progenies. In the second strategy, 110 S1 plants were grown under the same cage (one cage for RCC and one for EMZ), and in the presence of flies to ensure open pollination, and to facilitate outcrossing between plants. The following year, one seed per 110 G2 plants (from either RCC or EMZ) was grown in a cage (with flies) to obtain the G3 progenies (these plants are hereafter referred as "G3"). As a control of seed fertility for B. napus varieties, 50 plants of the spring variety Drakkar (provided by the CRB BrACySol, Rennes-Ploudaniel, France) were grown each year in another cage, but under the same conditions used for the populations obtained from open pollination.

DNA extraction and SNP genotyping

Genomic DNA of 18 RCC-G3 and 15 EMZ-G3 plants obtained from outcrossing (all presenting 38 chromosomes), as well as of the diploid parents ("RC34," "C1.3," "HDEM" and "Z1"), and the resynthesized S0 B. napus ("EMZ-S0" and "RCC-S0"), were extracted from young leaves using a Qiagen kit. For all these plants, genotyping was performed using the Illumina (http://www.illumina.com/) Brassica 60K Infinium SNP array [START_REF] Clarke | A high-density SNP genotyping array for Brassica napus and its ancestral diploid species based on optimised selection of single-locus markers in the allotetraploid genome[END_REF]. Hybridizations were performed according to the standard procedures provided by the manufacturer for each genomic DNA extracted. Three technical replicates were realized for each diploid parent and each resynthesized S0 B. napus.

Data analyses

The genotyping data from the Illumina Brassica 60K array were visualized using Genome Studio V2011.1 (Illumina, Inc., San Diego, CA), and processed with a manually adapted cluster file. The positions on the B. napus chromosomes of the SNP markers present on the array were obtained by blasting the 52,157 sequence contexts against the B. napus cv. Darmor reference genome sequence assembly (version 4.1: Chalhoub et al. 2014). Only the BLAST hits with a minimum of 90% global overlap and 90% identity were retained. Subsequently, only markers presenting no more than one BLAST hit on each subgenome were kept, enabling SNPs with the potential to hybridize at paralogous regions to be discarded.

Identification of structural rearrangements

To identify putative structural rearrangements in our G3 populations from resynthesized S0 B. napus, the markers that were homozygous and polymorphic (AA vs. BB) between the diploid parents and heterozygous (AB) in the S0 of each population were first determined. In this step, only the markers, for which identical genotype data were obtained for all technical replicates, were considered. Thereafter, a custom python script was applied to identify in each G3 plant, which of the polymorphic markers previously retained had lost one parental allele. The positions on the B. napus chromosomes of markers showing either additivity (both parental alleles) or loss of one parental allele were graphically represented for each G3 B. napus plant using CIRCOS [START_REF] Krzywinski | Circos: an information aesthetic for comparative genomics[END_REF]. By determining whether consecutive markers displayed the loss of one parental allele (from the same parent), we were able to identify large deletions and infer the size and position of each deleted region. The presence of an active centromere in the deleted region was then determined using the positions of the B. napus centromeres [START_REF] Mason | Centromere locations in Brassica A and C genomes revealed through half-tetrad analysis[END_REF]. In addition, we also determined whether deletions were present in the distal region of a chromosome arm: we considered the last 30% of each chromosome arm as the distal region.

Statistical analyses

To determine whether structural rearrangements were statistically prevalent on a specific B. napus subgenome, the proportion of each deleted region in each RCC-G3 and EMZ-G3 individual was estimated based on its subgenome size. Because the data did not follow a normal distribution, the nonparametric Wilcoxon rank-sum test [START_REF] Wilcoxon | Individual comparisons by ranking methods[END_REF] was applied independently for RCC-G3 and EMZ-G3 plants.

Similarly, the proportion of each deleted region in each RCC-G3 and EMZ-G3 was estimated per chromosome to determine whether each B. napus chromosome was equivalently subjected to structural rearrangements. Due to the non-normality of the data, the nonparametric Friedman rank-sum test [START_REF] Friedman | The use of ranks to avoid the assumption of normality implicit in the analysis of variance[END_REF] was applied. For all these analyses, only p-values ,0.05 were considered significant.

Similar statistical analyses were performed by considering the proportion of genes that were deleted in a subgenome and per chromosome.

Validation and identification of the underlying mechanism using BAC-FISH To validate some of the structural rearrangements identified using the SNP array, and to identify the mechanisms involved (i.e., deletion vs. translocation after HE), BAC-FISH combined with GISH-like experiments were performed for plants showing major structural modifications. In these plants, as well as in their diploid parents, chromosome preparations were performed according to procedures detailed in [START_REF] Ksiazczyk | Immediate unidirectional epigenetic reprogramming of NORs occurs independently of rDNA rearrangements in synthetic and natural forms of a polyploid species Brassica napus[END_REF]. For the BAC-FISH experiments, we identified BACs from B. rapa [START_REF] Mun | The first generation of a BAC-based physical map of Brassica rapa[END_REF]) and B. napus [START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF]) libraries that were present in the rearranged regions of some G3 individuals. Some of these BACs had previously been used in similar cytogenetic experiments (Xiong and Pires 2011), whereas others were identified in this study by developing primer pairs specific to the deleted region in a synthetic individual, and screening the B. napus BAC library [START_REF] Chalhoub | Efficient cloning of plant genomes into bacterial artificial chromosome (BAC) libraries with larger and more uniform insert size[END_REF]. All four BACs used in this study hybridized to one pair of chromosomes in each diploid species, and to two homeologous chromosome pairs (four signals) in the B. napus Darmor variety (see the positions of each of these BACs used for the BAC-FISH experiment in Supplemental Material, Table S1). For example, clone KBrB086G22 hybridizes to A02 chromosomes in B. rapa (two signals), C02 chromosomes in B. oleracea (two signals), and both A02 and C02 chromosomes in B. napus (four signals) (Xiong and Pires 2011). Thereafter, to identify the A or C chromosome to which each BAC hybridized in an allotetraploid B. napus, we used, on the same slide, the Bob014O06 BAC clone that hybridizes to C-genome chromosomes in B. napus (GISH-like) [START_REF] Szadkowski | The first meiosis of resynthesized Brassica napus, a genome blender[END_REF]. Each BAC present within a deleted region of a synthetic was labeled by random priming with biotin-14-dUTP (Invitrogen, Life Technologies), whereas the Bob014O06 BAC clone was labeled by random priming with Alexa 488-5-dUTP.

Biotinylated probes were immunodetected using Texas Red avidin DCS (Vector Laboratories), and the signal was amplified with biotinylated anti-avidin D (Vector Laboratories). The chromosomes were mounted and counterstained in Vectashield (Vector Laboratories) containing 2.5 mg/ml 49,6-diamidino-2-phenylindole (DAPI). Fluorescence images were captured using a CoolSnap HQ camera (Photometrics, Tucson, Ariz) on an Axioplan 2 microscope (Zeiss, Oberkochen, Germany), and analyzed using MetaVue (Universal Imaging Corporation, Downingtown, PA).

Meiotic behavior

The flower buds of only G3 plants obtained by open pollination (109 plants) were harvested, enabling the characterization of their meiotic behavior and the establishment of their chromosome number following the protocol described by [START_REF] Suay | Crossover rate between homologous chromosomes and interference are regulated by the addition of specific unpaired chromosomes in Brassica[END_REF]. Per plant, $20 pollen mother cells (PMCs) in metaphase I of meiosis were analyzed.

Assessment of fertility in resynthesized B. napus populations

To assess the fertility of the resynthesized B. napus plants (S0), the number of seeds per pollinated flower was assessed for S1, S2, and S3 plants obtained by controlled manual selfing. For the G1, G2, and G3 plants growing under cages, as well as for the B. napus var. Drakkar, the number of seeds per pollinated flower and per pod was assessed for each plant by counting the number of pods per 50 flowers on the first floral hampers, and the number of seeds per 50 pods.

Data availability

The positions on the B. napus chromosomes of each BAC used in the BAC-FISH experiments are shown in Table S1. The positions on the B. napus chromosomes of the different deletions observed in each synthetic B. napus individual investigated in this study are shown in Table S2. Finally, Table S3 contains the list of the primer pairs designed (names, sequences, and positions on B. napus chromosomes) to identify the presence of the putative genic conversion in some synthetic B. napus individuals.

RESULTS

Impact of open pollination vs. selfing on fertility

The evolution of fertility was evaluated in both "RCC" and "EMZ" progenies, and compared with the Drakkar variety. In all cases, the fertility of the resynthesized B. napus populations, which was assessed based on the number of seeds per 100 flowers, was drastically reduced compared with Drakkar (3-and 10-fold reduction in RCC-G1 and EMZ-G1, respectively) (Figure 2). We also observed that the fertility was significantly higher (t-test, p , 0.05) for plants obtained by open pollination except at the G2 generation for RCC, where no significant difference was observed (t-test, p = 0.815) (Figure 2). In all populations, fertility significantly decreased over subsequent generations (t-test, p , 0.05), except in self-fertilized EMZ plants. In the latter case, the lack of a decrease in fertility was due to the extremely low fertility of these plants since the first generation (Figure 2).

Chromosome number and meiotic behavior of G3 plants obtained via open pollination of resynthesized B. napus

The chromosome number, and meiotic behavior at metaphase I, of 109 RCC-G3 and EMZ-G3 individuals, as well as of the B. napus Drakkar variety, were determined. Most RCC (91%) and EMZ (95%) plants had 38 chromosomes, apart from a few numerical aneuploid plants that consistently presented one or two supernumerary chromosomes (i.e., 39 or 40 chromosomes). Meiosis was evaluated in 15-18 randomly chosen G3 plants obtained from open pollination (per population). Meiosis was affected in these plants, all of which presented 38 chromosomes, with only 30-70% of the PMCs having 19 bivalents, compared to 100% in the B. napus Drakkar variety (Table 1). All G3 individuals showed meiotic instability, and presented multivalents and/or univalents. We observed at least one quadrivalent in 4-48% PMCs of all the plants genotyped (apart from RCC-G3-111), and one univalent (from 1 to 6 per cell) in 5-70% of PMCs.

Genotyping data obtained from the 60K array, and the positions of the SNP markers on the B. napus chromosomes Genotyping data were obtained using the Brassica 60K SNP array for 18 RCC-G3 plants and 15 EMZ-G3 (2n = 4x = 38), as well as for their parental diploid parents ("HDEM," "Z1," "RC," and C1.3), and the two initial allotetraploid resynthesized S0 plants ("RCC-S0" and "EMZ-S0"). From this array, composed of 52,157 SNP markers, 15,058 hybridized in both B. rapa and B. oleracea (in either "HDEM" or "RC34" for B. oleracea, and in either "Z1" or "C1.3" for B. rapa), and thus hybridized to both homeologous chromosomes in our resynthesized B. napus. Using the recently released B. napus cv. Darmor genome sequence [START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF], we identified the positions of all the markers on the B. napus chromosomes. A total of 17,115 markers were specific to the C-subgenome, 12,969 were specific the A-subgenome, and 9504 hybridized to both the A-and C-subgenomes (homeologs).

The positions of these different SNP markers on B. napus chromosomes are presented in Figure S1. Using these stringent filtering BLAST parameters, we observed several BLAST hits on the same subgenome (paralogs) for ,8% of the markers, which were thus discarded from further analyses. Indeed, 534, 573, and 346 heterozygous markers were observed in the "HDEM," "RC34," and "Z1" DH lines, indicating that a low percentage of the SNP markers present on the array may hybridize to paralogous regions, as previously mentioned [START_REF] Clarke | A high-density SNP genotyping array for Brassica napus and its ancestral diploid species based on optimised selection of single-locus markers in the allotetraploid genome[END_REF].

Genotyping and identification of putative structural rearrangements in resynthesized B. napus populations

To identify SVs in each G3 individual from resynthesized B. napus, only the 60K SNP markers that hybridized to both diploid progenitors, and were polymorphic between them, were considered. For such markers, we verified that they were also heterozygous in resynthesized S0 B. napus. In total, 1757 markers were polymorphic between HDEM (B. oleracea) and Z1 (B. rapa), and heterozygous in EMZ-S0. A total of 2405 markers were polymorphic between RC (B. oleracea) and C1.3 (B. rapa), and heterozygous in RCC-S0, and most of these markers could be localized on B. napus chromosomes (Figure S2 and Figure S3). Overall, there was an average of one polymorphic marker every 243 kb for the RCxC1.3 (RCC) combination, and 316 kb for the HDEMxZ1 (EMZ) combination. These markers covered all B. napus chromosomes but with a lower density in the pericentromeric compared with distal regions. The lowest and highest marker densities were observed for the C01 (one every 452 kb), and the A03 (one every 108 kb) chromosomes for the RCC plants. For the EMZ plants, the lowest and highest marker density corresponded to C09 (one every 808 kb) and A08 (one every 153 kb). Using these polymorphic markers alone, we subsequently assessed whether one parental allele was missing from each of the 18 RCC-G3, and the 15 EMZ-G3 individuals obtained via open pollination (an example is provided for EMZ-G3-078 in Figure S4). In the RCC-G3 population, 0 (RCC-G3-075) to 245 SNP markers (RCC-G3-053 and RCC-G3-055) presented a loss of one of the two parental alleles, whereas between 12 (EMZ-G3-096) and 160 SNPs (EMZ-G3-088) demonstrated such a loss in the EMZ population. To identify losses of large genomic regions in an individual, we next determined whether consecutive markers revealed a loss of one parental allele (from the same parent) (Figure 3). Using this method, we identified between one and six large (mean: 3.06), and between zero and seven deleted regions (mean: 3.78) in RCC and EMZ individuals, respectively (Figure 4 and Table S2). In each population, few identical deletions (same borders) were present in several individuals, most likely resulting from the same initial rearrangements during the first meiosis. When considering the presence of identical deletions in several individuals resulting from a single event, a mean of 1.83 and 2.73 structural rearrangements were observed per RCC and EMZ individual, respectively (Figure 4).

The size of each deleted region was evaluated (Table S2), and ranged from 0.138 Mb (0.59% of A01) to 24.78 Mb (53.6% of C02: RCC-G3-077) in RCC-G3 plants, and from 0.048 Mb (0.12% of C01: EMZ-G3-092 and EMZ-G3-106) to 23.36 Mb (50.55% of C02: EMZ-G3-092) in EMZ G3 plants. Most of these deletions were in the distal region of a chromosome arm (79.41% in RCC and 67.4% in EMZ), and a few contained active centromeres (2.94% in RCC and 8.7% in EMZ) of the C01, C02, or C09 chromosomes. When cumulating the size of all deletions observed in the genome of a G3 plant, the individuals presenting the least and most important rearranged genomes corresponded to EMZ 072 (2.43 Mb: 0.37% of the Darmor genome) and EMZ 077 (35.00 Mb: 5.42% of the Darmor genome). Interestingly, no loss of any parental allele was identified for RCC-G3-075, despite the presence of 19 bivalents in only 57% of its PMCs (Table 1). To determine whether a subgenome was more prone to deletions, the total size of the deleted regions observed for each subgenome in a G3 population was determined. Subsequently, the proportion of deleted regions in a subgenome was calculated by taking into account the relative size of the A and C subgenome in B. napus. These analyses revealed significantly more deleted regions (total size) from the C than the A subgenome in both RCC (Wilcoxon ranksum test, W = 12,058, p-value = 3.49e205), and EMZ (Wilcoxon rank-sum test, W = 9123.5, p-value = 0.01823) (Figure 5). In addition, we observed that a larger number of C than A chromosomes were subjected to deletions. Indeed, the proportion of the deleted region per chromosome was not homogenous in RCC (Friedman rank-sum test, chi-squared = 133.5809, d.f. = 18, p-value , 2.2e216), or EMZ (Friedman rank-sum test, chi-squared = 52.1907, d.f. = 18, p-value = 3.497e205). Six and eight of the nine C chromosomes had a deleted fragment in the RCC and EMZ G3 populations, respectively, compared with three and five A chromosomes for RCC and EMZ. Deleted regions were never observed in the C07, A05, A06, A08, and A10 chromosomes in either of the genotyped RCC or EMZ-G3 individuals. Conversely, they were present in the C01, C02, C06, C09, A03, and A09 chromosomes of both populations, especially C02, C06, and C09. The C01 chromosome had the largest proportion of deleted regions in EMZ but not in RCC (Figure 5).

Identification of translocations (via HE) in resynthesized B. napus populations

The deletions (lost alleles) identified for each G3 plant may be the result of either a deletion or a translocation (via HE) that was fixed at the homozygous stage. To identify the type of rearrangement involved, we first examined the signal intensity of the SNP markers. Indeed, in the event of a translocation (after an HE), the signal intensity of the markers that showed a loss of one parental allele was hypothesized to be doubled for the other homeologous parental allele (Figure 3). To test this hypothesis, the relative copy number of such markers was estimated from the log 2 ratio of the signals observed in the individual that had lost a parental allele and the other diploid parent (data not shown). In addition, to evaluate the fluctuation in the hybridization signal for SNP markers presenting the same dosage, dominant markers (i.e., present in one diploid parent but absent in the other) were also used. No significant increase in copy number could be observed in any of the rearranged regions, suggesting that our data were not adapted to evaluate the allele dosage. Consequently, to unambiguously determine the type of rearrangement involved in these chromosome fragment losses, we performed BAC-FISH experiments using BAC probes from a deleted region of a G3 plant, and the Bob014O06 BAC clone [START_REF] Howell | Integration of the cytogenetic and genetic linkage maps of Brassica oleracea[END_REF] that specifically hybridizes to all C chromosomes in B. napus [START_REF] Szadkowski | The first meiosis of resynthesized Brassica napus, a genome blender[END_REF][START_REF] Suay | Crossover rate between homologous chromosomes and interference are regulated by the addition of specific unpaired chromosomes in Brassica[END_REF]. Due to the high sequence similarity between the A and C subgenomes, the BACs specific to a To identify putative structural rearrangements in each B. napus synthetic line (G3), only markers that were polymorphic between B. rapa and B. oleracea (i.e., "HDEM" and "Z1," or "RC34" and "C1.3"), and heterozygous in the first allotetraploid synthetic individual created ("EMZ-S0" or "RCC-S0") were considered. The position of these markers was then inferred using the B. napus cv. Darmor reference genome sequence assembly (version 4.1, Chalhoub et al. (2014)). This graphical representation shows an example of six markers that were polymorphic in the two diploid parental lines, which were present on C01 in B. oleracea (M1-M6), on A01 in B. rapa (M19-M69), and on the A01/C01 homeologous chromosomes (heterozygous markers) in B. napus S0. The loss of one parental allele from two consecutive markers (M1 and M2 on C01 herein) in a synthetic B. napus individual indicates the presence of a structural rearrangement resulting from either (A) a translocation (via homeologous exchange) or (B) a deletion. For these two cases, the reciprocal situation may also be observed deleted region also consistently hybridized to the homeologous chromosome in B. napus cv. Darmor, and thus provided four signals (Xiong and Pires 2011). Using a BAC clone present in a rearranged region of a G3 individual, two or four signals were expected in the case of a deletion or a translocation (via HE) that had been fixed at the homozygous stage, respectively. In all four cases tested (one EMZ and three RCC plants), four signals were identified in each G3 individual (Figure 6), indicating the presence of fixed translocations (via HE). These rearranged chromosomes harbored a dual color derived from both the C and A chromosomes. For example, the BAC specific to the C02 region deleted in RCC-G3-077 (Figure 6P) hybridized to four regions (red signals), and the Bob014O06 BAC clone specific to the C subgenome only hybridized to half of the C02 chromosome. These results indicated that the 25-Mb fragment lost on C02 was replaced by the homeologous A02 region. Similar observations were clearly observed between the A04 and C04 chromosomes in RCC-G3-055 (Figure 6H), for which a 6-Mb fragment from C04 was replaced with the homeologous A04 region. These BAC-FISH experiments not only facilitated the identification of the molecular mechanisms involved in the loss of certain genomic regions, but also permitted the observation of structural rearrangements that could not be identified using the 60K SNP array (indicated by stars in Figure 6). The rearrangements that were only identified by BAC-FISH may correspond to a translocation (via HE) present at the heterozygous stage.

The impact of structural rearrangements on gene content

We determined the effect of these large deletions on gene content using the B. napus Darmor reference genome composed of 101,040 genes [START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF]. A mean number of 1912 and 1689 genes were present in the deleted regions of the RCC and EMZ individuals, respectively (for details of the number of genes lost per deleted region, see Table S2), ranging from 0 to 4501 genes in RCC (4.45%: RCC-G3-053), and from 134 to 3595 genes in EMZ (3.56% EMZ-G3-088). As observed for the deletion size, there were more significant deletions in the C subgenome in both RCC (Wilcoxon rank-sum test, W = 12,052, p-value = 3.343e205) and EMZ (Wilcoxon rank-sum test, W = 9114.5, p-value = 0.01721). To assess the real impact of these deletions, the presence of a homeologous gene was determined using the homeologous gene relationship identified by [START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF]. In these homeologous regions, we observed that $75% of genes were present on both homeologs, whereas $15% of the deleted genes were not present on the homeologous chromosome (for details, see Table S2).

Putative gene conversion in resynthesized B. napus populations

In almost all G3 genotypes, we observed the presence of a single SNP marker that displayed a loss of one parental allele (mean of three SNPs per plant) on certain chromosomes. Such an allele loss could be derived from a deletion, a small translocation, or a gene conversion. To decipher the type of rearrangement involved, 10 primer pairs surrounding this type of SNP marker were designed. These primers (Table S3) were designed to consistently amplify homeologous regions of identical size (no indels). Direct sequencing of the PCR products obtained for the diploid parents revealed the presence of a polymorphism that could be detected (double peaks) in the resynthesized S0 B. napus. In contrast, the G3 allotetraploids presenting a loss of one parental allele did not exhibit any of these double peaks, validating the loss of one parental allele, but preventing the unequivocal determination of the origin of this deleted region (i.e., gene conversion, translocation, or deletion).

DISCUSSION

To mimic the early generations of a nascent polyploid, and to better understand the major challenges faced by B. napus after its formation, two resynthesized allotetraploid B. napus populations ("RCC" and "EMZ") were obtained for the first time by open pollination, enabling outcrossing. We evaluated their fertility, and compared it to that of the population obtained by manual selfing and to natural B. napus. Thereafter, in many G3 individuals obtained by open pollination, we evaluated the level of numerical aneuploids, identified SVs at the whole genome level, and determined their impact on the genome size, gene content, allelic diversity, and meiotic stability.

Aneuploidy and fertility of resynthesized B. napus populations

The level of aneuploid progenies was assessed in our resynthesized B. napus plants obtained via open pollination, and we observed a very low percentage of aneuploid progenies (,10%). This level of numerical aneuploids was much smaller than the $60% numerical aneuploid plants identified in resynthesized B. napus produced by controlled selfing (Xiong et al. 2011). In the latter study, there was no selection for plants with 38 chromosomes, similar to our populations obtained by open pollination. The level of aneuploidy observed in our populations obtained by open pollination was similar to those identified in resynthesized AT3 allotetraploid wheat (11.7%) presenting stable karyotypes (Zhang et al. 2013a), or in natural populations of Tragopogon neoallotetraploids (13.5% in T. mirus and 10% in T. miscellus) [START_REF] Chester | An assessment of karyotype restructuring in the neoallotetraploid Tragopogon miscellus (Asteraceae)[END_REF][START_REF] Chester | Patterns of chromosomal variation in natural populations of the neoallotetraploid Tragopogon mirus (Asteraceae)[END_REF]. Because of the negative impact of numerical aneuploidy on fertility (Xiong et al. 2011;Zhang et al. 2013b), the contribution of such aneuploids to the production of the following generation after outcrossing will decrease over time, and they will be less prevalent in the population.

The fertility of our two resynthesized B. napus populations (obtained after either open pollination or manual selfing) was assessed, and found to be very low compared to natural B. napus. It has been established that very frequent homeologous pairing occurs during the first meiosis of a resynthesized B. napus S0 plant [START_REF] Szadkowski | The first meiosis of resynthesized Brassica napus, a genome blender[END_REF], generating gametes that carry structural rearrangements. It is likely that some of these latter gametes are not viable, partially explaining the poor seed yield observed in resynthesized B. napus compared with natural oilseed rape. After controlled self-fertilization of resynthesized B. napus, structural rearrangements with more or fewer deleterious effects will be fixed, decreasing plant fertility. In contrast, open pollination under cages can generate progenies from either the selfing or outcrossing of plants carrying different rearrangements. This latter case will contribute to purging some of the deleterious modifications. Indeed, we observed that the decreased fertility was reduced in plants produced by open pollination compared to controlled selfing. Another factor affecting fertility is related to the origin of the diploid parental progenitors. Indeed, we observed a different level of fertility between our two populations (with EMZ being consistently less fertile than RCC plants). In our study, the B. oleracea and B. rapa individuals used to create the resynthesized B. napus belonged to different subspecies. In addition, the self-incompatibility alleles carried by the parental diploid species will also affect the fertility and the putative successful establishment of a novel polyploid. Indeed, controlled self-fertilization of our two populations was performed, and the EMZ plants were almost sterile, in accordance with the self-incompatibility of EMZ in contrast to the self-compatibility of RCC plants [START_REF] Hadj-Arab | Génétique de l'incompatibilité dans le genre Brassica L[END_REF].

To date, it is unknown whether the high fertility of natural B. napus resulted from the continuous selection of individuals with the highest karyotype stability and fertility, or to an unknown mechanism leading to the restoration of fertility soon after the formation of B. napus.

Structural rearrangements in resynthesized B. napus populations

Structural rearrangements were evaluated at the whole genome level in G3 resynthesized B. napus individuals using the Illumina 60K array. Using this array, a mean number of one SNP marker per 200-300 kb was polymorphic between the diploid parents used to create the resynthesized B. napus, enabling a relatively fine identification of SVs. We A, E, L, andM). The identification of a translocation between A01 and C01 in EMZ-G3-090 using Bna191G20 (in red) and Bob014O06 (in green) clones (A-D); between A04 and C04 in RCC-G3-055 using KBrH009I04 (in red), and Bob014O06 (in green) clones (E-H); between A07 and C06 in RCC-G3-072 using KBrB021P15 (in red) and Bob014O06 (in green) clones (I-L); and between A02 and C02 in RCC-G3-077 using KBrB086G22 (in red) and Bob014O06 (in green) clones (M-P). Additional translocations between the A and C chromosomes are indicated by a star ( Ã ). Chromosomes were counterstained with DAPI (blue). Bars represent 5 mm. A schematic representation of the expected structure of the translocated chromosomes (drawn to scale) is presented (D, H, L, and P).

observed that 32 of the 33 genetically studied G3 plants exhibited important SVs. We were able to determine that the identified fragment losses occurred significantly more on the C than the A-subgenome, which agreed with previous analyses of populations from resynthesized B. napus [START_REF] Song | Rapid genome change in synthetic polyploids of Brassica and its implications for polyploid evolution[END_REF][START_REF] Gaeta | Genomic changes in resynthesized Brassica napus and their effect on gene expression and phenotype[END_REF][START_REF] Szadkowski | The first meiosis of resynthesized Brassica napus, a genome blender[END_REF], or from B. napus varieties [START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF]. For the first time, by analyzing many synthetics derived from two different genetic backgrounds, we were able to determine that these losses may represent up to 5% of the whole genome.

To identify the molecular mechanisms involved in the losses of these genomic regions, BAC-FISH experiments were conducted using four different plants. In all cases, these losses were due to translocations (via HE). Even if all the identified fragment losses were the result of translocations, these structural rearrangements led to the reduced B. napus genome size because the lost C genomic regions were replaced by consistently smaller homeologous regions from the A genome. Among our plants obtained via open pollination, we were surprised to identify so many translocated regions that had been fixed at the homozygous stage after only three rounds of meiosis. This result indicates that the translocation of similar genomic regions was a highly frequent event. This observation is consistent with the high frequency of rearranged A01/C01 chromosomes in the gametes ($50%) of S0 resynthesized B. napus plants [START_REF] Szadkowski | The first meiosis of resynthesized Brassica napus, a genome blender[END_REF]. In our G3 populations, up to seven rearrangements and losses of up to 35 Mb were detected. However, it is important to note that all the losses observed in our material, and discussed thereafter, were most likely replaced by homeologous regions. Some chromosomes (A01, A09, C01, C02, C03, C04, C06, and CO9) presented genetic changes (losses) in individuals from both populations, while other chromosomes presented fragment losses in only one population (A02, A03, A07, A10, C05, and C08), or never presented any fragment losses in any synthetics. These results are in agreement with those of [START_REF] Gaeta | Genomic changes in resynthesized Brassica napus and their effect on gene expression and phenotype[END_REF], who observed that rearrangements were not randomly distributed across the genome in their aneuploid S5 plants from resynthesized B. napus. No fragment losses were observed in the A05, A06, A08, A10, or C07 chromosomes within our G3 populations. Interestingly, whole genome sequencing of five natural B. napus varieties also revealed an absence of fragment losses resulting from the homeologous exchange of all these chromosomes (except A05) [START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF]. The absence of structural variations on these chromosomes may be explained in several ways: (i) these chromosomes exhibit low synteny with their respective homeologous chromosomes (Xiong et al. 2011), (ii) there may be hot and cold regions of structural variation, and (iii) the loss of some chromosome regions may be deleterious. In addition to the nonrandom distribution of structural variations between chromosomes, most rearrangements occurred at the distal ends of chromosomes near the telomere or subtelomeric regions where most crossover events occur [START_REF] Nicolas | Non-random distribution of extensive chromosome rearrangements in Brassica napus depends on genome organization[END_REF]. Only a few rearrangements included the active centromeres of the C01, C02, and C09 chromosomes, which were consistently replaced by the corresponding homeologous centromeres. Overall, our results show that many regions of resynthesized B. napus genomes are highly rearranged, and are larger than those described for natural varieties of B. napus [START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF]. This increased genomic stability in natural B. napus may be the result of the continuous purging of the most detrimental rearrangements by human and natural selection [START_REF] Gaeta | Genomic changes in resynthesized Brassica napus and their effect on gene expression and phenotype[END_REF]). Another explanation relies on the presence of loci that suppress crossover formation between nonhomologous chromosomes at the time of B. napus formation, or soon after, enabling its successful establishment [START_REF] Jenczewski | Chromosomal and gene expression changes in Brassica allopolyploids[END_REF].

The SVs identified in our G3 plants may represent the tip of the iceberg, because only large translocations that were fixed at the homozygous stage, but not those fixed at the heterozygous stage, could be identified using the array (Figure 3). However, these latter rearrangements were observed in few G3 plants using GISH-like experiments (Figure 6). Due to the lower cost and common use of NGS technologies, a wider range of structural variants (indels, inversions, mobile element transpositions, or translocations) will be identified in the near future [START_REF] Saxena | Structural variations in plant genomes[END_REF][START_REF] Tattini | Detection of genomic structural variants from next-generation sequencing data[END_REF]. These technological advances will permit to finely determine how deeply the genome is rearranged, and will allow a better understanding of the molecular mechanisms involved in the evolutionary structural dynamics of nascent allopolyploid B. napus.

The impact of structural rearrangements on genome size, gene content, allelic diversity, and meiotic stability Among the G3 individuals studied herein, major losses of chromosome regions were observed. For the first time, we could determine the size and gene content of numerous resynthesized plants. We observed that a mean value of $4% of the whole genome size of B. napus was rearranged in a few generations. Overall, these rearrangements caused gene losses mainly from the C subgenome, with a loss of up to 10% of the genes belonging to the C subgenome. The impact of these deleted genes was most likely buffered by the replacement of most of these genes (77%) by their homeologs present on the A subgenome. Thus, if all the deletions identified in this study were due to translocations (via HE), then the overall gene content will not be dramatically decreased, but will strongly affect the genetic diversity in the rearranged regions. Duplication of a genomic region after HE may either have a detrimental, natural, or beneficial phenotypic effect, as exemplified by the increased seed yield [START_REF] Osborn | Detection and effects of a homeologous reciprocal transposition in Brassica napus[END_REF] or disease resistance [START_REF] Zhao | Quantitative trait loci for resistance to Sclerotinia sclerotiorum and its association with a homeologous non-reciprocal transposition in Brassica napus L[END_REF]. It has been shown that, for the main agronomical traits, some QTLs are often carried by homeologous regions with different genetic value in B. napus (i.e., quantitative blackleg resistance: Fopa Fomeju et al. 2014). Thus, the use of resynthesized B. napus individuals, which present a highly shuffled genome, will offer a new tactic to test the phenotypic impact of duplicating, or removing, a region carrying a QTL of agronomic interest.

For all genotyped G3 individuals, the meiotic stability was evaluated, which showed that it was strongly affected compared to natural B. napus. We compared the meiotic behavior determined in each G3 plant to the observed rearrangements (number, percentage, or position on the chromosome), but did not identify any clear correlation. Further studies are required to establish the relationship among the size, location, and number of translocations with meiotic stability of newly formed B. napus populations, and their impact on seed fertility.

In conclusion, our results, mimicking what probably occurred under natural conditions, highlight the extensive shuffling of the B. napus genome immediately after allopolyploidization. Overall, the identified structural variations that primarily result from translocations (via HE) led to a decrease in the size, gene content and allelic diversity of the genome. We show that open pollination of the resynthesized B. napus individuals within the initial generations enabled the natural selection of individuals with the most stable karyotype, and most likely purged some deleterious rearrangements. However, these resynthesized populations still presented unstable meiotic behavior and low fertility, indicating that further analyses of the structural and functional evolutionary dynamics of resynthesized B. napus populations are necessary to improve our current understanding of the key mechanisms involved in meiotic stabilization and the restoration of fertility.

Figure 1

 1 Figure 1 Schematic representing how the resynthesized B. napus populations used in this study were created. The names of the different B. oleracea and B. rapa varieties used to create these populations are indicated. The method (self-fertilization or open pollination) used to obtain the different generations is mentioned.

Figure 2

 2 Figure 2 Fertility of EMZ and RCC resynthesized B. napus populations obtained by either self-fertilization or open pollination, compared to the B. napus Drakkar variety. The fertility, corresponding to the number of seeds per 100 flowers (y-axis), was assessed for each population for three consecutive generations (Y1, Y2, and Y3). Each year, the fertility of the Drakkar variety was also determined by growing 50 Drakkar plants using the same growth conditions described for the resynthesized populations produced by open pollination.

Figure 3

 3 Figure 3 Methodology used to analyze the SNP genotyping data obtained from the Illumina Brassica 60K array and identification of structural rearrangements in each of the B. napus synthetic lines (third generation).To identify putative structural rearrangements in each B. napus synthetic line (G3), only markers that were polymorphic between B. rapa and B. oleracea (i.e., "HDEM" and "Z1," or "RC34" and "C1.3"), and heterozygous in the first allotetraploid synthetic individual created ("EMZ-S0" or "RCC-S0") were considered. The position of these markers was then inferred using the B. napus cv. Darmor reference genome sequence assembly (version 4.1,[START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF]). This graphical representation shows an example of six markers that were polymorphic in the two diploid parental lines, which were present on C01 in B. oleracea (M1-M6), on A01 in B. rapa (M19-M69), and on the A01/C01 homeologous chromosomes (heterozygous markers) in B. napus S0. The loss of one parental allele from two consecutive markers (M1 and M2 on C01 herein) in a synthetic B. napus individual indicates the presence of a structural rearrangement resulting from either (A) a translocation (via homeologous exchange) or (B) a deletion. For these two cases, the reciprocal situation may also be observed [i.e., duplication of a C01 region (A) or deletion of an A01 region (B)].

  Figure 3 Methodology used to analyze the SNP genotyping data obtained from the Illumina Brassica 60K array and identification of structural rearrangements in each of the B. napus synthetic lines (third generation).To identify putative structural rearrangements in each B. napus synthetic line (G3), only markers that were polymorphic between B. rapa and B. oleracea (i.e., "HDEM" and "Z1," or "RC34" and "C1.3"), and heterozygous in the first allotetraploid synthetic individual created ("EMZ-S0" or "RCC-S0") were considered. The position of these markers was then inferred using the B. napus cv. Darmor reference genome sequence assembly (version 4.1,[START_REF] Chalhoub | Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome[END_REF]). This graphical representation shows an example of six markers that were polymorphic in the two diploid parental lines, which were present on C01 in B. oleracea (M1-M6), on A01 in B. rapa (M19-M69), and on the A01/C01 homeologous chromosomes (heterozygous markers) in B. napus S0. The loss of one parental allele from two consecutive markers (M1 and M2 on C01 herein) in a synthetic B. napus individual indicates the presence of a structural rearrangement resulting from either (A) a translocation (via homeologous exchange) or (B) a deletion. For these two cases, the reciprocal situation may also be observed [i.e., duplication of a C01 region (A) or deletion of an A01 region (B)].

Figure 4

 4 Figure 4 Graphical representation of the deletions identified in each "RCC" (left) or "EMZ" (right) G3 plant. The B. napus chromosomes belonging to the A and C subgenomes are shown in blue and red, respectively. The size of each chromosome (in megabases) is indicated above each chromosome, and a ruler is drawn underneath each chromosome, with larger and smaller tick marks every 10 and 2 Mb, respectively. The locations of active centromeres (Mason et al. 2016) are indicated by green rectangles under each chromosome. Each circle represents one G3 plant. For each individual, a black rectangle represents a deleted region. The last inner circle corresponds to the homeologous relationship between the A and C chromosomes.

Figure 5 Figure 6

 56 Figure5Mean size of subgenome deletions (in megabases) in EMZ (A) or RCC G3 plants (C). The mean size of deletions per chromosome (in megabases) for EMZ (B) or RCC plants (D) is also presented. The presence of significant differences using the Wilcoxon rank-sum test[START_REF] Wilcoxon | Individual comparisons by ranking methods[END_REF]) is indicated by a star(s).

Volume 7 February 2017 | Structural Dynamic of Synthetic B. napus | 707

Volume 7 February 2017 | Structural Dynamic of Synthetic B. napus | 715

ACKNOWLEDGMENTS

We thank Reece Tollenaere for critically reading the manuscript, and the two anonymous reviewers for their careful evaluation of our manuscript and the many insightful comments and suggestions. We thank the Biological Resource Center (CRB BrACySol, Rennes-Ploudaniel, France) for providing seeds of B. oleracea, B. rapa, and natural B. napus genotypes. We would also like to thank all the technical staff of the experimental unit "la Motte" (INRA, UE30787, Le Rheu, France), and at the greenhouse, for management of the plant material (especially Laurent Charlon, P. Rolland, Jean-Paul Constantin, Jean-Marie Lucas, Frederic Letertre and C. Guerin). We thank the UMR 8199 genotyping service (Lille, France), which belongs to the "Federation de Recherche" 3508 Labex EGID (European Genomics Institute for Diabetes; ANR-10-LABX-46), and was supported by the ANR Equipex 2010 session (ANR-10-EQPX-07-01; "LIGAN-PM": http://www.good.cnrs.fr/). This work was made possible by financial support from ERA-CAPS13.017_Evo-Genapus.

LITERATURE CITED


Paolo Ballarini
email: paolo.ballarini@ecp.fr

Emmanuelle Gallet
email: emmanuelle.gallet@ecp.fr

Pascale Le Gall
email: pascale.legall@ecp.fr

Matthieu Manceny
email: matthieu.manceny@isep.fr

Formal analysis of the Wnt/β-catenin pathway through statistical model checking

Keywords: HASL Model Checking, Stochastic modelling, biological pathways, Wnt/β-catenin

The Wnt/β-catenin signalling pathway plays an important role in the proliferation of neural cells, and hence it is the main focus of several research aimed at understanding neurodegenerative pathologies. In this paper we consider a compact model of the basic mechanisms of the Wnt/β-catenin pathway and we analyse its dynamics by application of an expressive temporal logic formalism, namely the Hybrid Automata Stochastic Logic. This allows us to formally characterise, and effectively assess, sophisticated aspects of the Wnt/β-catenin pathway dynamics.

Introduction

Systems Biology [START_REF] Kitano | Foundations of Systems Biology[END_REF] is concerned with the development of formalisms for building "realistic" models of biological systems, i.e. models capable of reproducing wet-lab observations. A biological model consists of a set biochemical agents (i.e. species) whose interactions are expressed by a set of reaction equations. This leads to either a continuous-deterministic interpretation (i.e. in terms of a system of differential equations), or to a discrete-stochastic interpretation (i.e. in terms of a discrete-state stochastic process).

Stochastic modelling and systems biology. Within the discrete-stochastic semantics realm, which is what we consider in this work, molecular interactions are assumed to be of stochastic nature hence biochemical reactions occur according to probability distributions. In this case what modellers normally do is to generate one (or several) trajectory(ies) through stochastic simulation and observe the evolution of the species (under different model's configurations) in order to figure out how a given aspect of the model's dynamics is affected by the various elements of the model (i.e. what species/reactions is responsible for a given observed behaviour). Such an approach has two main advantages: its simplicity and its low computational cost (the runtime for generating a single trajectory or a normally small number of trajectories is very low even for large models). On the other hand the main disadvantage is that it is little formal, meaning that the modeller must draw conclusions based only on the observation of a single (stochastic) trajectory (or of a trajectory obtained by averaging a normally small number of trajectories).

Stochastic model checking and systems biology. Stochastic model checking [START_REF] Kwiatkowska | Stochastic model checking[END_REF] (SMC) is a formal technique that allows the modeller to formally express relevant properties in terms of a (stochastic) temporal logic and to assess them against a given stochastic model. This is achieved through an automatic procedure which can either provide an exact answer through exhaustive exploration of the model's state space (i.e. numerical model checking [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF]) or an estimated answer resulting from a finite sampling of the model's trajectory (i.e. statistical model checking [START_REF] Legay | Statistical model checking : An overview[END_REF]). SMC has at least two main advantages with respect to informal approaches: first it provides the modeller with a language for capturing relevant properties formally; second the answer it calculates (e.g. probability that a property is satisfied by the model) are either exact (i.e. they reflect the complete set of possible behaviours of the model) or are accurate estimates (i.e. calculated over a a sufficiently large sample of trajectories). The effectiveness of SMC in systems biology applications is demonstrated by an ever increasing number of publications, e.g. [START_REF] David | Runtime verification of biological systems[END_REF][START_REF] Heath | Probabilistic model checking of complex biological pathways[END_REF][START_REF] Ballarini | Expressive statistical model checking of genetic networks with delayed stochastic dynamics[END_REF].

β-catenin and the WNT pathway. In cellular biology signalling pathways are basic mechanisms responsible for controlling a cell's life-cycle. Simply speaking a signalling pathway represents a cascade of biochemical reactions which is triggered by a specific signal (i.e. type of molecules) whose presence, normally at the cell membrane, activates the cascade leading to the "transmission" of the signal inside the cell (i.e. cytosol and/or nucleus). In this paper we study a model of the Wnt/β-catenin pathway, a signalling pathway known to be involved in the pathological degeneration of neuronal cells [START_REF] Mazemondet | Quantitative and kinetic profile of wnt/β-catenin signaling components during human neural progenitor cell differentiation[END_REF].

Our contribution. In this work we present preliminary results of application of formal analysis, based on the so-called Hybrid Automata Stochastic Logic (HASL) statistical model checking, to a model of the Wnt/β-catenin pathway presented in [START_REF] Mazemondet | Elucidating the sources of β-catenin dynamics in human neural progenitor cells[END_REF]. In particular we show how one can define specific HASL formulae for assessing sophisticated characteristics of the Wnt/β-catenin pathway dynamics. This includes, for example, measuring the temporal location and the amplitude of transient peaks of nuclear β-catenin, exhibited by certain initial conditions, or assessing its oscillatory character resulting from other conditions. If in [START_REF] Mazemondet | Elucidating the sources of β-catenin dynamics in human neural progenitor cells[END_REF] the analysis of the Wnt/β-catenin model is simply done through plotting of simulated trajectories, here we move analysis to a higher and more formal level by demonstrating how, through model checking, one gains access to the analysis of sophisticated dynamical aspects of the Wnt/β-catenin pathway.

Paper organisation. We introduce the Wnt/β-catenin mechanism in Section 2 and describe the model presented in [START_REF] Mazemondet | Elucidating the sources of β-catenin dynamics in human neural progenitor cells[END_REF] which we have used for our analysis. In Section 3 we give a concise description of the HASL statistical model checking formalism. In Section 4 we present the results obtained by application of HASL model checking to the analysis of Wnt/β-catenin model. We wrap up the paper with some conclusive remarks and future perspectives in Section 5.

A model of the Wnt/β-catenin pathway

Neurodegeneration is the process of progressive lost of structure/function of neuronal cells (i.e. neurons) which is at the basis of many neurodegenerative diseases, such as, for example, the Parkison's disease, Alzheimer's disease and the Amyotrophic lateral sclerosis. Research in this field is particularly focused on the growth of in vitro population of neural cells that may potentially be used in replacement therapies for neurodegenerative diseases. Cultivated cells undergo so-called proliferation, a process of successive cell divisions and potential differentiation into neurones and glial cells.

The Wnt/β-catenin pathway is a signalling pathway known to be involved in the proliferation/differentiation of neural cells. Specific in vitro experiments [START_REF] Mazemondet | Quantitative and kinetic profile of wnt/β-catenin signaling components during human neural progenitor cell differentiation[END_REF] have exhibited a high activity of the Wnt/β-catenin pathway during the differentiation of ReNcell VM (RVM) cells, i.e. a type of cells derived from the brain of a fetus and that are believed to be an appropriate model for replacement therapies in neurodegenerative pathologies. The activity of the Wnt/β-catenin is summarised as follows: in absence of extracellular Wnt molecules (normally at cell's membrane), a degradation complex causes the phosphorylation and subsequent destruction of β-catenin located in the cell's cytosol (denoted β cyt ); on the other hand in presence of Wnt proteins, the degradation complex is inactivated resulting in accumulation of β cyt . Furthermore from the cytosol βcatenin undergoes a (reversible) relocation to the nucleus (denoted β nuc ) wherein it activates the expression of one component of its degradation complex, i.e. the Axin protein. The above described mechanism is captured by a core-version of the Wnt/β-catenin pathway model presented in [START_REF] Mazemondet | Elucidating the sources of β-catenin dynamics in human neural progenitor cells[END_REF]. This consists of the twelve biochemical reactions illustrated by equations [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF].

R1 :W nt k 1 -→ ∅ R7 : AxinP + βcyt k 7 -→ AxinP R2 :W nt + AxinP k 2 -→ W nt + Axin R8 : ∅ k 8 -→ βcyt R3 :AxinP k 3 -→ Axin R9 : βcyt k 9 -→ ∅ R4 :Axin k 4 -→ AxinP R10 : βcyt k 10 -→ βnuc R5 :AxinP k 5 -→ ∅ R11 : βnuc k 11 -→ βcyt R6 :Axin k 6 -→ ∅ R12 : βnuc k 12 -→ Axin + βnuc (1)
The model consists of three basic molecular species: the Axin protein, which can be either in normal (Axin) or phosphorylated (AxinP ) form, the Wnt protein (W nt) and the β-catenin which can be either located in the cytosol (β cyt ) or in the nucleus (β nuc ). Equations (1) account for the following aspects: two reversible events, i.e. the phosphorylation of Axin (reactions R 4 and R 3 ) and the relocation of β-catenin from/to cytosol/nucleus (reactions R 10 and R 11 ); the W nt enhanced de-phosphorylation of Axin (reaction R 2 ); the nuclear β-catenin (i.e. β nuc ) regulated expression of Axin (reaction R 12 ); the phosphorylated AxinP enhanced degradation of cytosolic β-catenin (i.e. β cyt ) (reaction R 7 )3 ; the constant (DNA regulated) expression of cytosolic β-catenin (i.e. β cyt ) (reactions R 8 ); the degradation of all species i.e. W nt (reactions R 1 ), Axin in either form (reactions R 5 or R 6 ) and β cyt (reaction R 9 ).

In this paper we focus on the discrete-stochastic interpretation of Equations (1), hence species populations are expressed in terms of number of molecules and reactions are of stochastic nature and are assumed to obey the mass action law (meaning that a reaction's rate is proportional to the current population of the reactants, except for R 8 whose rate is constant). With respect to the model configuration we consider two basic sets A and B of parameter values respectively taken from the sets 3 and 4 in [START_REF] Mazemondet | Elucidating the sources of β-catenin dynamics in human neural progenitor cells[END_REF], and indicated in Table 1 Throughout the remainder of the paper we will analyse the Wnt/β-catenin model by comparing the dynamics corresponding to the two parameter sets of Table 1. Following [START_REF] Mazemondet | Elucidating the sources of β-catenin dynamics in human neural progenitor cells[END_REF], we will also consider two variants of the basic model (1). The first variant, denoted Wnt-inject, represents a single injection of an extra amount (i.e. 1000) of W nt molecules in the system at a fixed delay d i . The second variant, denoted Wnt-doped, represents the presence of a doping mechanism that kicks in at a given delay d d and then it sustainably generates a fresh Wnt molecule at given frequency (assumed to be exponential distributed with parameter k d ).

Stochastic Petri Net model of the Wnt/β-catenin pathway

The COSMOS [START_REF] Ballarini | COSMOS: a statistical model checker for the hybrid automata stochastic logic[END_REF] model checker which we used for analysing the Wnt/β-catenin pathway model uses Generalised Stochastic Petri Net (GSPN) [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF] as modelling formalism. A GSPN model is a bipartite graph consisting of two classes of nodes, places and transitions. Places (circle nodes) may contain tokens (representing the state of the modelled system) while transitions (bar nodes) indicate how tokens "flow" within the net (encoding the model dynamics). The state of a GSPN consists of a marking indicating the distribution of tokens throughout the places (i.e. how many tokens each place contains). A transition is enabled whenever all of its input places contains a number of tokens greater than or equal to the multiplicity of the corresponding (input) arc. An enabled transition may fire consuming tokens (in a number indicated by the multiplicity 4 of the corresponding input arcs) from all of its input places and producing tokens (in a number indicated by the multiplicity of the corresponding output arcs) in all of its output places. Transitions can be either timed (denoted by thick bars) or immediate (denoted by thin bars). Timed transitions are associated to a probability distribution (e.g. Exponential, Uniform, Deterministic, etc).

In the context of this paper GSPN places represent biological species (and their marking the molecular population of a species), whereas timed transitions represent chemical reactions. For more details on GSPN we refer the reader to the literature [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF]. Figure 1 depicts the GSPN model encoding the Wnt/β-catenin chemical equations [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF]. The net contains a place for each species of the Wnt/β-catenin model and a transition for each reaction. Non filled-in transitions are exponentially distributed (with marking dependent rate) with rate-constant correspond-ing to that of either parameter set of Table 1. The sub-net enclosed in dashed line box (top left) has been added in order to add the Wnt-inject and Wnt-dope behaviour to the basic model. In order to study the behaviour of the Wnt-inject (Wnt-dope) variant it suffices to add one token in the initial marking of place W nt_inject (W nt_dope). Notice that the black filled-in transitions of the Wntinject/Wnt-dope subnet are associated to deterministic delays.

Observing β nuc dynamics on a single stochastic trajectory. Following [START_REF] Mazemondet | Elucidating the sources of β-catenin dynamics in human neural progenitor cells[END_REF] we first look at the dynamics of nuclear β-catenin (i.e. β nuc ) as observed along a single trajectory simulated over 24 hours (Figure 2 and Figure 3 with the units on x axis being minutes). Figure 2 compares the behavior of β nuc for the two parameter sets (Table 1) of the basic Wnt/β-catenin model in presence of Wnt (initial population of Wnt set to 1000). The interpretation in this case is quite straightforward: the presence of an initial Wnt signal (i.e. of 1000 molecules) triggers a (delayed) peak in β nuc which however quickly ends due to the steady degradation (and absence of reintegration) of Wnt. Eventually when all Wnt has faded away β nuc noisily converges to a certain level. Figure 3

HASL statistical model checking

The Hybrid Automata Stochastic Logic (HASL), introduced in [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF], extends Deterministic Timed Automata (DTA) logics for addressing Markov chain models [START_REF] Donatelli | Model checking timed and stochastic properties with CSL T A[END_REF][START_REF] Chen | Quantitative model checking of CTMC against timed automata specifications[END_REF], by employing Linear Hybrid Automata (LHA) as instruments for addressing a general class of stochastic processes, namely that of Discrete Event Stochastic Processes (DESP). An HASL formula φ ≡ (A, Z) consists of two elements: 1) A, a synchronising LHA (i.e. an LHA enriched with DESP state and/or event indicators) and 2) Z a target expression (see grammar [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF]) which expresses the quantity to be evaluated (either a measure of probability or, more generically, any real-valued measure).

Thus given a DESP model D and a formula φ ≡ (A, Z) the HASL model checking procedure employs stochastic simulation to samples trajectories of the synchronised process D × A, and then use (only) the paths selected by A (i.e. those paths of D ×A that reach an accepting location of A) for estimating the confidence-interval of the target measure Z. Such a procedure is implemented within the COSMOS [START_REF] Ballarini | COSMOS: a statistical model checker for the hybrid automata stochastic logic[END_REF] model checking framework, a tool which belongs to the fast expanding family of statistical model checkers (e.g. [START_REF] Boyer | Plasma-lab: A flexible, distributable statistical model checking library[END_REF][START_REF] Sen | VESTA: A statistical model-checker and analyzer for probabilistic systems[END_REF][START_REF] Håkan | Ymer: A statistical model checker[END_REF]). For practical reasons in HASL (and in particular within COSMOS) we employ GSPN as highlevel language for expressing a DESP. Below we give the definition of DESP and LHA and informally describe the synchronisation process (we refer the reader to [START_REF] Ballarini | HASL: an expressive language for statistical verification of stochastic models[END_REF] for a more formal treatment).

Definition 1 (DESP).

A DESP is a tuple D = S, π 0 , E, Ind, enabled, delay, choice, target where -S is a (possibly infinite) set of states, -π 0 ∈ dist(S) is the initial distribution on states, -E is a set of events, -Ind is a set of functions from S to R called state indicators (including the constant functions),

enabled : S → 2 E are the enabled events in each state with for all s ∈ S, enabled(s) = ∅. -delay : S × E → dist(R + ) is a partial function defined for pairs (s, e) such that s ∈ S and e ∈ enabled(s). where dist(A) denotes the set of distributions whose support is A.

-choice : S × 2 E × R + → dist(E)
Dynamics of a DESP. A configuration of a DESP consists of a triple (s , τ , sched) with s being the current state, τ ∈ R + the current time and sched : E → R + ∪ {+∞} being the function that describes the occurrence time of each scheduled event (+∞ if an event is not yet scheduled). The evolution (i.e simulation) of a DESP D can be informally summarised in terms of an iterative procedure consisting of the following steps (assuming (s , τ , sched) is the current configuration of D): 1) determine the set E m of events enabled in state s and with minimal delay δ m ; 2) select the next event to occur e next ∈ E m by resolving conflicts (if any) between concurrent events through probabilistic choice according to choice(s, E m , τ ); 3) determine the new configuration of the process resulting from the occurrence of e next , this in turns consists of three sub-steps: 3a) determine the new state resulting from occurrence of e next , i.e. s = target(s, e next , δ m ); 3b) update the current time to account for the delay of occurrence of e next , i.e. τ = τ +δ m ; 3c) update the schedule of events according to the newly entered state s (this implies setting the schedule of no longer enabled events to +∞ as well as determining the schedule of newly enabled events by sampling through the corresponding distribution). The above procedure maps directly on GSPN models, in which case the set of states S corresponds to the set of possible markings of a GSPN, the events E correspond to the (timed) transitions of a GSPN, and the remaining elements (i.e. delay, choice and target) are determined by the semantics of GSPN (i.e. the so-called token game).

Definition 2. A synchronised Linear Hybrid Automaton is a tuple A = E, L, Λ, I, F, X, flow, → where:

-E is a finite alphabet of events; -L is a finite set of locations; -Λ : L → P rop is a location labelling function; -I ⊆ L is the initial locations; -F ⊆ L is the final locations; -X = (x 1 , ...x n ) is a n-tuple of data variables; flow : L → Ind n associates an n-tuple of indicators with each location (projection flow i denotes the flow of change of variable x i ).

-→⊆ L × (Const × 2 E ) (lConst × { }) × Up × L is the set of edges of the LHA
where denotes the disjoint union, Const and lConst denotes the set of possible constraints, respectively left closed constraints, associated with A (see description below) and Up is the set of possible updates for the variables of A. Furthermore A fulfils the following conditions.

-c1 (initial determinism): ∀l = l ∈ I, Λ(l)∧Λ(l ) ⇔ false. This must hold whatever the interpretation of the indicators occurring in Λ(l) and Λ(l ). -c2 (determinism on events):

∀E 1 , E 2 ⊆ E : E 1 ∩ E 2 = ∅, ∀l, l , l ∈ L, if l γ,E1,U ----→ l and l γ ,E2,U
-----→ l are two distinct transitions, then either Λ(l) ∧ Λ(l ) ⇔ false or γ ∧ γ ⇔ false. Again this equivalence must hold whatever the interpretation of the indicators occurring in Λ(l), Λ(l ), γ and γ .

-c3 (determinism on ): ∀l, l , l ∈ L, if l γ, ,U ---→ l and l γ , ,U ----→ l are two distinct transitions, then either Λ(l) ∧ Λ(l ) ⇔ false or γ ∧ γ ⇔ false.

-c4 (no -labelled loops): For all sequences ----------→ l n such that l 0 = l n , there exists i ≤ n such that E i = .

l 0 γ0,E0,U0 ------→ l 1 γ1,E1,U1 ------→ • • • γn-1,En-1,Un-1 -
Synchronisation of LHA and DESP. The role of a synchronised LHA A is to select specific trajectories of a corresponding DESP D while collecting relevant data (maintained in the LHA variables) along the execution. For the sake of brevity we omit the formal semantics of the product process D × A in this paper, but we provide an intuitive description of it.

A state of the D ×A process is described as a triple (s, l, ν) where s is the current state of the DESP, l the current location of the LHA and ν : X → R the current valuation of the LHA variables. The synchronisation starts from the initial state (s, l, ν), where s is an initial state of the DESP (i.e. π 0 (s) > 0), l is an initial location of the LHA (i.e. l ∈ I) and the LHA variables are all initial set to zero (i.e. ν = 0) 5 . Notice that, by initial determinism, for every s ∈ S there is at most one l ∈ I such that s satisfies Λ(l). From the initial state the synchronisation process evolves through transitions where each transition corresponds to traversal of either a synchronised or an autonomous edge of the LHA (notice that because of the determinism constraints of the LHA edges at most only one autonomous or synchronised edge can ever be enabled in any location of the LHA. Furthermore if an autonomous and a synchronised edge are concurrently enabled the autonomous transition is taken first). If in the current location of the LHA (i.e. location l of the current state (s, l, ν) of process D×A) there exists an enabled autonomous edge l γ, ,U ---→ l , then that edge will be traversed leading to a new state (s, l , ν ) where the DESP state (s) is unchanged whereas the new location l and the new variables' valuation ν might differ from l, respectively ν, as a consequence of the edge traversal. On the other hand if an event e (corresponding to transition s e -→ s ) triggered by process D occurs in state (s, l, ν), either an enabled synchronous edge l γ,E ,U ----→ l (with e ∈ E ) exists leading to new state (s , l , ν ) of process D × A (from which the synchronised process will proceed) or the system goes to a dedicated rejecting state ⊥ and the synchronisation halts (indicating rejection of the trace).

HASL expressions. The second component of an HASL formula is an expression related to the automaton. Such an expression, denoted Z, is based on moments of a path random variable Y and defined by the grammar [START_REF] Baier | Model-checking algorithms for CTMCs[END_REF].

Z ::= c | E[Y ] | Z + Z | Z -Z | Z × Z | Z/Z Y ::= c | Y + Y | Y × Y | Y /Y | last(y) | min(y) | max(y) | int(y)| avg(y) y ::= c | x | y + y | y × y | y/y (2)
Z represents the expectation of an arithmetic expression based on LHA data variables and which uses path operators such as: last(y) (i.e. the last value of y along a synchronising path), min(y) (max(y)) the minimum (maximum), value of y along a synchronising path), int(y) (i.e. the integral over time along a path) and avg(y) (the average value of y along a path). In recent updates the COSMOS model checker [START_REF] Ballarini | COSMOS: a statistical model checker for the hybrid automata stochastic logic[END_REF] has been enriched with operators for assessing the Probability (Cumulative) Distribution Function (PDF/CDF) of the value that an expression Y takes at the end of a synchronising path. This requires specifying a discretised support of Y through the following syntax: Z = P DF (Y, s, l, h) which means that the probability of Y to take value in any sub-interval of fixed width s corresponding to the partition of the considered [l, h] support of Y is going to be evaluated (assuming that [l, h] is discretised in hl/s sub-intervals). Example Figure 4 shows a simple example of LHA that synchronises with the Wnt/β-catenin GSPN model of Figure 1. It uses three data variables: a clock t (storing the simulation time), a counter n (counting the number of occurrences of reactions R 2 ) and a variable b which keeps track of the population of β nuc . In the initial location l 0 the clock variable t grows with constant flow ṫ = 1, whereas n and b flows is null. On occurrence of R 2 the top synchronising self-loop edge on l 0 is traversed hence n is incremented whereas on occurrence of any other reaction the bottom self-loop on l 0 is traversed, hence n is not updated. On the hand b is updated with the current value of β nuc on occurrence of any reaction. As soon as N occurrences of R 2 have been observed the autonomous edge l 0 → l 1 is traversed and synchronization halts (reaching of accepting location l 1 ). Below few examples of complete HASL formulae composed with the LHA of Figure 4.

Model analysis through HASL formulae

In order to analyse the dynamics of the Wnt/β-catenin model we define a number of HASL formulae dedicated to capturing specific dynamical aspects of the GSPN model in Figure 1.

4.1 Measuring the maximal peaks of β nuc resulting from an unsustained Wnt signal.

Both the Wnt-basic and Wnt-inject models are designed to study the behaviour of the Wnt/β-catenin pathway in presence of an unsustained Wnt signal: i.e. a given amount of initial Wnt signal is present in the system but is steadily being consumed (reaction R 1 ) without being reintegrated (Wnt-basic) or being reintegrated once after a delay d i (Wnt-inject). The effect of a non-reintegrated Wnt signal results is the production of a single peak of β nuc (Figure 2) whereas a single, delayed, reintegration of 1000 Wnt molecules produces a second, shifted peak (Figure 3 left) in the population of β nuc . We introduce some HASL formulae for formally measuring the time location and the amplitude of such β nuc transient peaks. Observe that the analysis of a stochastic model through observation of a single simulated trajectory (as proposed in [START_REF] Mazemondet | Elucidating the sources of β-catenin dynamics in human neural progenitor cells[END_REF]) is in general little informative l and even more so in this case as repeated trajectories of the wntb model exhibit a rather large variance. In the light a more formal approach is vital to obtain a meaningful analysis.

Automaton A peaks . The LHA in Figure 5 is conceived for locating the maximal and minimal peaks along an alternating trace of a given observed species, in this case β nuc . The automaton uses a number of data variables (Table 2) and is dependent on two configuration aspects: the setting of a parameter δ (the chosen noise level, see below) and the partition of the event set E = E +βnuc ∪E -βnuc ∪ E =βnuc where E +βnut , E -βnuc and E =βnut are the events yielding respectively: an increase of β nuc , a decrease of β nuc and having no effect on β nuc population. Specifically, for model [START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF] we have:

E +βnuc = {R 10 }, E -βnuc = {R 11 } and E =βnuc = {R 1 , R 2 , R 3 , R 4 , R 5 , R 6 , R 7 , R 8 , R 9 , R 12 }.
A peaks consists of an initial location start, a final location end, and 4 intermediate locations (Min,Inc, Max and Dec) where the actual analysis of the synchronised trajectory takes place. From start (on entering of which the initial amount of β nuc is stored in x) the processing of the simulated trajectory leads either to Min or Max depending if the an increase (decrease) of β nuc above (below) the chosen level of noise δ is observed (i.e. A peaks copes fine both in the case that the observed species initially increases or decreases). Once in location Min (Max) the behaviour of the automaton depends on the type of observed event. If an event e ∈ E +βnuc (e ∈ E -βnuc ) is observed then location noisyInc (noisyDec) is entered indicating that β nuc has increased (decreased) although the increase (decrease) has not (yet) exceeded δ (with respect to the most recent detected minimum (maximum) previously stored in x). On the other hand if while in Min (Max) an event e ∈ E -A (e E +βnuc ) is observed, then this means that the current value of β nuc went below (above) the previously detected minimum (maximum) hence x must be updated with the newly found (potential) minimum (maximum) x := β nuc . Finally an occurrence of any event e ∈ E =βnuc while in Min or Max is simply ignored. From noisyInc (noisyDec) the processing of input trace may lead back to Min (Max) if β nuc re-decreases (re-increases) below (above) x (hence requiring an update x := β nuc ) or it may lead to Max ( Min)as soon as β nuc has increased (decreased) above (below) the noise level (i.e. x > β nuc -δ). The autonomous edge noisyInc→Max (noisyDec→Min) is traversed as soon as the value stored in x corresponds to an actual minimum (maximum) along the processed β nuc trace (i.e. this is the case when the current value of β nuc gets δ molecules far away from that stored in x). Hence when traversing noisyInc→Max (noisyDec→Min) we are sure that x contains a minimum (maximum) thus its value and its occurrence time are stored in the n th element of the array xmin[n] := x (xmax[n] := x), respectively tmin[n] := t (tmax[n] := t). The processing terminates (entering of End from any other location) as soon as the simulation time is t = T at which point all detected maxima and minima are stored in A peaks data variables. of β nuc peaks). We performed a number of experiments aimed at addressing this aspect, specifically we assessed φ xmax and φ tmax against different instances of the Wnt/β-catenin model with delayed Wnt re-injection (i.e. the Wnt-inject model) where each instance corresponds to a different value of k 1 (the decay rate Wnt). Figure 6 displays the results concerning the evaluation of the first and second peak of β nuc . They indicate that both the average height (left) and the average occurrence time (right) of the first and second peaks of β nuc decrease as the Wnt decay rate k 1 increases6 . Figure 6 left also shows that the second peak of β nuc (induced by Wnt re-injection) has, on average, a smaller amplitude than the first one and with a roughly constant difference of about 10% less between the two except for a k 1 = 0.1 for which the first and second peak's amplitude differs of about 5%.

Measuring the PDF of occurrence time of β nuc peaks. Figure 7 displays the PDF of the first (left) and second (right) peak of β nuc obtained by evaluation of φ P DF max against the Wnt-inject model (parameter set B). Both PDF curves exhibit a slight long-tail character with the majority of points being to the right of the maximum likely occurrence time.

Conclusion

We have presented a formal study of a stochastic model of the Wnt/β-catenin pathway, a biological mechanism with a relevant role in controlling the life-cycle of neuronal embryonal cells. This model has been previously considered [START_REF] Mazemondet | Elucidating the sources of β-catenin dynamics in human neural progenitor cells[END_REF] however it was analysed only informally, i.e. through observation of simulated trajectories. By means of a powerful formalism (i.e. HASL model checking) we formally characterised and accurately assessed a number of relevant aspects of the Wnt/β-catenin dynamics. In particular in this work we have focused on studying the effects induces on nuclear β-catenin (a basic element of the Wnt/β-catenin pathway) by the presence of a degrading (possibly reintegrated) Wnt signal. That included measuring of the average value and the PDF of the occurrence time and the amplitude of the β nuc peaks resulting from a transitory Wnt signal. We plan to evolve this preliminary study in several directions, including the formal analysis of the effects induced by a sustained Wnt signal, as well as the analysis of the dynamics of β nuc over a population of asynchronously evolving cells.

Fig. 1 .

 1 Fig. 1. GSPN model corresponding to Equations (1) of the Wnt/β-catenin pathway.

Fig. 2 .

 2 Fig. 2. Dynamics of βnuc along a 24 hours single trajectory of the Wnt/β-catenin pathway model with parameter set A (left) and set B (right)

Fig. 3 .

 3 Fig. 3. Dynamics of βnuc along a single trajectory of the Wnt/β-catenin pathway model with re-injection of Wnt at t = 450 minutes (left) and with doping of Wnt started at t = 150 minutes (right) and parameter set B

  is a partial function defined for tuples (s, E , d) such that E ⊆ enabled(s) and such that the possible outcomes of the corresponding distribution are restricted to e ∈ E .target : S×E× R + → S is a partial function describing state changes through events defined for tuples (s, e, d) such that e ∈ enabled(s).

Fig. 4 .

 4 Fig. 4. Simple example of LHA that synchronises with the Wnt/β-catenin GSPN model of Figure 1: the automaton selects paths containing N occurrences of reaction R2.

-φ 1 ≡

 1 (A, E[last(t)]): representing the average time for observing N occurrences of R 2 . -φ 2 ≡ (A, E[max(b)]): representing the maximum population reached by β nuc within the first N occurrences of R 2 . -φ 3 ≡ (A, P DF (last(t), 0.1, 0, 10)): representing the PDF of the delay for observing N occurrences of R 2 (computed over the interval [0, 10] with a discretisation step of 0.1)

1 ]: 1 E

 11 , ( x > n ) , { x + n ,>,{x := n } E = n ,>,; n ,>,{x := n } E n , > , ; ],(x> n + )^(t<T), {x := n , xmax[n]:=x, tmax[n]:=t, ¬Up? n++} ] t = T ),; ],(t = T ),; ],(t = T ),; ],(t = T ),; ],(x> n + ),{x := n , Up := ?} {x := n } ],(x> n ),{x := n , Up := >} ],(x< n )^(t<T), {x := n , xmin[n]:=x, tmin[n]:=t, ¬Up? n++}

Fig. 5 .

 5 Fig. 5. A β_peaks : an LHA for locating the maximal peaks (up to noise level δ) of βn .

RTable 2 .Fig. 6 .

 26 Fig. 6. Average value (left) and occurrence time (right) of the first and second βnuc peak for the Wnt-inject model in function of Wnt decay rate (highlighted points correspond to k1 original value as inTable 1 set B).

Fig. 7 .

 7 Fig. 7. The PDF of the time of occurrence of the first βnuc peak (left) and second peak (right) for the Wnt-inject model (with re-injection of 1000 Wnt molecules at t = 450 minutes).

Table 1 .

 1 . Parameter

	initial populations (mol.)		rate constants (mol. • min -1 )
	par. id Set A Set B par. id Set A	Set B par. id Set A	Set B
	nβcyt 11145 12989	k1	0.6	0.27	k6	2.4 • 10 -3 4.48 • 10 -3
	nβnuc 4532	5282	k2	10	20	k7	3 • 10 -4 2.1 • 10 -4
	nAxin 144	252	k3	0.03	0.03	k8	420	600
	nAxinP 125	219	k4	0.03	0.03	k9 1.13 • 10 -4 1.13 • 10 -4
	nWnt 1000	1000	k5 4.48 • 10 -3 4.48 • 10 -3 k10	0.0549	0.0549
			k11	0.135	0.135	k12	2 • 10 -4	4 • 10 -4

sets for stochastic interpretation of Wnt/β-catenin pathway model given by equations

[START_REF] Ajmone Marsan | Modelling with Generalized Stochastic Petri Nets[END_REF]

notice that βnuc dependent Axin expression and AxinP enhanced βcyt degradation determine, de facto, a negative feedback loop between β-catenin and the Axin protein.

Notice that because of the "initial-nondeterminism" of LHA there can be at most one initial state for the product process.

results computed with confidence level 99% and interval-width of 1% of the estimated measure.