Almudena Cánovas-Molina

Monica Montefalcone

Mauatassem Giorgio Bavestrello

Billah Masmoudi

Anne Haguenauer

Pachka Hammami

Lamya Chaoui

Mohammed Hichem Kara

Didier Aurelle
email: didier.aurelle@univ-amu.fr

Mauatassem Billah Masmoudi

From depth to regional spatial genetic differentiation of Eunicella cavolini in the NW Mediterranean

Keywords: Eunicella cavolini, gorgonian, connectivity, spatial scale, Mediterranean, genetic structure

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

Introduction

In a context of growing pressures affecting the marine environments, estimating connectivity among populations is essential since the scale of larval dispersal is crucial for understanding species reorganization and populations evolution after large-scale disturbances [START_REF] Bellwood | Confronting the coral reef crisis[END_REF][START_REF] Salm | Marine protected area planning in a changing climate. Coral reefs and climate change: science and management[END_REF]. Size and interconnectivity of populations are important parameters to estimate for conservation strategies, since small and isolated populations are vulnerable to inbreeding depression, limited spread of beneficial alleles [START_REF] Alleaume-Benharira | Geographical patterns of adaptation within a species' range: interactions between drift and gene flow[END_REF]. and loss of genetic diversity, which might reduce their evolutionary potential and increase their risk of extinction [START_REF] Saccheri | Inbreeding and extinction in a butterfly metapopulation[END_REF]. Population genetics describes connectivity by comparing allele frequencies among spatially discrete subpopulations. High levels of genetic similarity between populations can indicate gene flow over time, usually through larval dispersal, whereas significant differentiation between populations can correspond to significant and persistent barriers to larval exchange. Nevertheless, the observed level of genetic differentiation depends on the interaction between gene flow, genetic drift and mutation (depending on the time scale or molecular markers).

For most benthic marine species with complex life cycles, the exchange of individuals among populations occurs primarily through the pelagic larval stage [START_REF] Cowen | Larval dispersal and marine population connectivity[END_REF]. These exchanges depend on biotic and abiotic factors including spawning outputs, larval dispersal, habitat availability, hydrodynamics and trophic interactions [START_REF] Padrón | Modelling the effect of demographic traits and connectivity on the genetic structuration of marine metapopulations of sedentary benthic invertebrates[END_REF][START_REF] Pineda | Fine-scale genetic structure and inferences on population biology in the threatened Mediterranean red coral, Corallium rubrum[END_REF].

Evidence from direct and indirect approaches using geochemical and genetic techniques suggests that the genetic structure of marine populations, while often low, ranges from fully open to fully closed [START_REF] Cowen | Larval dispersal and marine population connectivity[END_REF].

Understanding the scale of population connectivity is fundamental for the development of suitable conservation strategies and for the design of marine reserves [START_REF] Almany | Connectivity, biodiversity conservation and the design of marine reserve networks for coral reefs[END_REF][START_REF] Palumbi | Marine reserves and ocean neighborhoods: the spatial scale of marine populations and their management[END_REF][START_REF] Shanks | Propagule dispersal distance and the size and spacing of marine reserves[END_REF]. Decisions predicated on the assumption that larvae disperse at long distance may lead to erroneous actions if this assumption is wrong [START_REF] Cowen | Connectivity of marine populations: open or closed?[END_REF]. Studies regarding connectivity conducted in the Mediterranean Sea showed the existence of different patterns of genetic differentiation, ranging from complete panmixia to an east-west or north-south genetic partitioning. These findings warn against any generalization based on the presently available information and strongly indicate the need to obtain genetic connectivity data from a wider array of species, including keystone species such as habitat formers as well as rare or threatened species [START_REF] Paterno | Population genomics meet Lagrangian simulations: Oceanographic patterns and long larval duration ensure connectivity among Paracentrotus lividus populations in the Adriatic and Ionian seas[END_REF].

Among the benthic sessile invertebrates of the Mediterranean, gorgonian octocorals play an important role, contributing to the enhancement of the overall structural complexity, biomass, and species diversity of habitats by their three-dimensional growth pattern [START_REF] Ballesteros | Mediterranean coralligenous assemblages: a synthesis of present knowledge[END_REF][START_REF] Gili | Study and cartography of the benthic communities of Medes Islands (NE Spain)[END_REF][START_REF] True | Étude quantitative de quatre peuplements sciaphiles sur substrat rocheux dans la région marseillaise[END_REF]. Furthermore, organisms living in habitats characterized by the presence of gorgonian corals may experience a sort of buffer zone where environmental alterations occur more slowly and within narrower ranges with respect to the surrounding 'unforested' environment [START_REF] Ponti | Ecological shifts in Mediterranean coralligenous assemblages related to gorgonian forest loss[END_REF].

The slow population dynamics of gorgonians makes them susceptible to a wide range of direct or indirect anthropogenic stressors: mechanical damage (mainly caused by fishing and unregulated recreational diving), pollution and mucilaginous algal aggregates [START_REF] Bavestrello | Damage by fishing activities in the Gorgonian coral Paramuricea clavata in the Ligurian Sea[END_REF][START_REF] Coma | Long-term assessment of the patterns of mortality of a temperate octocoral in protected and unprotected areas: a contribution to conservation and management needs[END_REF][START_REF] Mistri | Effects of a mucilage event on the Mediterranean gorgonian Paramuricea clavata. Short term impacts at the population and colony levels[END_REF], biological invasions [START_REF] Cebrian | Exploring the effects of invasive algae on the persistence of gorgonian populations[END_REF], bacterial pathogens and mass mortality events linked to thermal anomalies [START_REF] Cerrano | A catastrophic mass-mortality episode of gorgonians and other organisms in the Ligurian Sea (North-western Mediterranean), summer 1999[END_REF][START_REF] Garrabou | Mass mortality in NW Mediterranean rocky benthic communities: effects of the 2003 heat wave[END_REF]. For example, during the 2003 mortality event, the yellow gorgonian Eunicella cavolini was affected from the coasts of Provence in south-east France to the gulf of Naples in south Italy, with a percentage of affected colonies reaching 50.8% [START_REF] Garrabou | Mass mortality in NW Mediterranean rocky benthic communities: effects of the 2003 heat wave[END_REF]. This species, along with other anthozoans, is also impacted by other human pressures such as, for example, the sideeffects of fisheries (Cánovas-Molina et al., 2016a;[START_REF] Topçu | Composition and abundance of octocorals in the Sea of Marmara, where the Mediterranean meets the Black Sea[END_REF]. Considering these threats and the need to extend the knowledge on population dynamics and connectivity for conservation and ecosystem-based management purposes, several studies have been performed on the genetic structure of Mediterranean species from the coralligenous assemblages, most of them on octocorals. These previous genetic studies revealed strong population genetic structure potentially because of short-distance larval dispersal for some coralligenous species such as the red coral Corallium rubrum [START_REF] Abbiati | Proceedings of the 1st Mediterranean Symposium on the Conservation of the Coralligenous and Other Calcareous Bio-Concretions[END_REF][START_REF] Boavida | A well-kept treasure at depth: precious red coral rediscovered in Atlantic deep coral gardens (SW Portugal) after 300 years[END_REF]Costantini et al., 2007a, b;[START_REF] Ledoux | Genetic survey of shallow populations of the Mediterranean red coral Corallium rubrum (Linnaeus, 1758): new insights into evolutionary processes shaping nuclear diversity and implications for conservation[END_REF], the red gorgonian Paramuricea clavata [START_REF] Arizmendi-Mejía | Combining genetic and demographic data for the conservation of a Mediterranean marine habitat-forming species[END_REF]Mokhtar-Jamaï et al., 2011), the yellow gorgonian E. cavolini [START_REF] Masmoudi | Contrasted levels of genetic diversity in a benthic Mediterranean octocoral: Consequences of different demographic histories?[END_REF], the only symbiotic Mediterranean gorgonian Eunicella singularis (Constantini et al., 2016;[START_REF] Pey | Thermal threshold and sensitivity of the only symbiotic Mediterranean gorgonian Eunicella singularis by morphometric and genotypic analyses[END_REF] or the sponge Crambe crambe [START_REF] Duran | Strong population structure in the marine sponge Crambe crambe (Poecilosclerida) as revealed by microsatellite markers[END_REF]. Regarding E. cavolini, the first molecular study gave little information on the genetic structure of this species because of a lack of informative markers [START_REF] Calderón | Finding the relevant scale: clonality and genetic structure in a marine invertebrate (Crambe crambe, Porifera)[END_REF]. The use of microsatellite markers in E. cavolini allowed the identification of different genetic clusters in the Mediterranean corresponding to the main geographical regions, and with different levels of genetic diversity: the allelic richness of eastern Mediterranean populations (i.e. from Turkey) was two-fold lower than in Algeria [START_REF] Masmoudi | Contrasted levels of genetic diversity in a benthic Mediterranean octocoral: Consequences of different demographic histories?[END_REF]. The genetic structure in the region of Marseille (France) indicated a differentiation between most populations [START_REF] Masmoudi | Contrasted levels of genetic diversity in a benthic Mediterranean octocoral: Consequences of different demographic histories?[END_REF]. Although no differentiation was observed between depths (between 20 m and 40 m depth) at a given site [START_REF] Masmoudi | Contrasted levels of genetic diversity in a benthic Mediterranean octocoral: Consequences of different demographic histories?[END_REF], significant differences in thermo-tolerance were demonstrated at the bathymetrical scale [START_REF] Pivotto | Highly contrasted responses of Mediterranean octocorals to climate change along a depth gradient[END_REF].

In the present work, we focused on the yellow gorgonian Eunicella cavolini (Koch, 1887), an octocoral species widely distributed in the Mediterranean, which can be observed from less than 10 m to more than 100 m depth [START_REF] Sini | The yellow gorgonian Eunicella cavolini: demography and disturbance levels across the Mediterranean Sea[END_REF]. The genetic structure of this species has been studied mainly at the basin-wide Mediterranean spatial scale [START_REF] Masmoudi | Contrasted levels of genetic diversity in a benthic Mediterranean octocoral: Consequences of different demographic histories?[END_REF]. Definition of the boundaries of Marine Protected Areas (MPAs) should consider the connectivity among habitats for management and conservation of marine species [START_REF] Palumbi | Marine reserves and ocean neighborhoods: the spatial scale of marine populations and their management[END_REF], and this requires a more precise study of connectivity at both the regional and the local spatial scales. In this sense, NW Mediterranean offers an interesting context for genetic structure studies at different spatial scales where the yellow gorgonian forms populations from 15 m depth reaching here the maximum densities of the Mediterranean sites studied by [START_REF] Sini | The yellow gorgonian Eunicella cavolini: demography and disturbance levels across the Mediterranean Sea[END_REF]. The absence or low genetic differentiation between shallow and deep E. cavolini observed at Marseille (SE France) is of great importance considering the impact of climate change and the possibility of exchanges between depths with different thermal conditions [START_REF] Pivotto | Highly contrasted responses of Mediterranean octocorals to climate change along a depth gradient[END_REF], and it should be tested in other areas as well. Located 300 km apart in the same Ligurian Sea, Portofino (NW Italy) offers a perfect location for the genetic study of Eunicella cavolini populations at regional, local and depth scale.

Given the lack of information on the effects of the spatial scale on genetic structure of this species, our aims were: i) to examine the genetic differentiation of E. cavolini at the regional scale in the Ligurian Sea, by focusing on the differences between two adjacent regions 300s kilometres apart, i.e. Marseille in SE France and Portofino in NW Italy, where a number of MPAs have already been established; ii) to analyze the genetic differentiation at the local scale 10s kilometres apart, i.e. within the areas of Marseille and Portofino, in order to acquire insights into the dispersal abilities of this species at small scale;

iii) to investigate the vertical connectivity between populations by including samples from different depths. We tested whether the absence of a vertical genetic differentiation previously observed in Marseille can be extended in other areas of the Ligurian Sea.

Materials and methods

Samples collection

Colonies of Eunicella cavolini were sampled by scuba divers removing a 3 cm branch (usually terminal) using sea-snips (n=628 colonies in total; see Tab. 1). In order to have colonies separated by few and hundreds of km and at different depths, samples were collected few kilometres apart at five sites in the MPA of Portofino (Genoa, NW Italy) and at ten sites in the area around Marseille (SE France) (Fig. 1). Samples around Marseille were taken at two depths (20 and 40 m) in Riou, Veyron and Méjean sites and at one depth in the other seven sites. Nine colonies were sampled on artificial reefs (REE), which were installed in 2007-2008. Samples at Portofino were collected at two depths (20 m and 40 m), except for Punta Cervara (PAR), where E. cavolini only occurs at 13-18 m depth. Geographical distances between sites in the Marseille area ranged from 2 km to 100 km whereas in the Portofino area they ranged from 1 km to 6 km (Fig. 1). Distances between Marseille and Portofino ranged from 381 km to 274 km. Sample size was even in all sites, around

DNA extraction and microsatellite genotyping

Total genomic DNA was extracted using QIAamp ® DNA Mini Kit (Qiagen ® ) according to manufacturer's instructions. Individuals were genotyped at six microsatellite loci: C21, C30, C40, S14, EVER007, EVER009 [START_REF] Holland | Isolation and characterization of fourteen microsatellite loci from the endangered octocoral Eunicella verrucosa (Pallas 1766)[END_REF]Molecular Ecology Resources Primer Development Consortium et al., 2010;[START_REF] Pivotto | Highly contrasted responses of Mediterranean octocorals to climate change along a depth gradient[END_REF]. The loci were amplified through PCR Multiplex using the Type-it ® Microsatellite PCR Kit (Qiagen ® ). PCR amplifications were done in a final reaction volume of 10 μL containing 1 μL of DNA, 5 μL of MasterMix, 3 μL of H2O Milli-Q® and 1 μL of primer mix (with each primer at 2 µM).

Electrophoresis was carried out at the IMBE laboratory on an ABI 3130 genetic analyzer and at the Genotyping and Sequencing facility at Montpellier (Plateforme Génotypage -Séquençage, Université de Montpellier) on an ABI 3500xL Genetic Analyzer (Applied Biosystems), using

GeneScan LIZ 600 (Applied Biosystems) as the internal size standard. Four reference individuals were included in each migration for cross-validation. Electropherograms were analysed with GeneMapper ® version 4.0 (Applied Biosystems).

Genetic diversity and structure

We used MICRO-CHECKER v.2.2.3 (Van Oosterhout et al., 2004) to check for scoring errors owing to stutters, large allele dropout and to estimate null allele frequencies. We used GENCLONE 2.0 (Arnaud-Haond and [START_REF] Arnaud-Haond | GENCLONE: a computer program to analyse genotypic data, test for clonality and describe spatial clonal organization[END_REF] to test for the occurrence of repeated multilocus genotypes (MLGs). Repeated MLGs were omitted from the following analyses.

Linkage disequilibrium was tested among all pairs of loci at each site with a permutation test (n = 1000) using GENETIX v.4.05 [START_REF] Belkhir | GENETIX 4.05, Logiciel Sous Windows TM Pour La Génétique Des Populations. Laboratoire, Génome, Populations, Interactions[END_REF]. The f estimator of FIS [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] was computed for each population using GENETIX. Departures from panmixia were tested with an exact implemented in GENEPOP 4.5.1 with default parameters [START_REF] Rousset | Genepop'007: a complete reimplementation of the Genepop software for Windows and Linux[END_REF].

We analysed the genetic diversity for each population by computing observed (Ho) and [START_REF] Nei | Analysis of gene diversity in subdivided populations[END_REF] unbiased expected heterozygosity (He) with GENETIX. Allelic richness [Ar (g)] and private allelic richness [Ap (g)] were estimated with a rarefaction procedure using the HP-RARE software [START_REF] Kalinowski | HP-RARE 1.0: a computer program for performing rarefaction on measures of allelic richness[END_REF] with the minimum number of genes set to 18. FreeNA was used to estimate FST for each pair of populations using the ENA correction described in [START_REF] Chapuis | Microsatellite null alleles and estimation of population differentiation[END_REF], to provide more precise estimates of FST considering the presence of null alleles. As a complementary measure of genetic differentiation, we computed the estimator of actual differentiation Dest [START_REF] Jost | GST and its relatives do not measure differentiation[END_REF] with the online Software for the Measurement of Genetic Diversity [START_REF] Crawford | SMOGD: software for the measurement of genetic diversity[END_REF]SMOGD 1.2.5). Exact tests of genic differentiation were performed with GENEPOP with default parameters.

We estimated contemporary connectivity among populations with BAYESASS 3.0 [START_REF] Wilson | Bayesian inference of recent migration rates using multilocus genotypes[END_REF]. This software was used to assess the proportions of individuals assigned to their hypothesized population and to the other populations. As the S14 locus showed strong departures from panmixia, we performed this analysis without this locus. The current implementation of BAYESASS is limited to 20 populations. We therefore performed two analyses: one for Portofino with 9 populations, and another for Marseille with 13 populations. The added benefit of this method over a standard genetic distance method is the ability to estimate the direction of gene flow.

To ensure consistent and accurate estimates, we varied seed numbers and explored the number of runs needed (i.e. burn-in) to be discarded before data collection. The burn-in length was set after log-likelihood values peaked. We also varied seed numbers and mixing parameters to get acceptance rates between 40% and 60% [START_REF] Wilson | Bayesian inference of recent migration rates using multilocus genotypes[END_REF]. Delta values of m = 0.05, P = 0.15, and F = 0.3 yielded an average number of changes in the accepted range. For each analysis, we performed five independent (i.e. with different seeds) replicate runs (their average value was then used) of the algorithm for 500 000 iterations, and a burn-in of 200 000 generations was determined an appropriate number of runs for convergence and stabilization of posterior probabilities with a sampling frequency of 2000. As a complementary analysis to BAYESASS, we used the GENECLASS2 software for the assignment of individuals [START_REF] Piry | GENECLASS2: A software for genetic assignment and first-generation migrant detection[END_REF]. We used the criterion of [START_REF] Rannala | Detecting immigration by using multilocus genotypes[END_REF], and the computation of the probability according to [START_REF] Paetkau | Genetic assignment methods for the direct, real-time estimation of migration rate: a simulation-based exploration of accuracy and power[END_REF] with 1000 simulated individuals, and a type I error of 0.01. The colonies sampled on artificial reefs in Marseille (REE) were also analyzed separately with GENECLASS2:

we used the other populations from Marseille as reference samples to study the potential origin of colonies on the reefs. Considering the strong departure from panmixia observed for the S14 locus,

we did not use it with GENECLASS2.

The Bayesian method implemented in STRUCTURE v.2.2 [START_REF] Falush | Inference of population structure using multilocus genotype data: Linked loci and correlated allele frequencies[END_REF][START_REF] Falush | Inference of population structure using multilocus genotype data: dominant markers and null alleles[END_REF][START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF] was used to evaluate the number of genetic clusters (K) in our data set from individual genotypes without prior information on their geographical locations. Again, as the S14 locus showed strong departures from panmixia, we repeated the analysis without this locus.

STRUCTURE was launched using the admixture model, with correlated allele frequencies among clusters and the recessive allele option to cope with null alleles [START_REF] Falush | Inference of population structure using multilocus genotype data: dominant markers and null alleles[END_REF]. Ten independent runs were performed for each K using 500000 iterations and a burn-in period of 200000. Following the recommendations of [START_REF] Jakobsson | Genotype, haplotype and copy-number variation in worldwide human populations[END_REF] and [START_REF] Rosenberg | Genetic structure of human populations[END_REF] for cases of large datasets, we first studied the results of STRUCTURE at several small values of K during a first round of analysis, and then we looked for additional substructure on partitioned datasets during a second round. For the first round, the whole dataset (622 individuals) was used. K was set to vary between 1 and 6. As the analysis of the whole dataset indicated the presence of two clusters corresponding to the two sampling regions, we performed a second analysis by separating the Marseille sites from the Portofino sites. The choice of the best K values was based on the result of STRUCTURE HARVESTER [START_REF] Earl | STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF] with the Evanno's criterion (Evanno et al., 2005), and we retained the values which seemed the most informative, i.e. those with the lowest levels of individual admixture: we therefore discussed K = 2 and 3 for the global analysis and K = 2 for the analysis within regions. The merging of the different runs and the visualization of the results were done with the CLUMPAK software [START_REF] Kopelman | Clumpak: a program for identifying clustering modes and packaging population structure inferences across K[END_REF].

We tested for isolation by distance (IBD) by examining the correlation between pairwise FST (1-FST) values and the logarithm of the geographical distances [Ln(d)] between populations [START_REF] Rousset | Genetic differentiation and estimation of gene flow from FStatistics under isolation by distance[END_REF] with IBD Web Software (IBDWS 3.16;[START_REF] Jensen | Isolation by distance, web service[END_REF]. Geographical distances were measured using GOOGLE EARTH 6.2.2.6613 (available at https:// www.google.fr/intl/fr/earth/), considering the minimum distances by sea between sites. IBD was assessed independently within Marseille and Portofino, as well as across the entire Liguro-Provençal area in the Ligurian Sea. The significance of the correlation between the two distance matrices was tested using a Mantel test (n = 1000) in IBD Web Software.

Whenever necessary, corrections for multiple tests were performed with the False Discovery Rate (FDR) method of [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF].

Results

Genetic variability

All loci were polymorphic in all populations except for locus EVER007 in PFD and C30 in PTD.

The total number of alleles ranged between 4 for EVER007 and 16 for S14, with a mean of 10.5 alleles per locus. Observed heterozygosity varied from 0 for EVER007 in PFD and C30 in PTD to 0.86 for EVER009 in PAR and PFD. Expected heterozygosity ranged from 0 for EVER007 in PFD and C30 in PTD to 0.86 for S14 in VED. MICRO-CHECKER suggested the presence of null alleles for S14, as indicated by the general excess of homozygotes for most allele size classes.

Global significant linkage disequilibrium (LD) between the loci C21 and S14 was detected considering all samples only (P < 0.05). When considering each population separately, significant LD among pairs of loci was observed in nine populations, but these were not the same pairs of loci in the different populations. After FDR correction, significant LD was only detected in VES between loci EVER009 -C30. Seven cases of duplicated MLGs were observed (Table S1). In two cases, the repeated MLGs were observed in different populations (PTD vs SGD, MJS vs CAV).

Population genetic structure

The mean global FST was 0.07 (± 0.04), the mean FreeNA corrected FST was 0.07 (± 0.03), and the mean Dest was 0.17 (± 0.08). Pairwise multilocus FST between all populations ranged from 0.002 between MJS and MEJ to 0.166 between PAR and SAM (Tab. S3). Pairwise FreeNA corrected FST ranged from 0.01 (PFS vs PFD, CSS vs CSD, CSD vs SGD, RIS vs RID, VES vs SOM, MEJ vs SAM) to 0.17 (PAR vs SAM). Pairwise Dest ranged from 0 (VED vs VES) to 0.38 (CSD vs POU) (Tab. S4). After the FDR correction, all pairwise exact tests of genic differentiation were significant apart from CSS vs CSD and RIS vs RID (two comparisons between depths; Tab. S3).

Apart from depth, significant comparisons were found with a minimum distance of 0.99 km in Portofino and 2.95 km in Marseille. When locus S14 was omitted, two additional differentiation tests were non-significant: VES vs SOM and PFS vs PFD.

The correlation between FST (1 -FST) and Ln(dist) was significant using geographical distances (rLn(dist) = 0.596; P< 0.001), confirming the occurrence of an IBD model of gene flow (Fig. 2).

Within areas, correlations were also significant for the Portofino (rLn(dist) = 0.625; p < 0.001) and Marseille areas (rLn(dist) = 0.436; p < 0.001) (Fig. 2B and2C). The IBD slope was similar in the two areas (0.019 in Portofino and 0.015 in Marseille). When using FST corrected by FreeNA, IBD was significant at the global level (p < 0.001), within Marseille (p = 0.04) and within Portofino (p = 0.002). Again, the IBD slope was similar in the two areas (0.032 at Portofino and 0.035 at Marseille; data not shown). The results of the genetic clustering using STRUCTURE indicated that the highest Delta (K) value for the whole dataset was 2 (data not shown). With K = 2, the clustering showed a clear distinction between Marseille and Portofino, and K = 3 did not provide additional information (Fig. 3). The addition of locus S14 to the analysis gave a very similar result (data not shown). Each of these two regional clusters was then subjected to a second round of STRUCTURE. The highest Delta (K)

value was 3 at Marseille and 4 at Portofino. Nevertheless, inside each region, the results did not indicate the presence of any additional cluster: in both cases, all individuals were equally shared among clusters for all K (data not shown).

Fig. 3. Population structure as inferred by the STRUCTURE analysis. Each individual is represented by a vertical line partitioned into K-coloured segments that represent the individual's membership fraction in K clusters. Each population is delineated by black vertical lines and named as in Table 1.

Migration analysis using the BAYESASS analysis indicated similar levels of self-recruitment in the two regions: the proportion of individuals derived from the source population varied between 0.69 (PFD) and 0.80 (PAR) at Portofino, and between 0.68 (CAV) and 0.80 (POU) at Marseille (Tab.

S5

). The estimates of migration rates between populations were generally low and with quite large standard errors. The assignment analysis of GENECLASS2 indicated that the probability of assignment (averaged over all individuals of the population) was higher for populations from their region of origin than for another region: the intra-region assignment probability was 0.31 (± 0.10)

for Portofino and 0.32 (± 0.10) for Marseille (Table S6). The within-region highest assignment probabilities corresponded to the population where individuals were sampled: these retention probabilities were very similar among populations and regions: 0.46 (± 0.02) at Portofino and 0.45 (± 0.03) at Marseille. The assignment probability of Portofino individuals to Marseille populations was 0.11 (± 0.07) and 0.16 (± 0.07) in the other direction. The assignment probabilities of E.

cavolini colonies sampled on artificial reefs near Marseille did not allow the identification of a single main origin, and for a given colony the probabilities were often similar among different candidate populations (Table S7). Only the colony 14B displayed a much higher probability of origin for one population (Saména) than for the other populations. The colonies 15J, 15M, 15O and 15P displayed low assignment probabilities (< 0.10) for all the populations tested here.

Discussion

Genetic diversity

The analysis of six microsatellite loci evidenced low genetic diversity for Eunicella cavolini (He by population over loci = 0.51) when compared with other Mediterranean octocorals: Corallium rubrum studied in the same area (He by population over loci = 0.77; [START_REF] Ledoux | Genetic survey of shallow populations of the Mediterranean red coral Corallium rubrum (Linnaeus, 1758): new insights into evolutionary processes shaping nuclear diversity and implications for conservation[END_REF] and Paramuricea clavata (He by populations over loci = 0.74; Mokhtar-Jamaï et al., 2011). There could be several possible and non-mutually exclusive reasons. i) the microsatellites used here were less variable because of their sequence or origin: S14 and EVER007, and 009 correspond to crossamplification from E. singularis and E. verrucosa respectively; ii) the species is less variable than other octocoral species, likely due to the impact of past bottlenecks; iii) a geographical effect: the comparison of populations from Marseille with other Mediterranean areas indicated that in the NW Mediterranean, the genetic diversity of E. cavolini is lower than in the SW Mediterranean populations, but higher than in the E Mediterranean [START_REF] Masmoudi | Contrasted levels of genetic diversity in a benthic Mediterranean octocoral: Consequences of different demographic histories?[END_REF]. Nevertheless, even in the area of the highest genetic diversity for E. cavolini, the expected heterozygosity was somewhat lower than in the two others Mediterranean octocorals (He = 0.68 in Algeria; [START_REF] Masmoudi | Contrasted levels of genetic diversity in a benthic Mediterranean octocoral: Consequences of different demographic histories?[END_REF]. A comparative study of the demographic history of these species with similar markers (such as Single Nucleotide Polymorphisms, SNPs) could help in understanding the origin of these differences.

4.2 Genetic structure at large spatial scale A pattern of isolation by distance was observed at the regional spatial scale with significant genetic structuring found between the two regions. Values of FST for the Liguro-Provençal geographical area ranged between 0.002 and 0.166, with a mean value of 0.09. If we make comparisons with other octocorals, FST for the Liguro-Provençal area was between 0.07 and 0.14 with a mean value of 0.09 for P. clavata (Mokhtar-Jamaï et al., 2011), and between 0.09 and 0.2 with a mean value of 0.13 for C. rubrum [START_REF] Ledoux | Genetic survey of shallow populations of the Mediterranean red coral Corallium rubrum (Linnaeus, 1758): new insights into evolutionary processes shaping nuclear diversity and implications for conservation[END_REF]. Such comparisons between species can be impacted by different levels of heterozygosities among markers, which influences the highest FST values (Jost et al. 2018). As previously mentioned, Hexp is lower in E. cavolini based on our markers, than in C.

rubrum and P. clavata. Despite this, FST is lower for E. cavolini than for the other species. It thus shows that genetic differentiation within the Liguro-Provençal area is higher for C. rubrum than for E. cavolini and P. clavata. Such differences in the levels of genetic differentiation between species could be due to differences in dispersal abilities or in local effective sizes. Although all brooders are more likely to exhibit restricted gene flow (Carlon et al., 1999), C. rubrum has a mean dispersal ability of 10s of meters or less (Constantini et al., 2007a), whereas for P.clavata it would be around 1 km (Mokhtar-Jamaï et al., 2011) and for E. cavolini our results suggest a mean dispersal ability of less than 1 km. Nevertheless this may be at odd with experimental data on larval phase duration in the red coral (Martinez-Quintana et al., 2015), and this does not preclude the possibility of exchanges at higher distances (Ledoux et al., 2010b), even if their frequency remains unknown.

The structure pattern revealed by STRUCTURE showed a distinction between Portofino and Marseille. This pattern is also supported by the higher values of FST between areas (average FST= 0.092) than within areas (average FST= 0.063 for Portofino and average FST= 0.069 for Marseille),

and by the marked differences in assignment probabilities between regions compared to within regions. These results indicate a limited gene flow between the two regions for the yellow gorgonian, which could be linked to intrinsic limits to dispersal at such distances (the two regions are 273 km apart), and/or to the effect of a barrier to gene flow. The genetic analyses of yellow gorgonian populations in other intermediate positions between Marseille and Portofino, could help in testing whether there is a genetic barrier between the two regions, or if the observed differences are induced by sampling along a geographical gradient of genetic differences (i.e. IBD, [START_REF] Aurelle | Interplay between isolation by distance and genetic clusters in the red coral Corallium rubrum: insights from simulated and empirical data[END_REF]. The asymmetry of assignment probabilities between the two regions suggests a higher gene flow from east to west that is consistent with the main direction of currents in this area [START_REF] Berline | A connectivity-based eco-regionalization method of the Mediterranean Sea[END_REF]. Nevertheless, the confidence intervals are overlapping and more precise estimates would be useful here to test this hypothesis.

Regarding the larval biology of E. cavolini, very little information is available, but it may contribute to the observed genetic structure. Eunicella cavolini is also a brooding species, which implies more reduced larval dispersal abilities. Such low larval dispersal is observed in other coralligenous species, where the studies available to date showed significant differentiation from a few kilometres or tens of cm [START_REF] Abbiati | Proceedings of the 1st Mediterranean Symposium on the Conservation of the Coralligenous and Other Calcareous Bio-Concretions[END_REF]Calderon et al., 2007;[START_REF] Goffredo | Genetic differentiation among populations of the Mediterranean hermaphroditic brooding coral Balanophyllia europaea (Scleractinia: Dendrophylliidae)[END_REF][START_REF] Ledoux | Genetic survey of shallow populations of the Mediterranean red coral Corallium rubrum (Linnaeus, 1758): new insights into evolutionary processes shaping nuclear diversity and implications for conservation[END_REF]. In contrast to these cases, another Mediterranean benthic species, the sea urchin Paracentrotus lividus, displayed a much larger genetic structure [START_REF] Paterno | Population genomics meet Lagrangian simulations: Oceanographic patterns and long larval duration ensure connectivity among Paracentrotus lividus populations in the Adriatic and Ionian seas[END_REF]. It would therefore be interesting to study on more species the ecological or evolutionary factors leading to the selection of different dispersal strategies in these species. Indeed multispecific approaches can reveal common barriers to gene flow and identify cold and hot spots of genetic diversity (Cahill et al., 2017).

Genetic structure at small spatial scale

A pattern of isolation by distance was observed also at the local spatial scale with significant genetic structuring found within the two regions. The populations from each area, Marseille and Portofino, correspond to a single genetic cluster according to STRUCTURE. The additional genetic clusters identified within regions were not informative and corresponded to the limits of inference of this method. Within regions, most populations appeared genetically differentiated except for comparisons between depths (see below). Significant genetic differentiation was found between sites few kilometres apart. The pattern of differentiation within regions did not correspond to different clusters, but to the isolation by distance of differentiated populations. The use of FST and Dest gives different and complementary information on genetic differences between populations (Jost et al., 2018). Here at different scales, DEST estimates where higher than Fst: DEST then indicates an important allelic differentiation between populations, higher than the differentiation estimated from fixation index alone (Jost et al., 2018). At a local scale, and on the basis of DEST, one can identify populations with particularly high differentiation: such as Pouard in Marseille, or Punta Faro Deep in Portofino. These sites did not correspond to particularly distant populations, and the origin of these differences remains to be studied (see below). More generally, in the short term, the interaction between gene flow and genetic drift precludes the homogenization of genetic diversity at the metapopulation scale.

The BAYESASS software was used to estimate contemporary connectivity among sites. The inferred non-migration rates fell mostly within two classes (around 0.67 and 1). Such values correspond to the bounds of the prior distribution and could indicate a bias in the inference [START_REF] Meirmans | Nonconvergence in Bayesian estimation of migration rates[END_REF], despite the use of five independent runs. The Sarranier, 3PP cave and Punta

Cervara populations displayed higher non-migration rates (above 0.9), which would agree with their relative geographical isolation, but this was also observed for a few a priori less isolated populations such as Somlit and Saména. The GENECLASS2 software was then used to obtain another estimate of contemporary migration. With this method, the low probabilities of assignment to the population of origin could be the result of a quite low within region genetic structure. As a comparison, much higher levels of self-recruitment were inferred with this method for P. clavata, including for some sites studied here around Marseille (Mokhtar-Jamaï et al., 2011). This confirms that the dispersal abilities of E. cavolini are higher than those of P. clavata. Within a given region, sporadic gene flow could take place and contribute to limiting the observed differentiation compared to other octocoral species. In this case, the use of additional genetic markers could be tested to obtain higher probabilities of assignment. (2015), who observed that a recently founded population may correspond to several source populations.

A study combining more markers (such as SNPs, [START_REF] Paterno | Population genomics meet Lagrangian simulations: Oceanographic patterns and long larval duration ensure connectivity among Paracentrotus lividus populations in the Adriatic and Ionian seas[END_REF] and integrating the putative larval characteristics with local hydrodynamics, could help to better understand the observed structure and to estimate connectivity patterns (e.g. [START_REF] Thibaut | Connectivity of populations of the seaweed Cystoseira amentacea within the Bay of Marseille (Mediterranean Sea): genetic structure and hydrodynamic connections[END_REF]. In addition, the possibility of hybridization between E. cavolini and the other species of the genus Eunicella (i.e. E. singularis and E. verrucosa) occurring in the same area (Aurelle et al., 2018) could also shape the observed genetic structure and diversity, depending on the density of the different species and on the frequency of genetic incompatibilities. These processes are only beginning to be considered in this context, but would require an in-depth study.

Genetic structure between depths

Over seven comparisons between depths, only in two sites we found no genetic difference with depths. This is partly contradictory with the lack of genetic differentiation observed between depths in the Marseille area [START_REF] Masmoudi | Contrasted levels of genetic diversity in a benthic Mediterranean octocoral: Consequences of different demographic histories?[END_REF], where seven loci (six loci in the present paper) and a permutation test (exact test in the present paper) were used. These differences suggest that the differentiation between depths is quite low, but near significance depending on sampling effort or drift effects, and that the exact test is more powerful than the permutation test. to water stratification or larval properties) interacting with hydrodynamics might explain the differences observed between these two species and E. cavolini. For the congeneric species E.

singularis, a strong restriction to gene flow was observed between 30 m and 40 m depth [START_REF] Costantini | Limited genetic connectivity between gorgonian morphotypes along a depth gradient[END_REF]. We did not observe this for E. cavolini, but enlarging this study over a wider bathymetrical range might be interesting.

The pattern of genetic structure along the depth gradient can be compared with the pattern of adaptation of E. cavolini. Environmental variables, such as light and temperature, vary along the depth gradient and can induce local adaptation [START_REF] Webster | Who should pick the winners of climate change?[END_REF]. Thermo-tolerance experiments showed that the response to thermal stress is very different between depths in E.

cavolini, with a lower thermo-tolerance in deep individuals [START_REF] Pivotto | Highly contrasted responses of Mediterranean octocorals to climate change along a depth gradient[END_REF]. Such differential phenotypic buffering is observed in a context of low and non-significant genetic differentiation for a studied site (e.g. in the site of Riou, and see also [START_REF] Pivotto | Highly contrasted responses of Mediterranean octocorals to climate change along a depth gradient[END_REF]. For E. cavolini, the adaptation to the different thermal regimes could correspond to acclimatization or to a combination of genetic and non-genetic effects. In addition, the reduced number of loci analysed here does not allow us to reject the possibility of genetic adaptation determined by other genomic regions. More generally, E. cavolini can be found in a wide range of light conditions (Cánovas Molina et al., 2016b), which could point to the general physiological plasticity of the species.

The role of deep environments as possible refugia from disturbances for shallow marine species (i.e. the deep refugia hypothesis) has been widely recognized [START_REF] Bongaerts | Assessing the 'deep reef refugia' hypothesis: focus on Caribbean reefs[END_REF]; however, the genetic differentiation of E. cavolini with depth, although low, questions the possibility for deep populations to act as a regular source of new larvae in the event of mass mortality of shallow colonies due to thermal anomalies [START_REF] Cerrano | A catastrophic mass-mortality episode of gorgonians and other organisms in the Ligurian Sea (North-western Mediterranean), summer 1999[END_REF][START_REF] Garrabou | Mass mortality in NW Mediterranean rocky benthic communities: effects of the 2003 heat wave[END_REF]. Moreover the possible adaptation of octocorals to different depths (Ledoux et al., 2015), including in E. cavolini [START_REF] Pivotto | Highly contrasted responses of Mediterranean octocorals to climate change along a depth gradient[END_REF], could influence the evolution of populations along this depth gradient.

Conclusion

Our results confirm the strong genetic structure observed in the Mediterranean octocorals studied so far, which contrasts with what has been observed for other anthozoan species (e.g. Alcyonium digitatum in the northeast Atlantic; [START_REF] Holland | Contrasting patterns of population structure and gene flow facilitate exploration of connectivity in two widely distributed temperate octocorals[END_REF]. The genetic structure in Eunicella cavolini is lower than in Paramuricea clavata and Corallium rubrum, the other two typical benthic species occurring in coralligenous assemblages of the Mediterranean Sea. The high differentiation observed at the large spatial scale of 100s of kilometres, i.e. between the two regions, argues for the importance of local management of such benthic species. For example, a network of connected MPAs focusing on this species, as well as on other species of the coralligenous assemblages with similar dispersal abilities, should ideally consider these spatially restricted dispersal abilities. The weaker but significant genetic structure observed in E. cavolini at the small spatial scale of kilometres, i.e. among sites within the same region, suggests that gene flow can be important in the short-term evolution of populations. The identification of sites with important genetic differentiation could be taken into account in conservation. Regarding the response to climate change, the observed diversity in thermo-tolerance levels is correlated with the genetic structure estimated with microsatellites.

  Commenté [DA1]: c'est un essai, à voir si c'est mieux? Commenté [DA2]: Made?

Fig. 1 .

 1 Fig. 1. Geographical location of the 15 sites in the Ligurian Sea where populations of Eunicella cavolini were sampled and analyzed in the present study. A: ten sampling sites in the area of Marseille (SE France); B: five sampling sites in the area of Portofino (NW Italy). Refer to Tab. 1 for the labels of each site.

Fig. 2 .

 2 Fig. 2. Isolation by distance (IBD) pattern for all populations (A), taking into account the minimum geographical distances between populations. Separate analyses of IBD for Portofino populations (B) and Marseille populations (C). The graphs show the linear regression of the genetic distance measured as FST (1 -FST) over logarithms of the geographical distance (m).

Table 1 .

 1 Information on the samples collection of Eunicella cavolini in Marseille (France)

	Commenté [DA3]: Suggestion: "In the area of Marseille, a new Marine Protected Area has recently been designed. It is then important to study its potential connectivity with other populations and MPAs along the coast.'" (we can also indicate that there are almost no other populations on the west side...)

  Two repeated MLGs were observed on artificial reefs and four MLGs were observed in the CSS Deviations from Hardy-Weinberg equilibrium Significant multilocus heterozygote deficiencies were observed in ten populations: PFS, CSS, CSD, SGS, SGD, PTS, SAR, 3PP, POU and VED (Tab. S2). When not considering the locus S14, significant heterozygote deficiencies was observed only in VED and 3PP (Tab. S2). For EVER007, EVER009 and C21, the null hypothesis of panmixia was rejected for two populations. For C30 and C40, the null hypothesis of panmixia was not rejected for any population.Multilocus FIS values (f) ranged between -0.03 for RID and 0.27 for 3PP, with a mean value equal to 0.10. Considering each locus separately, f values ranged from -0.34 for C40 in SOM to 0.78 for S14 in 3PP (Tab. S2).

	and REE to 0.38 in PFS (Tab. 2). No difference of expected heterozygosity was found between the areas of Portofino (He= 0.52 ± 0.02) and Marseille (He=0.51± 0.02). He f Ar(18) Ap(18) 3.2 Ho SOM 0.46 (0.18) 0.49 (0.24) 0.07 3.65 (1.64) 0.05 (0.12) SAM 0.43 (0.26) 0.46 (0.27) 0.08 3.62 (1.83) 0.07 (0.18) CAV 0.48 (0.18) 0.51 (0.24) 0.08 4.22 (1.68) 0.17 (0.18) MJS 0.51 (0.33) 0.49 (0.33) -0.02 3.99 (2.18) 0.14 (0.23) MEJ 0.48 (0.31) 0.50 (0.34) 0.06 4.25 (2.36) 0.13 (0.31) REE 0.50 (0.29) 0.51 (0.25) 0.06 3.60 (1.55) 0 (0) Mean Value 0.47 (0.04) 0.51 (0.03) 0.10 3.88 (0.24) 0.08 (0.09)
	population. At the population level, expected heterozygosity ranged from 0.58 for VED to 0.46 for SAM. Ar (18) values ranged from 3.41 in PAR to 4.25 in MEJ. Ap (18) values were scattered from 0 in RIS PAR 0.51 (0.26) 0.53 (0.24) 0.04 3.41 (1.59) 0.07 (0.18) PFS 0.48 (0.22) 0.55 (0.27) 0.16 4.16 (1.90) 0.38 (0.55) PFD 0.52 (0.31) 0.54 (0.31) 0.06 3.66 (1.96) 0.01 (0.02) CSS 0.47 (0.11) 0.55 (0.24) 0.15 4.10 (2.14) 0.03 (0.08) CSD 0.45 (0.24) 0.52 (0.29) 0.14 4.06 (2.19) 0.06 (0.15) SGS 0.42 (0.24) 0.49 (0.30) 0.17 4.06 (2.19) 0.03 (0.05) SGD 0.44 (0.31) 0.50 (0.35) 0.13 3.95 (2.18) 0.07 (0.13) PTS 0.42 (0.29) 0.49 (0.33) 0.16 3.78 (2.11) 0.01 (0.02) PTD 0.46 (0.29) 0.47 (0.31) 0.05 3.63 (2.01) 0.22 (0.25) SAR 0.42 (0.20) 0.47 (0.25) 0.13 3.96 (2.08) 0.03 (0.08) 3PP 0.37 (0.27) 0.50 (0.31) 0.27 3.63 (1.56) 0.13 (0.14) RIS 0.50 (0.28) 0.53 (0.32) 0.09 3.96 (2.08) 0 (0) RID 0.55 (0.32) 0.52 (0.31) -0.03 3.86 (1.80) 0.06 (0.12) POU 0.42 (0.30) 0.51 (0.35) 0.22 3.82 (2.01) 0.05 (0.11) VES 0.48 (0.30) 0.50 (0.32) 0.05 3.78 (2.12) 0.03 (0.05) VED 0.52 (0.21) 0.58 (0.20) 0.13 4.16 (2.16) 0.11 (0.13)

Table

2

. Measures of genetic diversity for 22 populations of Eunicella cavolini based on six microsatellites loci. Ho: observed heterozygosity; He: unbiased expected heterozygosity; f:

[START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF]

estimator of multilocus FIS; Ar(18) and Ap(18): allelic and private allelic richness, respectively, with rarefaction for a corresponding sample size of 18. Standard deviations over loci are in brackets. Values in bold are significant at the 0.05 level after FDR correction.

  [START_REF] Padrón | Modelling the effect of demographic traits and connectivity on the genetic structuration of marine metapopulations of sedentary benthic invertebrates[END_REF]. Different levels of genetic drift could induce contrasted levels of genetic structure, either because of different current population sizes or different demographic histories. Regarding the colonization of artificial reefs near Marseille, several populations could correspond to potential sources. The very low assignment probabilities for four colonies indicate that non-sampled populations may have contributed to this new population. This possibility of connection with different populations may seem at odds with the observed differentiation for most populations. Nevertheless, gene flow in populations might be limited by established colonies (e.g. through space limitation), which is not the case in new populations. Recolonization might then not be limited in an area with enough surrounding populations. As similar results was obtained on P. clavata by Arizmendi-Meja et al.

	structure is a pattern shaped by the interaction between physical drives (e.g. currents, barriers) and
	biological processes, such as larval biology and demography (e.g. fluctuations in population density due to habitat fragmentation and recruitment;	Commenté [DA6]: Jost, L., Archer, F., Flanagan, S., Gaggiotti, O., Hoban, S., & Latch, E. (2018). Differentiation measures for conservation genetics. Evolutionary Applications.
	Some of the studied sites could be more isolated because of the distance (e.g. Sarranier and 3PP
	cave), a particular position behind a headland (e.g. Punta Cervara), or caused by oceanographic
	features, which still remain to be investigated (Pouard, Punta Faro Deep). Population genetic

  Two others Mediterranean octocoral species, Paramuricea clavata (Mokhtar-Jamaï et al., 2011) and Corallium rubrum (Constantini et al., 2011; Ledoux et al., 2010) showed a generalized genetic differentiation with depths in the area of Marseille. Intrinsic biological differences (e.g. spawning date compared

Commenté [DA4]: Evanno, G.,Regnaut, S., & Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular ecology, 14(8), 2611-2620.

Commenté [DA7]: Arizmendi-Mejía, R., Linares, C., Garrabou, J., Antunes, A., Ballesteros, E., Cebrian, E., ... &Ledoux, J. B. (2015). Combining genetic and demographic data for the conservation of a Mediterranean marine habitatforming species. PLoS One, 10(3), e0119585.

Acknowledgments

We are grateful to Alice Oprandi, Simone Musumeci (DiSTAV, Università degli Studi di Genova), and F. Zuberer (Centre d'Océanologie de Marseille) for field collection support. We thank Sandrine Ruitton for the indication of the sampling of E. cavolini on artificial reefs. We thank the molecular biology service of the IMBE. This study was funded by the Università degli Studi di Genova (A. Cánovas-Molina PhD fellowship), by the project PRIN 2010-11 Bioconstructions funded by the Italian Ministero dell'Istruzione, dell'Università e della Ricerca and by the Algerian -French collaboration program Tassili n°12MDU853. The project leading to this publication has received funding from European FEDER Fund under project 1166-39417 I state in behalf of all co-authors, that no conflict of interest exists in this work.

Supplementary data:

Supplementary material is available at the online version of the manuscript.

S1


An Ha Truong
email: truonganha87@gmail.com

Piera Patrizio

Sylvain Leduc

Florian Kraxner

Minh Ha-Duong

To appear in Journal of Cleaner Production Reducing emissions of the fast growing Vietnamese coal sector: the chances offered by biomass co-firing

Keywords: Co-firing, emission reduction, bioenergy, rice residues, greenhouse gas emissions, spatial explicit exploration; Highlights

Vietnam's Power Development Plan 7A authorized many new coal power plants projects, implying an increase of greenhouse gases emissions from 90 MtCO2eq/year today to 360 MtCO2eq/year in 2030. How could co-firing technology -that is the partial substitution of coal by biomass-contributes to mitigate that problem? In this study, we assess the costs and potentials of co-firing rice residues in present and planned coal power plants in Vietnam using a spatially explicit optimization model: BeWhere, adapted as recursive annual dynamic. We found that, the cost of CO2 emissions is the key parameter determining at what level the technology is used. A cost of CO2 emissions of 8 $/tCO2 mobilizes the maximum technical potential of the rice straw and husk domestic resource, with an annual emission reduction of 28 MtCO2eq/year by 2030. At this level, biomass co-firing contributes to an 8% emission reduction in the coal power sector with the abatement cost of 137 Million USD.

Introduction

In 2018, coal contributes for about 40% of world's electricity production. However, global trend of coal power plant development is going downhill with decline of 0.6% for 2014 and 2.8% for 2015 [1]. This trend reflects the rise of renewable energy in global energy agenda, not only because of the climate goals set by almost all nations in Paris Agreement but also because of the increasing competitiveness of renewable energy technologies. Vietnam committed to reduce its emissions by 8% in 2030 compared to Business As Usual (BAU) and up to 25% with international support as its Nationally Determined Contribution (NDC) [START_REF] Monre | Intended Nationally Determined Contribution of Vietnam[END_REF]. The emission projection in Vietnam's NDC is 784 MtCO2eq in 2030 in BAU scenario. The Vietnam Green Growth Strategy also set target for energy activities to reduce GHG emission by 20% to 30% by 2030 compared to BAU. These ambitions show the willingness of the Vietnamese government to fight climate change. Nevertheless, these ambitions are not compatible with the vision of National Power Development Plan (PDP7A) approved by the Government in 2016.

In 2017, twenty coal power plants were in operation in Vietnam with a total installed capacity of 13 GW [START_REF]Energy Statistics Vietnam[END_REF]. Together they provided 33.5% of total electricity supply [START_REF] Evn | Vietnam electricity annual report[END_REF]. By 2030, the current PDP 7A authorizes to add 43 GW of coal-based generation capacity [START_REF]The adjustments of 7th Vietnam Power Development Planning[END_REF]. According to this plan, new coal power plants are set to open every year from 2018 until 2030. More than 50% of total electricity production could come from coal-fired stations by 2030. The main argument for increasing coal in power mix is to meet the power demand of a fast-growing economy at a reasonable price.

During the past decade, Vietnam has experienced the highest increase of greenhouse gas emission among Association of South East Asian Nations (ASEAN) countries: total GHG emission and emission per capita increased three folds in ten years while carbon intensity over GDP increased by 48% [START_REF]Energy Statistics Vietnam[END_REF]. The Initial Biennial Updated Report of Vietnam to UNFCCC [START_REF] Monre | The initial biennal updated report of Viet Nam to the UNFCCC[END_REF] projects that the GHG emissions of the energy sector will rise from 141 MtCO2eq in 2010 to 381 MtCO2eq in 2020 and to 648 MtCO2eq in 2030. This sector accounts for more than 85% of total national GHG emissions. In 2014, total emissions of all operating coal-based power plants in Vietnam is estimated at 36.8 Mt [START_REF] Nguyen | Determination of GHG control measures in coal power sector and proposed implementation pathway[END_REF]. With the increase of coal capacity, GHG emissions from coal power sector is estimated at 176.9 Mt in 2020 and 490.6 Mt in 2030 [START_REF] Nguyen | Determination of GHG control measures in coal power sector and proposed implementation pathway[END_REF], which is an order of magnitude higher than the 2014 level.

Setting a national emission reduction target at 25% in a context where the planned emissions from coal increase ten times raises a tough challenge. These coal emissions account for 63% total national GHG emissions in BAU scenario in 2030. This issue of balancing between electricity supply and climate change goal is under the umbrella of well-known energy trilemma concept, which defined by the World Energy Council [START_REF]Policies for the Future, 2011 Assessment of Country Energy and Climate Policies[END_REF] as "the conflicting goals that governments face in securing energy supplies, providing universal energy access and promoting environmental protection". Therefore, it is urgent to assess all kind of technologies that could help Vietnam to harmonize the emission target with the power development plan and work towards the sustainable development of the country.

Biomass co-firing has long been recognized as a way to reduce the greenhouse gas emission from coal burning power plants [START_REF] Agbor | A review of biomass co-firing in North America[END_REF][START_REF] Sahu | Coal-biomass co-combustion: An overview[END_REF][START_REF] Andrić | Environmental performance assessment of retrofitting existing coal fired power plants to co-firing with biomass: carbon footprint and emergy approach[END_REF][START_REF] Morrison | Life cycle assessment of co-firing coal and wood pellets in the Southeastern United States[END_REF] and as the most cost-effective approach to biomass utilization by the electric utility industry [START_REF] Agbor | A review of biomass co-firing in North America[END_REF][START_REF] Sahu | Coal-biomass co-combustion: An overview[END_REF][START_REF] Andrić | Environmental performance assessment of retrofitting existing coal fired power plants to co-firing with biomass: carbon footprint and emergy approach[END_REF][START_REF] Morrison | Life cycle assessment of co-firing coal and wood pellets in the Southeastern United States[END_REF]. Co-firing can utilize the existing coal power plants equipment, such as boilers and fuel feeding chains. This makes investment costs for cofiring the lowest among biomass power technologies [START_REF]IRENA. Biomass for Power Generation[END_REF]. Moreover, co-firing has higher conversion efficiency than dedicated biomass-only power plants. Efficiency of biomass combustion to electricity in co-firing plants could be 33% -45% [START_REF] Agbor | A review of biomass co-firing in North America[END_REF][START_REF] Sahu | Coal-biomass co-combustion: An overview[END_REF][START_REF]Biomass for Power Generation and CHP[END_REF], depending on plant size, technology and biomass share and quality, while it is only about 30% in dedicated biomass power plant [START_REF]Biomass for Power Generation and CHP[END_REF]. With all these advantages, co-firing is an attractive option in Vietnam where coal-based electricity generation expands rapidly and a large biomass resource base is available.

Co-firing is a mature technology and has been applied at industrial scale in many countries [START_REF] Sahu | Coal-biomass co-combustion: An overview[END_REF][START_REF] Irena | Biomass Co-firing-Technology Brief[END_REF]. The technology, however, is yet to be used much in Vietnam. The PDP 7A mentions co-firing biomass with coal but does not state clearly the roadmap to develop co-firing technology in the country. Aware of domestic coal shortage in near future, the government has funded a research project in pilot scale to mix imported coal with domestic coal in a coal power plant [START_REF] Nguyen | Study on co-firing domestic coal with imported coal to improve the fuel utilization efficiency in coal power plants in Vietnam[END_REF]. Nevertheless, there is a lack of local expertise on co-firing biomass with coal. To this end, the aim of our study is to explore quantified scenarios for implementing co-firing at the national scale, considering the cost and greenhouse gas mitigation. To our best knowledge, this study is the first attempt to evaluate co-firing in Vietnam at national scale.

Studies show that co-firing using domestic biomass resources provide a better results in emission reduction compared to using imported resources [START_REF] Morrison | Life cycle assessment of co-firing coal and wood pellets in the Southeastern United States[END_REF][START_REF] Tsalidis | Life cycle assessment of direct co-firing of torrefied and/or pelletised woody biomass with coal in The Netherlands[END_REF]. As Vietnam has substantial biomass potential, in this study we look only at domestic biomass resource as co-firing fuel. An estimation by [START_REF] Arvo | Development of biomass fuel chains in Vietnam[END_REF] shows that total solid biomass fuel potential in Vietnam by 2010 was 1 346 PJ, about half of total energy consumption of the country in 2009. According to this estimation, agriculture produced almost 70% of the total solid biomass fuel potential and the rest came from firewood and wood residue. Rice residues (straw and husk) have the largest potential among agricultural residues which account for 30% of total solid biomass potential in term of energy. Rice is the most important crop with annual production in 2016 of 44 Mt [START_REF]Statistical Handbook of Vietnam 2016[END_REF]. For 1 t of rice produced, there will be 1 t rice straw and 0.2 t rice husk generated. With heat value of straw about 10.9-12.6 MJ/kg and 10.8-12.1 MJ/kg for husk [START_REF] Arvo | Development of biomass fuel chains in Vietnam[END_REF], the total potential of rice residue can be estimated at 480-554 PJ per year and 95-106 PJ for straw and husk respectively.

In this study, we explore co-firing technology as a potential pathway to reduce carbon emissions from the fast-growing Vietnamese coal sector. We first assess the economic feasibility of implementing co-firing in Vietnam at national scale and then the potential emission reductions of co-firing technology deployment. To do so, we adapt the BeWhere model [START_REF][END_REF] in a recursive fashion, to determine a technically optimal scenario for Vietnam based on the biomass availability in the country and different development trajectories of the electricity sector. The results will get insight into the following aspects concerning co-firing practices: which coal plants should do co-firing, how much biomass will be used, what type of biomass, and when should it be implemented. Then, we evaluate the influence of different policies that affect the deployment of co-firing and the associated greenhouse gas emission reductions.

Method: the BeWhere model

Spatial explicit optimization is a method commonly used to optimize the allocation of resources. This method has been widely used for various types of resources such as solar, wind and especially for optimization of biomass supply chains. A recent study [START_REF] Nurariffudin | Spatial optimisation of oil palm biomass co-firing for emissions reduction in coal-fired power plant[END_REF] used spatial optimization method to assess the emission reduction from co-firing oil palm biomass with coal in Malaysia. In our study, we select BeWhere -a techno-economic engineering model for renewable energy systems optimization. BeWhere model integrates Geographical Information System (GIS) and Mixed-Integer Linear Programming (MILP) to identify the localization, size and technology of the renewable energy system that should be applied in in a specific region.

BeWhere has been used to optimize bioenergy system in regional scale and country scale [START_REF] Arvo | Development of biomass fuel chains in Vietnam[END_REF][START_REF]Statistical Handbook of Vietnam 2016[END_REF][START_REF][END_REF][START_REF] Nurariffudin | Spatial optimisation of oil palm biomass co-firing for emissions reduction in coal-fired power plant[END_REF]. This is the first study using BeWhere for co-firing in Vietnam. We develop the model further as a yearly recursive dynamic model. Vietnam has a national plan to develop coal power stations describing their locations, capacity and the year when each of them will be operated. Accordingly, we consider the location and size of existing and future coal power plants as exogenous, and only look at the types of cofiring technology that could be deployed, based on their investment costs, resource cost and hypothetical emission reduction benefits. We adopt a recursive dynamic approach and we analyze the deployment of co-firing technologies until 2030. These following assumptions have been made:  Coal power plants are built each year until 2030 as stated in the approved PDP 7A. We assume that the Plan with be implemented as scheduled. Therefore, our aim is not to identify/optimize the location and size of the new plants to be built but to see which plant uses co-firing, and when.  We assume the absence of any regulation that mandates coal power plants to generate part of their electricity from renewable energy such as the Renewable Standard Portfolio. Power plants' owner decide based on the profitability of co-firing option given an hypothetical internalized carbon cost, instead of an overall emission reduction target set by the government and allocated to each individual plant. Thus, the model represents decentralized decision about co-firing adoption.  We assume that the input costs and output prices will remain unchanged in the future: long term supply and sales contracts with no risks or revisions. It means once a plant has installed co-firing, its technological and economic conditions are fixed for the rest of its lifetime. Thus, recursive approach is more suitable since it is about taking investment decision as input costs and output prices will not be changed in the future [START_REF] Babiker | Forward-looking versus recursive-dynamic modeling in climate policy analysis: A comparison[END_REF].

In our recursive-dynamic approach, optimization is run every year for a set of candidate coal power plants for co-firing. Set of plants for the starting year includes all the existing plants.

Results indicate which plant will do co-firing with which technology and at which mixing ratio.

In the year Y, the set of candidate coal power plants includes the new plants scheduled to open according to the Power Development Plan and the plants in year Y-1 that is not switched to cofiring (Equation 1). Biomass resource is distributed among the candidate plants by the optimization. The amount of biomass available in year Y is determined by Equation 2.

Equation 1

CandidatePlants Y = ProjectedPlants Y + ( CandidatePlants Y -1 -CofirePlants Y -1 ) Equation 2 BiomassAvailable Y =BiomassGrowth Y -BiomassConsumed Y -1
In this study, we integrate a value of CO2 emissions in the optimization objective to represent policy incitation to co-firing in coal power sector. GAMS's objective (Function 1) is to minimize the production cost of electricity of the system, as in Equation 3, with consideration of avoided cost of emission reduction. Emissions include emission from electricity generation from coal combustion, emission from electricity from co-firing with biomass and emission from biomass transportation (Equation 4). Emission= Emission Coal + Emission Cofire + Emission Biomasstransport

Data and parameters

Biomass resource

In this study, we choose two types of raw materials: rice straw and rice husk as these rice residues are the most abundant types of biomass resources in Vietnam [START_REF] Arvo | Development of biomass fuel chains in Vietnam[END_REF] and most of them, especially rice straw, is not yet utilized for energy purposes but discarded as waste and usually burned in the fields [START_REF] Nurariffudin | Spatial optimisation of oil palm biomass co-firing for emissions reduction in coal-fired power plant[END_REF]. Spatial data on rice production in Vietnam in shapefile format is taken from the study of [START_REF] Laborte | RiceAtlas, a spatial database of global rice calendars and production[END_REF].

From this data, information on rice straw and rice husk availability is derived using Residue to Product Ratio (RPR). In this study, RPR of straw is 1 : 1, which means one ton of rice harvested will produce one ton of straw. RPR of husk is 1 : 0.2 [START_REF] Nguyen | Identification of biomass market opportunities in Vietnam[END_REF]. Biomass availability data is disaggregated into a grid of 0.2 x 0.2 degree resolution that covered all Vietnam and distributed to grid points where there is cropland for rice (Figure 1). Spatial data on administrative areas and cropland for Vietnam is obtained from [START_REF]Free Spatial Data | DIVA-GIS[END_REF].

Coal power plants

The study considers all existing and planned coal-fired utilities listed in the revised National Power Development Plan. Data on coal power plants was extracted from the revised PDP7 dated 18 March 2016 and from GreenID, Vietnam Sustainable Energy Alliance. These data include name, year of operation, capacity and location. Existing coal power plants' efficiencies are taken from [START_REF]Development of greenhouse gas control measures in coal-fired thermal power and proposed schedule to apply these measures[END_REF]. Most of existing coal power plants in used in Vietnam are deploying subcritical Pulverized Coal technology with efficiency range from 28% to 38% [START_REF] Laborte | RiceAtlas, a spatial database of global rice calendars and production[END_REF]. For the planned coal power plants, we assume the average efficiency at 38% and capacity factor of 0.7. Electricity generation of existing plants is taken from the coal database of GreenID, Vietnam Sustainable Energy Alliance [START_REF]Coal Power Plant Database[END_REF]. Location and capacity of coal power utilities are shown in

Co-firing technologies

We consider four existing co-firing technologies with different biomass mixing ratio and cost. These technologies are: (i) direct co-firing using a single boiler, blending biomass together with coal in coal miller; (ii) direct co-firing, milling biomass separately then inject directly to the burner; (iii) indirect co-firing with a gasifier, which converts biomass into synthetic gas (or syngas) that feeds the burner; and (iv) parallel co-firing in which a separate boiler is used for biomass and the generated steam is then mixed with steam from coal boiler [START_REF] Irena | Biomass Co-firing-Technology Brief[END_REF]. These technologies are both demonstrated in commercial scale or pilot scale. Direct co-firing is the simplest and cheapest option since it does not require major change in coal boiler to adapt to biomass penetration. However, the mixing ratio, percentage of biomass that replace coal on a heat basis, for this technology is 3-5% [START_REF] Agbor | A review of biomass co-firing in North America[END_REF][START_REF] Irena | Biomass Co-firing-Technology Brief[END_REF]. With separate injection option, more investment cost is required for new biomass miller, conveyor and injector. Indirect and parallel co-firing require separate gasifier of boiler for biomass, thus increase the investment cost significantly. However, the co-firing rate is much higher compare to direct co-firing technology, up to 20% according to [START_REF] Irena | Biomass Co-firing-Technology Brief[END_REF].

Technical assumptions

From rice data, we estimate the theoretical potential of rice straw and rice husk in Vietnam at 470 PJ per year and 95 PJ per year, respectively. This estimation is based on the annual rice production of Vietnam, the Residue to Product Ratio (RPR) and heat value of rice straw and rice husk. RPR of straw is 1:1 and of husk is 0.2:1. Heat value of straw is 11.7 MJ/kg [START_REF] Truong | Feasibility and sustainability of co-firing biomass in coal power plants in Vietnam[END_REF] and of husk is 11.9 MJ/kg [START_REF] Arvo | Development of biomass fuel chains in Vietnam[END_REF].

The technical potential, however, is lower than theoretical. This is because not all residues produced are collected and there are other demands for straw and husk rather than co-firing such as mushroom cultivation or cattle fed. In our study we assume fifty percent of straw and husk produced is available for co-firing.

We also assume that once a coal power plant retrofit for co-firing, they will ensure a long-term biomass supply every year until 2030. Thus, biomass availability for a certain year will exclude the part of biomass that used for co-firing in the previous years.

For our model input, assumptions on investment cost and the upper bound of biomass mixing ratio for each technology can be found in Table 1. Baseline emission factor for coal power plants in taken at 1.056 kg/MWh [START_REF] Finenko | Future Carbon Dioxide Emissions from Vietnam's Coal Power Generation Sector[END_REF].

Economic assumptions

We estimate the cost of producing electricity with pure-coal technology (Equation 5) and cofiring technologies (Equation 6) base on three elements: investment cost, operation and maintenance cost and fuel cost. The investment cost is annualized assuming a discount rate of 7% and plants' life time of 25 year. The assumptions used to calculate cost of electricity are listed in Table 1. [START_REF] Lüschen | Economic viability of biomass cofiring in new hard-coal power plants in Germany[END_REF][START_REF] Bxd | Decision[END_REF] Coal price ($/GJ) 3.72 [START_REF]Transparency in Electricty and Petrol business[END_REF] Straw price ($/GJ) 3.18 [START_REF] Truong | Feasibility and sustainability of co-firing biomass in coal power plants in Vietnam[END_REF] Husk price ($/GJ) 1.67 [START_REF] Arvo | Development of biomass fuel chains in Vietnam[END_REF] Fix O&M coal ($/kW/y) [START_REF] Babiker | Forward-looking versus recursive-dynamic modeling in climate policy analysis: A comparison[END_REF].31 [START_REF] Broadman | Logistics, Costs, and GHG Impacts of Utility-Scale Cofiring with 20% Biomass[END_REF] Variable O&M coal ($/ kWh) 0.0048 [START_REF] Broadman | Logistics, Costs, and GHG Impacts of Utility-Scale Cofiring with 20% Biomass[END_REF] Fix OM co-fire ($/kW/ y) 32.24 [START_REF] Broadman | Logistics, Costs, and GHG Impacts of Utility-Scale Cofiring with 20% Biomass[END_REF] Variable OM co-fire ($/kWh) 0.006 [START_REF] Broadman | Logistics, Costs, and GHG Impacts of Utility-Scale Cofiring with 20% Biomass[END_REF] Truck Transportation cost (USD/t/km) 0.088 [START_REF] Tai | Transportation price for truck[END_REF]

Sensitivity Analysis -uncertainty of fuel costs

Thermal power plants economics are sensitive to fuel costs. To see how fuel cost influences the model's results, we run local sensitivity analysis for eight different scenarios when cost of coal, straw and husk vary by ±10% of the reference value (see Table 2).

We run the model for each cost scenario with different social carbon value ranging from 0 to 30 USD/tCO2eq with a 1-USD-step to see how the emission reduction curve change under different fuel costs scenarios.

Results

Carbon value and emission reduction

We simulated scenarios with social carbon value (as used in Function 1) ranging from 1 to 30$/tCO2eq with a 1-USD step. Figure 3 plots the total emission reduction from all coal power plants by 2030 from these 30 social carbon value scenarios. As expected, results show that emission reduction from co-firing increases with social carbon value. The result is a staircase shaped curve, as the model essentially solves from discrete set of investment options: there are four technologies and two types of biomass.

Results show a 4 MtCO2eq at a zero social carbon value. Co-firing appears to be interesting even in the absence of emission reduction incentives.

Under our reference costs and technological performance assumptions, co-firing technology demonstrates its potential of reducing carbon emission from coal power plants in Vietnam.

When social carbon value is less than 8 $/tCO2eq, the emission reductions obtained by 2030 is only 4 MtCO2eq, about 1.1% compared to baseline emissions as can be seen in Figure 3.

When carbon value is 8 $/tCO2eq, emission reductions reach 27.3 MtCO2eq (7.6% of baseline emission). This jump is explained by the selection of raw material by the system, to be discussed in 4.2 below. From 8 $/tCO2eq, the emission reductions plot remains flat. This indicates that increasing carbon value above 8 $/tCO2eq does not have much impact to emission reduction. However, higher carbon price still has impact to the selection of co-firing technology by the system as described in 4.3 When social carbon value is not accounted (Figure 4 left chart), coal power plants will first select rice husk as co-firing feedstock as its market price is two times cheaper than that of rice straw as shown in Table 1. In the first year, 43 PJ of total 47 PJ technical potential husk is used to generate 14.5 PJ of electricity. The total share of biomass in power generation from coal power plants in Vietnam by 2030 in this case is only 1.2%.

Biomass utilization

From 2017 onward, the amount of biomass used for co-firing remains unchanged, meaning no more co-firing capacity is added. This is because without any incentive, it is not profitable for the plants to co-firing using rice straw, a feedstock that is twice more expensive than husk.

Co-firing with straw becomes competitive and used by the power plants when carbon value is at 8 $/tCO2eq (Figure 4, right chart). At this level of subsidy, co-firing can contribute to a 7.6% share to power generation by 2030 using 256 PJ of biomass. With higher availability than husk, rice straw can provide biomass for more plant to deployed co-firing in later years. The amount of straw co-fired by the system grows over years since new plants can only use straw for co-firing when husk is used up. However, rice residue is not sufficient for all power plants to do co-firing. This can be seen from Figure 5. When carbon value is 8 $/tCO2eq, biomass is used up by 2020 resulted in a flat line from this year onward. When carbon value is as high as 20 $/tCO2eq, rice residue run out even earlier in 2019 as more biomass is mobilized in the first year when higher biomass mixing ratio technology is deployed (as described in section 4.3).

Technology deployment

Two out of four co-firing technologies are selected: direct co-firing with separate biomass injection and indirect co-firing Figure 6. Although direct co-firing with fuel blend has lower investment rate than direct co-firing with separate injection, the coal substitution percentage is much lower. Thus, the power generation cost of fuel blend technology is slightly higher than that of separate injection technology. As shown in the left chart of Figure 3, the increase in carbon value higher than 8 $/tCO2eq does not change the total emission reduction of the system. However, the configuration of co-firing technology of the system does change with the increase of social carbon value.

When there is no carbon value, separate injection technology is selected to co-fire rice husk.

There is not enough rice husk produced in Vietnam for co-firing in all existing coal power plants, only 14 out of 28 existing plant can do co-firing in the first year and no more co-firing capacity added in the following years (Figure 6a). This is also the same situation for the system when carbon value is lower than 8 $/tCO2eq.

When carbon value is at 8 $/tCO2eq, all coal power utilities deploy direct co-firing technology with separate injection until year 2020 (Figure 6b), at which point no more plant is able to cofire with straw and husk because of resource depletion and the coal plants added after 2020 still run on coal only. After this year, the system keeps the configuration of co-firing with direct separate injection technology with the same biomass mixing ratio when carbon value varies from 8 to 15 $/tCO2eq.

Since indirect co-firing is a more expensive technology than direct co-firing, it is only deployed when carbon value exceed the thread hold of 16 $/tCO2eq (Figure 6c). However, only one plant uses indirect co-firing in this carbon value scenario. The number of plant using indirect co-firing technology increases with carbon value until the carbon value reaches 20 $/tCO2eq (Figure 6d).

Because the technology allows plants to co-fire at 20%, more straw and husk are utilized at the beginning. This results in rice residue resource depletion in 2019, a year earlier than in the scenarios where direct co-firing with separation injection technology is mostly used.

(a) (

Sensitivity analysis

Results of sensitivity analysis are presented in four charts in Figure 7 with different scenarios plotted with baseline case to see how the changes in different fuel impact the deployment of cofiring technology. For example, Figure 7d shows emission reduction curve (similar curve as in Figure 3) of baseline scenario in comparison to the emission reduction curves in S1 and S2 scenario where husk price increased and decreased by 10% respectively. Figure 7c shows the emission reduction curve when straw price varies. Figure 7b shows emission reduction curves for coal price sensitivity analysis and Figure 7a compares baseline scenario with "best case" scenario for co-firing, where coal price increased by 10% and straw and husk price both decreased by 10%, and "worst case" scenario where coal price decreased by 10% and biomass price increased by 10%. The results of sensitivity analysis (Figure 7) of fuel price (coal price, straw price and husk price) show that the co-firing system is sensitive to fuel price. However, the sensitivities vary from fuel to fuel.

Result shows that the system is least sensitive to husk price as the 3 lines ovelapped in Figure 7d. This is because husk price is about two times lower than straw and coal, thus the cost to purchase husk for co-firing is accounted for much smaller part in total fuel cost of coal power plants. Moreover, rice husk availability can only satisfy a small portion of electricity production (Figure 4, left chart). Thus, the impact of husk price change is much less important than changes in coal price or straw price.

The system has the same sensitivity to coal price and straw price. Straw price lowered by 10% (S4) and coal price increased by 10% (S1) can shift the emission reduction line to the left, which means it can help the system to reach the same level of emission reduction as in baseline case at lower carbon value. In these scenarios, carbon value to reach maximum emission reduction is about 4 $ lower than baseline case. On the other hand, straw price increased by 10% or coal price decreased by 10% both shift the emission reduction line to the right where higher carbon value is needed to reach the same emission reductions.

Discussion

We found that even without carbon value incentive, co-firing with rice husk is cheaper than using coal in coal power plants. The total cost of electricity generation of the system is reduced by 36 M$ compared to the baseline total cost of 20 389 M$ when co-firing husk. This is a 0.2% savings, which may not be enough to trigger industrial investment in a technology never used before in the country. The barrier is not only financial risk. Inadequate market structures and contracts issues alone can prevent a technology like co-firing to arise.

The limited rice husk resource can only afford co-firing in 14 plants out of total 96 coal power plants at a 15% mixing ratio. In this case, the amount of emission reduced is only 1.1% compared to baseline emission of coal power plants without co-firing. We interpret that without incentives, co-firing can be considered as a way to valorize rice residue, but is not likely to lead to large greenhouse gas emission reduction from coal power plants in Vietnam.

We found that with incentives, co-firing could have a significant role in reducing greenhouse gas emissions from coal-based power stations, up to 7.6% below baseline. This finding corroborates with the lessons from international experience: in all other countries, co-firing needed incentives to be deployed.

In our study, Vietnam coal sector can achieve 8% emission reduction already with a carbon value at 8 $/tCO2eq. Emission reductions remain at this level even when carbon value is as high as 30 $/tCO2eq. When comparing the price of carbon with the emission reduction it can deliver, we found that carbon value at 8 $/tCO2eq is the most interesting cost-effective. In this scenario, the total additional cost of the system for co-firing is 137 M$.

When carbon price is at 20 $/tCO2eq, most of the plants switch from direct to indirect co-firing. However, the total emission reduction is at the same level as when all co-firing plants deploy direct technology. The interpretation of that result is that indirect co-firing technologies are not competitive. Among the two co-firing technologies selected by the system, direct co-firing with separate injection offers a good compromise between mixing ratio and economic efficiency. This allows the technology to be deployed at moderate carbon value (8 $/tCO2eq). Indirect cofiring allows higher mixing ratio but required large amount of subsidy to be utilized. The additional cost for producing electricity of the system in case indirect co-firing is deployed is 206 M$. For the short and medium-term, direct co-firing is the better option for coal power plants in Vietnam when they want to generate part of their power production from renewable resources.

To complement the techno-economic analysis conducted so far, we offer a few remarks about opportunities and barriers for co-firing in Vietnam from macro-economic point of view.

First the recursive-dynamic model predicts that only rice residues are not enough for all coal power plants in Vietnam to co-fire at high biomass percentage. This is shown in Figure 4 for carbon value at 8 $/tCO2, in which from 2021 rice residue is not enough for all plants to co-fire at 15%. This limit can be relaxed if we consider international trade, or other biomass types for co-firing such as forest residue or biomass pellets. There is an active international market in biomass pellet for energy. Vietnam has production capacities and exported to South Korea and Japan, for example, with 2.5 million metric ton exported in 2017 [START_REF] Strauss | Vietnam Wood Pellet Exports Rapid growth and low prices. Can it last? What is the Future[END_REF].

Second, carbon value and biomass price interact to influence the choice of raw materials and cofiring technology. If raw materials are cheaper, then less incentive are needed for co-firing to be applied. Conversely, the emission reduction curve shifts to higher carbon price when biomass becomes more expensive. Markets are dynamic and uncertain, but there is a known pattern: when a biomass waste is valorized, it becomes a resource, and its price increases according to the law of offer and demand. For example in 2016 we observed that the demand for rice husk in Thailand market led to increases in price in the neighboring country of Cambodia [START_REF] Nguyen | A critical look at rice husk gasification in Cambodia: technology and sustainability[END_REF].

Third, sensitivity analysis shows that the adoption of co-firing is sensitive to coal price. Coal price in Vietnam for electricity generation is controlled by the government, who also regulates electricity retail price. Keeping both the price of input (coal) and output (electricity) low allows serving affordable electricity to the population, while controlling the debt level of the utilities. However, this low prices policy can be seen as a barrier to the adoption of renewable electricity technologies including co-firing.

Fourth, while our model uses a carbon value to internalize the climate policy constraints, incentives are not limited to price instruments. International experience shows that a renewable portfolio standard, which prescribes power generation companies to provide a fraction of green electricity in their mix, can lead to the adoption of co-firing.

Fifth, agricultural residues are a domestic resource and contribute to the energy security of supply. The latest Coal Development Plan approved by the Prime Minister predicts that coal demand for electricity generation will increase from 33 Mt in 2016 to 64 Mt in 2020 and to 131 Mt in 2030 [START_REF] Minister | Approval of revised Master plan for Coal sector development in Vietnam to 2020 with vision to[END_REF]. Supply side, the total coal production of the country is planned to grow to 57 Mt in 2030. We doubt that this will much be produced, since recent past statistics rather show that Vietnam passed peak coal already. In any case, many coal power plants will rely on imported coal. The 2030 plan is to import more than 75 Mt of coal, about 60% of coal demand for power plants. This will make the national economy more sensitive to international coal price variation. Co-firing would help to mitigate the vulnerability to international coal markets.

It should be noted that some issues might occur when blending biomass with coal in co-firing applications such as explosion risk from biomass storage, low grindability of biomass, fouling and corrosion of boilers [START_REF] Sahu | Coal-biomass co-combustion: An overview[END_REF][START_REF] Roni | Biomass cofiring technology with policies, challenges, and opportunities: A global review[END_REF][START_REF] Basu | Biomass co-firing options on the emission reduction and electricity generation costs in coal-fired power plants[END_REF]. Most of the technical issues result from the differences in combustion properties of coal and biomass, thus the performance of co-firing plants is largely affected by the type of feedstock adopted in the combustion process. For instance, rice residues are much more problematic to co-fire than conventional wood in term of grindability and fouling.

Co-firing thresholds that are considered in this study (5%, 15%, 20% and 20% coal replaced by biomass in term of heat for direct co-firing with fuel blend, direct co-firing with separate injection, indiect and parallel co-firing, respectively) should also be checked on a case-by-case basis, by performing detailed analysis of the fuel compositions, and/or pilot and full-scale testing.

Conclusion

We explored the potential of co-firing rice residues in coal power plants to reduce greenhouse emission reduction in Vietnam as one approach to tackle energy trilemma in Vietnam. The BeWhere model produced a marginal abatement cost curve which has two domains.

A/ At low carbon values, co-firing with rice husk blend in the grinder can profitably be used. This leads to a 1.1% sectoral emission reduction.

B/ At higher carbon values, rice straw is used in addition to husk. We estimate the economic potential to a 8% sectoral emissions reduction. Increasing the carbon value above further does not bring more co-firing capacity online, because the rice residue resource is limited.

The critical carbon value where rice straws is used is 8 $/tCO2eq in the central case. It depends on the prices of coal and straw. This carbon value may look high in a domestic policy context given that Vietnam is a lower middle-income country. However, it may be worth a look in the context of an international company seeking to reduce its global emissions.

Achieving ambitious national emission reduction policy targets -25% below baseline with international support by 2030-requires a portfolio of measures. In the power generation sector, besides demand management, measures include using other energy sources such as natural gas, wind and solar radiation to avoid building new coal power plants. But given the current trends, coal is to remain a major energy source in Vietnam in 2030. Mitigation policies will also include measures like increasing plants' thermal efficiency, biomass co-firing, and possibly negative emission technologies such as carbon capture and sequestration. In this portfolio of measures, co-firing has the potential to contribute a significant part to the emission reduction effort of Vietnam: a 8% reduction of the largest emission sector. With the difficulty that cofiring might fall under a global policy ban on anything coal-related, the technology is mature and appropriate to the development level of the country. Co-firing can utilize the existing infrastructure of coal power plants and domestic biomass resources. This is not a silver bullet, but an interestingly low-hanging fruit.

Figure 1 .

 1 Figure 1. Rice cultivation area (left), rice production by provinces (middle) and discretization of rice production into 0.2x0.2 degree grid (right).

Figure 2 .

 2 Figure 2. Existing coal power plants are mostly located in the Northeast of Vietnam where the largest coal reserves are located. In the future, many more coal power plants are planned in the South to satisfy the demand of the region. These plants are planned near the coast and ports to use imported coal. Recently built and planned coal power plants are switching to Circular Fluidized Bed technology with higher efficiency and more flexibility in terms of fuel used. Cofiring has been demonstrated to work with both coal-fired technologies. There are no plans to deploy the most expensive and efficient technologies such as ultra-supercritical or integrated gasification combined cycle (ICCG).

Figure 2 .

 2 Figure 2. Existing (left) and planned (right) coal power utilities in Vietnam

Figure 3

 3 Figure 3 Impact of carbon value to emission reduction from co-firing. The chart shows the emission reduction of allpower plants in 2030 when carbon value increase from 0 to 30 $/tCO2eq.

Figure 4

 4 Figure 4 shows the influences of carbon value on the penetration of biomass in coal power plants in two different carbon value scenarios.

Figure 4 .Figure 5 .

 45 Figure 4. Share of rice residue in power generation from coal power plants at different carbon value: 0 $/tCO2eq (left), 8 $/tCO2eq (right)

Figure 6 .

 6 Figure 6. Co-firing technology deployment at different carbon value.

Figure 7 .

 7 Figure 7. Sensitivity analysis of fuel costs impact to emission reduction. The three lines are overlapped in Figure 7d.

Table 1

 1 Generation Cost purecoal =C Inv +C OM + P Coal × Q Coal Summary of the technical and economic assumptions on the plants for reference scenario

	Equation 5
	Equation 6
	Generation Cost cofire =C Inv ' +C OM '	+P BM ×Q BM + P Coal ×Q Coal '
	Where
	C Inv is annualized investment cost for coal power plant
	C Inv '	is annualized investment cost of co-firing plant (including investment cost to build the
		coal power plant and additional cost to retrofit the plant for co-firing)
	C OM is operation and maintenance cost of coal power plant
	C OM '	is operation and maintenance cost of co-firing plant
	P Coal is price of coal
	P BM is price of biomass
	Q Coal is quantity of coal used in coal power plant
	Q Coal '	is quantity of coal used in co-firing plant and Q Coal -Q Coal '	is the amount of coal
		substituted by biomass

Q BM is quantity of biomass used in co-firing plant

Table 2

 2 Scenarios for sensitivity analysis

	Scenario	Baseline	S1	S2	S3	S4	S5	S6	S7	S8
		case
	Coal price	3.72	+10%	-10%					+10%	-10%
	Straw price 3.18			+10%	-10%			-10%	+10%
	Husk price	1.67					+10%	-10%	-10%	+10%

Acknowledgement

Part of the research was developed in the Young Scientists Summer Program at the International Institute for Applied Systems Analysis, Laxenburg (Austria) with financial support from the Vietnam National Member Organization. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.