Nicolas Dietrich

Kritchart Wongwailikhit

Mei Mei

Feishi Xu

Francisco Felis

Abderrahmane Kherbeche

Gilles Hébrard

Karine Loubiere

Karine Loubière

Using the "red bottle

Keywords: Upper-Division Undergraduate, Demonstrations, Chemical Engineering, Hands-On Learning/Manipulatives, Kinetics, Industrial Chemistry, Mechanisms of Reactions, Transport Properties, Organic Chemistry

come

INTRODUCTION

The "blue bottle" experiment, in which a transparent liquid turns blue when shaken and becomes colorless when left at rest, is a very popular reaction for both demonstration purposes and practical work in educational activities. It requires three chemicals: glucose, sodium hydroxide and methylene blue (dye) in water. The origins of this experiment are not well known. The first references to reactions between sugar and this dye (synthesized by Heinrich Caro in 1876) date back to the 1920s [START_REF] Lane | Estimation of Sugar in Urine by Means of Fehling's Solution with Methylene Blue as Internal Indicator[END_REF][START_REF] Dubois | Colorimetric Method for Determination of Sugars and Related Substances[END_REF][START_REF] Dutton | Methylene Blue -Reduction and Oxidation[END_REF] . The name "Blue Bottle Experiment" was first popularized by Campbell 4 in 1963 who demonstrated the occurrence of a 4-step reaction scheme between glucose, dissolved air (oxygen) and dye, and proposed a global equation of the reaction of glucose and air with the dye, regenerated in the reaction system, as a catalyst. Many variations and modifications of the pioneer experiment have been put forward. Researchers [START_REF] Cook | The Blue Bottle Experiment Revisited: How Blue? How Sweet?[END_REF] proposed the use of other sugars and dyes to achieve different reaction rates and colors or a "green" version [START_REF] Wellman | Greening the Blue Bottle[END_REF] based on the use of vitamin and avoiding the use of sodium hydroxide. Some of these systems can be implemented to display other colors, like the "red bottle" experiment, sometimes called the "Vanishing Valentine" [START_REF] Limpanuparb | The Colorful Chemical Bottle Experiment Kit: From School Laboratory To Public Demonstration[END_REF][START_REF] Rajchakit | Greening the Traffic Light: Air Oxidation of Vitamin C Catalyzed by Indicators[END_REF] , involving resazurin, or a sequence of three colors yellow/red/green called the "multi-colored bottle" using indigo carmine [START_REF] Hile | The Blue Bottle Revisited[END_REF] . These reactions have been used as an interesting activity on reaction mechanism identification [START_REF] Limpanuparb | Blue Bottle Experiment: Learning Chemistry without Knowing the Chemicals[END_REF] . Recently, alternatives compositions have been proposed [START_REF] Staiger | Variations on the "Blue-Bottle" Demonstration Using Food Items That Contain FD&C Blue #1[END_REF][START_REF] Rajchakit | Rapid Blue Bottle Experiment: Autoxidation of Benzoin Catalyzed by Redox Indicators[END_REF] by adding benzoin in the liquid phase in order to catalyze the reaction, and suggesting other molecules as the dye: safranin (red), phenanthrenequinone (orange), erioglaucine (light blue) and tetrazolium chloride (red). This reaction has been largely used as a demonstration but, recently, some authors have used it as an accurate tool [START_REF] Paul | Reaction Systems for Bubbly Flows[END_REF][START_REF] Rüttinger | Laser-Induced Fluorescence in Multiphase Systems[END_REF] to characterize the mass transfer at the gas-liquid interface in chemical engineering systems. The first works devoted to mass transfer around bubbles flowing in a millimetric channel have been reported by [START_REF] Dietrich | A New Direct Technique for Visualizing and Measuring Gas-Liquid Mass Transfer around[END_REF][START_REF] Yang | Visualization and Characterization of Gas-Liquid Mass Transfer around a Taylor Bubble Right after the Formation Stage in Microreactors[END_REF][START_REF] Yang | Local Investigations on the Gas-Liquid Mass Transfer around Taylor Bubbles Flowing in a Meandering Millimetric Square Channel[END_REF][START_REF] Krieger | Local Mass Transfer Phenomena and Chemical Selectivity of Gas-Liquid Reactions in Capillaries[END_REF][START_REF] Kováts | Optical Determination of Oxygen Mass Transfer in a Helically-Coiled Pipe Compared to a Straight Horizontal Tube[END_REF] . Later, the mass transfer was studied in a thin gap cell [START_REF] Kherbeche | Multi-Scale Analysis of the Influence of Physicochemical Parameters on the Hydrodynamic and Gas-Liquid Mass Transfer in Gas/Liquid/Solid Reactors[END_REF][START_REF] Xu | Fast Determination of Gas-Liquid Diffusion Coefficient by an Innovative Double Approach[END_REF][START_REF] Felis | Experiments on Reactive Mass Transfer around an Oxygen Bubble Rising Freely in a Confined Cell Using Colourimetric Methods[END_REF] and in a bubbly column [START_REF] Dietrich | Visualisation of Gas-Liquid Mass Transfer around a Rising Bubble in a Quiescent Liquid Using an Oxygen Sensitive Dye[END_REF][START_REF] Xu | Mass Transfer and Diffusion of a Single Bubble Rising in Polymer Solutions[END_REF] .

Thanks to a specific calibration procedure [START_REF] Dietrich | A New Direct Technique for Visualizing and Measuring Gas-Liquid Mass Transfer around[END_REF] and by considering a simple mass balance in the liquid phase [START_REF] Yang | Local Investigations on the Gas-Liquid Mass Transfer around Taylor Bubbles Flowing in a Meandering Millimetric Square Channel[END_REF] , the oxygen flux transferred and the gas-liquid mass transfer coefficients could be quantified in different types of systems. However, all these studies have required expensive highspeed or high-resolution image acquisition systems and specific, heavy image processing that are quite impossible to implement in a chemistry lab for students. Recently, a simple demonstration [START_REF] Piccione | Direct Visualization of Scale-Up Effects on the Mass Transfer Coefficient through the "Blue Bottle" Reaction[END_REF] has been proposed to illustrate the gas-liquid mass transfer in a stirred tank at different scales;

for that, they have used the "blue bottle" reaction to visualize the rate of mass transfer and how this rate was affected by the size of the vessel, the air flow rate and the rotational speed. The time required to obtain a steady color change could be analyzed and discussed in relation with the classical gas-liquid mass transfer correlations [START_REF] Roustan | Transferts gaz-liquide dans les procédés de traitement des eaux et des effluents gazeux[END_REF] in order to illustrate the effects of the operating conditions. These authors did not note that this time could be also directly linked to the overall gas-liquid mass transfer coefficient, and that a direct and simple calculation would be possible if the "red bottle" experiment was used instead of the "blue bottle" thanks to its faster kinetic properties [START_REF] Yang | Optical Methods to Investigate the Enhancement Factor of an Oxygen-Sensitive Colorimetric Reaction Using Microreactors[END_REF] .

The objective of the present paper is to propose a demonstration that enables the gas-liquid mass transfer to be visualized and, with the observation of the color change from colorless to red, the overall gas-liquid mass transfer coefficient to be determined. For illustrative purposes, this is applied in a small bubble column where bubbles are generated by sparging.

PRINCIPLE

This method is performed with the reaction system involving resazurin with sodium hydroxide and glucose in a small bubbly reactor (250 mL). In bubbly reactors, aeration provides a perfect mixing to the liquid phase (i.e. the concentration is uniform in the liquid phase), then, based on an oxygen mass balance in the liquid phase, one can write: (1) with the mass flux of oxygen transferred per unit of volume in g/s/m 3 , is the mass concentration of oxygen in g/m 3 , kL.a is the volumetric gas-liquid mass transfer coefficient in s -1 , E is the enhancement factor and the oxygen saturation concentration in g/m [START_REF] Dutton | Methylene Blue -Reduction and Oxidation[END_REF] .

In this study, a chemical reaction involving resazurin with sodium hydroxide and glucose in a small bubbly reactor (250 mL) has been considered. This reaction is fast (typically few milliseconds [START_REF] Yang | Optical Methods to Investigate the Enhancement Factor of an Oxygen-Sensitive Colorimetric Reaction Using Microreactors[END_REF] ): the concentration of oxygen, , remains always null ( ) as it is instantaneously consumed by the dye (dihydroresorufin) [START_REF] Dietrich | A New Direct Technique for Visualizing and Measuring Gas-Liquid Mass Transfer around[END_REF] . Some authors [START_REF] Dietrich | A New Direct Technique for Visualizing and Measuring Gas-Liquid Mass Transfer around[END_REF][START_REF] Yang | Local Investigations on the Gas-Liquid Mass Transfer around Taylor Bubbles Flowing in a Meandering Millimetric Square Channel[END_REF] have also demonstrated that the enhancement factor, E, is equal to 1 in the classical bubble column condition (namely, when liquid-side mass transfer coefficients kL are about 10 -4 m/s, see the supplementary section). As shown by Equation 2, two moles of dye react with one mole of oxygen:

O2 + 2 Dihydroresorufin (colorless) → 2 Resorufin (pink/red) + 2 H2O (2)
As a consequence, and because no oxygen is present in the liquid at the beginning, one can write:

(

Where , , are respectively the molar concentrations in oxygen, dihydroresorufin (colorless) and resorufin (colored) in mol/m 3 . By coupling Equation 1 & 3, the following equation is obtained:

By integrating Equation 4 from t = 0 to a time t, the following equation is deduced: (5) with M the molar mass (either of oxygen or dye) in g/mol and the dye concentration in g/m 3. When all the colorless reactant (dihydroresorufin) is consumed (complete conversion), the solution reaches its maximum color intensity, corresponding to a concentration where the red dye (resorufin) is equal to the initial concentration of the colorless dye (dihydroresorufin), namely . In such conditions, it becomes possible to calculate the volumetric gas-liquid mass transfer coefficient directly from the knowledge of the initial concentration of dye and the measured aeration time, as illustrated in Equation 5. For make the calculation easier, the initial concentration of the dye could fulfill the following condition:

Knowing that Mdye = 229.19 g/mol, MO2 = 32 g/mol and = 8.15 mg/L [START_REF] Dietrich | A New Direct Technique for Visualizing and Measuring Gas-Liquid Mass Transfer around[END_REF] , then this concentration corresponds to = 116.7 mg/L or = 0.51 mmol/L. In this specific condition, the calculation of the volumetric gas-liquid mass transfer is equal to the inverse of the coloration time as:

It can be noted that this concentration is only recommended for an easy and direct calculation of the mass transfer coefficient (Equation 7); the concentration of the dye can be indeed increased in order to make the experiment longer, thus increasing the precision of the measurement. The concentration can be also reduced in order to reduce the cost of the experiments.

It is interesting to observe that this experiment is very close to the sulfite static method classically used for the measurement of mass transfer coefficient [START_REF] Painmanakul | Effect of Different Contaminants on the α-Factor: Local Experimental Method and Modeling[END_REF][START_REF] Jamnongwong | Experimental Study of Oxygen Diffusion Coefficients in Clean Water Containing Salt, Glucose or Surfactant: Consequences on the Liquid-Side Mass Transfer Coefficients[END_REF] . In the latter method, for a given aeration time, the sodium sulfite (Na2SO3) reacts with the oxygen transferred into the liquid phase by the bubbles generated. The concentration of the remaining sulfite is determined by the iodometric method (oxidation of sulfite by iodine, and titration of excess iodine by sodium thiosulfate). As the initial oxygen concentration is kept at zero (nitrogen flushing), kL.a is then calculated from a very similar mass balance linking the aeration time to the amount of sulfite reacted:

As a consequence, this technique with the "red bottle" proposed in this paper can be seen as a visualization extension of the classical titration with sulfite.

It is also important to keep in mind that, at present, the technique can be applied only with resazurin. Indeed, the knowledge of additional data (kinetic constant and order, diffusion coefficients in basic solution of oxygen and dye) is needed to validate that the chemical regime is still fast with respect to the mass transfer process. This has been demonstrated for the resazurin system [START_REF] Yang | Local Investigations on the Gas-Liquid Mass Transfer around Taylor Bubbles Flowing in a Meandering Millimetric Square Channel[END_REF] but not yet for the "blue bottle" experiment and other systems. If this condition is not reached (ttransfer=1/kL.a >> treaction), the enhancement factor E can become higher than 1, meaning that the methods cannot be used as previously described and further calculations are needed.

MATERIAL

Experiments were carried out in a glass bubble column (0.045 m in diameter, 0.40 m high, 0.23 m height of liquid), submerged in a water bath regulated in temperature (Figure 1). The air sparger was a multiple orifice plate with two orifice sizes: 0.5 and 0.8 mm, having an interval of 1 cm between adjacent orifices. The gas was provided from a pressurized gas network in the laboratory. The gas superficial velocity (flowrate over the column section) was set to 0.6 cm/s, corresponding to a flow rate of 0.6 L/min. The column was filled with 250 mL of the colorimetric solution. The concentration of the dye was chosen at 116.7 mg/L in order to use the simple equation 7. In the presence of sodium hydroxide (20 g/L), and glucose (20 g/L), the resazurin (blue) was first reduced to resorufin (pink/red), which was then reduced to dihydroresorufin (colorless) as shown by Equation 2. The presence of glucose and sodium hydroxide in larger 165 amount than the dye made possible this reaction to be reversible (the resorufin being transformed back into dihydroresorufin), with a slow kinetics (few minutes). However, this reversibility remains limited by the initial amount of glucose, which tends to oxidize after several hours. For demonstrative purposes, the volumetric overall mass transfer coefficients (kL.a) associated with the train of bubbles generated were also measured using the dynamic method [START_REF] Wongwailikhit | Gas Sparger Orifice Sizes and Solid Particle Characteristics in a Bubble Column-Relative Effect on Hydrodynamics and Mass Transfer[END_REF] involving oxygen probes (Ox 50, Unisense, Denmark) and also using a high-speed monochromatic camera (Vision Research, Miro -M110, USA) as depicted in Figure 1. These elements are not required when proposing this demonstration to students; they are here used only to validate the technique.

DEMONSTRATION

The procedure consists in filling the vessel with the "red bottle" solution of potassium hydroxide, glucose, and resazurin (at the concentrations reported previously) and by sparging it with nitrogen gas until the pink/red color disappears. Then, nitrogen is replaced by air, at the desired flow rate, making the pink/red color reappears. The time where the color level remains unchanged can be measured with a stopwatch.

A typical sequence is represented in Figure 2. As the reaction is reversible, if the air is stopped, the solution becomes colorless again, and the experiment can be run another time (this operation could also be performed with a nitrogen flush). Thus, the student can experiment with many operating conditions in a short time. The authors recommend nevertheless operating the experiment less than 3 hours after making up the dye solution in order to minimize the glucose oxidation reaction. Sequential runs under different conditions can demonstrate the effect of orifice size on the mass transfer, and the effect of adding a packing system in the column (Figure

HAZARDS

Sodium hydroxide and its solutions are corrosive. Protective gloves and a laboratory coat should be worn together with safety eyewear. Dissolution of solid sodium hydroxide in water is a highly exothermic reaction and a large amount of heat liberated can pose a threat to safety through the possibility of splashing. Personal protective equipment such as dust mask, eyeshields and gloves should be used for the manipulation of the resazurin dyes.

RESULTS AND DISCUSSION

In order to validate the technique, the images shown in Figure 2 were analyzed by implementing an image processing software (Matlab, Mathworks, USA). The averaged grey level of the liquid part is plotted versus time in Figure 3.a. Initially, no color is observed (the grey level is maximal), because the oxygen concentration in the vessel is null. At t=0, air begins to be injected into the vessel and, as the reaction starts, the solution becomes pink/red and darker, corresponding to a decrease of the grey level. After more than a minute, the solution has turned pink/red and the grey level becomes constant as all the dye has reacted. It can be deduced that an aeration time of 70 s is required to reach 95% of the grey level value decrease. With this criterion, the precision of the measurement was estimated at  5 s. It can be observed in Figure 3.a that curve is not linear contrary to what is expected from equation 5. This could be explained by the fact that the bubbly column does not exactly behave as a perfectly mixed reactor, especially the parts close to the wall, as shown in Figure 2 The same analysis of the aeration time required to turn the liquid in the column red was performed in presence of a grid packing system located in the column. As depicted in Figure 3.b, the time obtained was 74 s, which is very close to the value found without packing.

These experiments were reproduced for a second orifice size (i.e., 0.8mm) at the same gas flowrate and aeration times of 85 s and 100 s were obtained with and without packing. The associated volumetric mass transfer coefficients were calculated with equation 7, and the values were reported in Table 1. In Table 1, the results were also compared with the ones deduced from the classical dynamic method (with oxygen probes). One can observe that a very good agreement (3% of difference) is obtained between two techniques, thus validating the visualization approach proposed. In addition, it can be noted that, for the small orifice, the packing does not have a notable impact on the mass transfer coefficient (2% of difference) but, for the bigger orifice, the packing increases the mass transfer coefficient significantly (15% of difference). This could be explained by the occurrence of bubble breakup on the packing grid, leading to a decrease in bubble size and an increase in interfacial area a (the ratio between the total bubble surface and the liquid volume, in m -1 ). Finally, the increase of the orifice size decreases the mass transfer coefficient strongly (25% of difference), which can be explained by the formation of bigger bubbles. These types of experiments, in different operating conditions, represent an interesting base to discuss with students on what parameters control the gas-liquid mass transfer phenomena and how the latter can be improved.

DEMONSTRATION WITH STUDENTS

The same set of experiments has been displayed to each member of a volunteer panel of 35 students and professors at the University of Toulouse. These sequences corresponded to the experiment performed in the device at two different orifice sizes, with and without packing. The panel was asked the simple question: "What is the time corresponding to the final change of color of the liquid?" The results obtained were averaged and the mass transfer coefficient corresponding to the panel as a whole has been added to Table 1. The difference between the three techniques was very small (less than 3%), but it can be observed that the panel always underestimated the mass transfer, and so overestimated the aeration time needed to reach the steady red color. This phenomenon could be explained by the persistence of vision occurring when an observer watched the bubble column for a long time. However, the overestimation was small enough to allow the experiment to be run with only the eyes of the students, without sophisticated and expensive cameras or probes. Finally, this activity is a striking example how a well-known chemical reaction could help in the quantification of the chemical engineering operation of gas-liquid mass transfer.

CONCLUSION

In this paper, an experimental setup and specific procedure have been proposed to adapt the well-known "red bottle" reaction so as to visualize and quantify chemical engineering principles related to mass transfer phenomena. The main advantages are the low cost of the materials and the ease of use. Contrary to the classical techniques, it requires no oxygen probes, no titration and no calibration procedure. The experiment is as simple as a colorimetric titration method. The new calculation method proposed, based on the measurement of the aeration time, requires only a stopwatch and, thanks to the reversible properties of the reaction, multiple runs of the visualization can be performed in succession. This new demonstration offers thus a striking, colorful visualization of the effects of various operating parameters (ie. orifice size, gas flow rate and packing) on mass transfer, but a number of other conditions can easily be compared, such as the effect of vessel size or of the sparger design. The demonstration can also be combined with a simple acquisition video with the objective to dissociate the effect of the interfacial area, a (linked to the bubble size and gas hold-up) and the mass transfer coefficient, kL. This experiment also constitutes a good introduction to discuss the competition between mass transfer and reaction. A further opportunity, not developed in this work, would be to perform comparisons of the recent adaptations of the blue bottle reaction system to "environment-friendly systems" 6 , reducing the hazardous conditions or systems catalyzed [START_REF] Rajchakit | Rapid Blue Bottle Experiment: Autoxidation of Benzoin Catalyzed by Redox Indicators[END_REF] .

SUPPORTING INFORMATION

Extensive discussion and detailed experimental methods (pdf file).

Figure 1 .

 1 Figure 1. Experimental setup.

Figure 2 .

 2 Figure 2. Evolution of the color in the column in function of the time. The superficial velocity of gas is 0.6 cm/s and the orifice diameter is 0.5 mm. The solution is getting more red in function of the time. These photographs are used to measure the gray light intensity over time.

  and Figure 3.b.

Figure 3 .

 3 Figure 3. (a) Evolution of the averaged grey level of the liquid in the bubble column as a function of time. The superficial velocity of gas is 0.6 cm/s and the orifice diameter is 0.5 mm. (b) A typical sequence for the experiments with a grid solid packing system. The superficial velocity of gas is 0.6 cm/s and the orifice diameter is 0.5 mm

Table 1 .

 1 Experimental values of the mass transfer coefficient.

	Operating condition	Orifice size	Solids	Panel kL.a (s -1 )	Image processing kL.a (s -1 )	Probe technique kL.a (s -1 )
			No solids	0.0134 0.003	0.0142 0.005	0.0135 0.001
	0.6 L/min	0.5 mm	Packing	0.0132 0.003	0.0136 0.005	0.0135 0.001
	(0.6 cm/s)		No solids	0.0975 0.003	0.0100 0.005	0.0102 0.001
		0.8 mm	Packing	0.0110 0.003	0.0117 0.005	0.0122 0.001

AUTHOR INFORMATION

Corresponding Author

Nicolas DIETRICH *E-mail: nicolas.dietrich@insa-toulouse.fr Personal website: ndietrich.com ORCID: orcid.org/0000-0001-6169-3101 Note: The author declares no competing financial interest.

The authors declare no competing financial interest.


Corentin Le Hesran
email: corentin.le-hesran@insa-lyon.fr

Anne-Laure Ladier

Valérie Laforest

Valérie Botta

Using flow

Keywords: Scheduling, Flow control, Waste minimization

Introduction

In the arising context of resource scarceness and environmental issues, new tools are being developed in order to promote greener manufacturing methods. Part of this endeavor involves flow assessment, which studies how resource flows (be they materials or energy) circulate within a production system and how they are consumed at the operational level. From a decision-maker's perspective, knowledge regarding the cost and environmental impacts of the various flows is important in order to consider trade-offs, especially since the real cost of waste flows tends to be severely underestimated (ADEME 2016). From a research perspective, solving shop-floor scheduling problems at the operational level has traditionally been done bearing only economic objectives in mind, such as production costs or profit [START_REF] Fang | A New Approach to Scheduling in Manufacturing for Power Consumption and Carbon Footprint Reduction[END_REF][START_REF] Giret | Sustainability in manufacturing operations scheduling: A state of the art review[END_REF]. A review of scheduling problems involving waste minimization is done in Le [START_REF] Hesran | Data gathering surveys[END_REF]. Providing environmental information linked to operational parameters would facilitate the integration of environmental aspects into the modeling of scheduling problems. Thus, this works focuses on the following research question: how can flow assessment support the identification and modeling of waste-minimizing shop-floor scheduling problems? To answer this question, a methodology is proposed and tested on an application case of hubcap manufacturing. Firstly, an overview of the current literature on flow assessment is done and the methodology steps are described. Secondly, a practical case is presented. Finally, conclusions are drawn in the last section.

Literature review

Several flow assessment methodologies were reviewed and their focus on environmental or economic assessment and decision level (strategic, tactical, operational) determined. The analysis of these methodologies, which come from various sources [START_REF] Jasch | The use of Environmental Management Accounting (EMA) for identifying environmental costs[END_REF][START_REF] Jasch | Environmental and Material Flow Cost Accounting: Eco-efficiency in industry and science[END_REF][START_REF] Cagno | Eco-efficiency for sustainable manufacturing: An extended environmental costing method[END_REF][START_REF] Gould | A Material Flow Modelling Tool for Resource Efficient Production Planning in Multi-product Manufacturing Systems[END_REF][START_REF] Vinodh | Life cycle assessment integrated value stream mapping framework to ensure sustainable manufacturing: a case study[END_REF][START_REF] Schubert | Implementation of energy-related aspects into model-based design of processes and process chains[END_REF][START_REF] Smith | Steps towards sustainable manufacturing through modelling material, energy and waste flows[END_REF] and are grouped in Figure 1, highlights the lack of environmentally oriented flow assessment tools for production scheduling.

Only five consider at the same time economic and environmental criteria as well as the operational decision level -and while those five studies do consider operational parameters in their approach, scheduling is not explicitly considered as an improvement lever. Only [START_REF] Despeisse | Sustainable manufacturing tactics and cross-functional factory modelling[END_REF] (Integrated FM) include production schedules in their model, and only [START_REF] Gould | A Material Flow Modelling Tool for Resource Efficient Production Planning in Multi-product Manufacturing Systems[END_REF] (MFAM) use production scheduling to improve the environmental performance, although economic performance is not considered in the results. However, although none of the reviewed methodologies are directly fit to answer our research question, they provide insights regarding the problem at hand which can be used to build our required framework. For instance, Life Cycle Assessment (LCA) has well-tried guidelines for defining the perimeter of a study and allocation methods for environmental impact assessment. Similarly, Material Flow Cost Accounting (MFCA) is broadly used for material flow inventory applications, both on environmental and economic aspects. From the operational point of view, the Input Throughput Output (ITO) method is effective in describing process parameters.

Thus, while flow assessment is a promising tool in order to promote waste-minimization through scheduling, no dedicated methodology for this specific purpose exists as yet.

Figure 1 -Methodologies grouping according to their included criteria and decision-level

Therefore, a new approach for facilitating the identification and modeling of waste-minimizing scheduling problems is proposed.

Methodology overview

Based on the results of the literature review and the insight provided by the existing studies, a four step methodology is devised. Rather than creating an entirely new framework, it combines relevant parts of current flow assessment tools in order to answer our research question. The four steps are represented in Figure 2. Step 1 serves to define the scope of the study, system and subsystems considered. Substep 1 is carried out similarly to the methodology proposed by the Environmental management -Life cycle assessment -Principles and framework (n.d.) for LCA, and the following items should be clearly identified and defined: boundaries and function of the system; functional unit; objective of the study; data requirements; assumptions and limitations regarding the study. Then, Substep 2 focuses on finding the most relevant subsystems -i.e. the groups of processes which have independent scheduling problems -to focus on. All the quantity centers (transformation, storage or transportation processes) are characterized using the ITO methodology [START_REF] Schubert | Implementation of energy-related aspects into model-based design of processes and process chains[END_REF]). Thus, their operating parameters related to scheduling, costs or waste are identified. Then, an estimation of the environmental impact is done for each subsystem using an LCA software. Adequate environmental indicators should be chosen depending on the waste type, surrounding ecosphere of the system and decision-maker preferences. The economic cost is determined using the Activity Based Environmental Costing (ABEC) method [START_REF] Cagno | Eco-efficiency for sustainable manufacturing: An extended environmental costing method[END_REF]. All situations where scheduling within a subsystem can influence waste generation should be listed, as these are the levers that will be used later on to improve the subsystem performance. The subsystems are then ranked by the decision makers, and the ones considered most important are further studied in the following steps.

In Step 2, a parametric flow inventory is carried out. Using the ISO 14051 (2011) framework, all the input, intermediary and output flows of the subsystems are defined. Using the operating parameters identified in the previous step, each flow can be quantified as a function of these parameters and the primary input flows entering the subsystem.

Proceeding in a downhill manner results in final output flows expressed according to the parameters of all the quantity centers they have crossed, or "parametric assessment".

Step 3 is the material flow assessment, where the waste output flows are characterized. Their environmental and economic evaluation is carried out, and their parametric representation is studied to identify how each parameter affects it. This serves the dual purpose of finding which flow or process is responsible for waste generation/cost and in which measure, as well as identifying all the parameters that can influence both scheduling and the quantity of waste.

Finally, in step 4 the scheduling problem is identified and modeled. From step 3, we know which scheduling parameters can be influenced to improve the system, i.e. what the decision variables for the scheduling problem are. Results from steps 1 and 2 regarding the system and subsystems yield information on the workshop configuration (α) and the constraints on production (β). Additionally, knowledge on the different waste flows and their associated costs gives us information on the possible environmental and economic improvements, i.e. the objective functions of the scheduling problem (γ). It then becomes possible to determine the problem's three field notation (α|β|γ) which is the common representation for scheduling problems [START_REF] Pinedo | Scheduling: Theory, Algorithms, and Systems[END_REF] as well as the decision variables that are used to solve it.

Case study

To illustrate the methodology, a practical application is carried out. This case involves a French hubcap production plant which includes raw plastic reception, oven drying, injection moulding, painting, quality control and expedition. The different steps of the methodology are successively applied in order to validate its feasibility.

Step 1: Study scope The production of the hubcap manufacturing plant ranges from raw materials reception and preparation to the expedition of finished products. In Substep 1, the product system consists in the whole production site, including all storage facilities for materials, products and waste. In order to gather information on the system, a survey was sent to the person responsible for production (Le Hesran, Corentin 2019), and a plant visit conducted. The production is composed of products grouped into three main families: plastic pieces, unicolor hubcaps and bicolor hubcaps. Additionally, stringent requirements placed on automotive parts suppliers places each batch of hubcap under a hard due date constraint. The functional unit (FU) chosen is the production of one day's worth of hubcaps, as the production schedule is determined on a daily basis. An average of 250 workdays per year was used in this study. Such a functional unit combines scheduling (through the daily planning of production) and waste generation (represented by the daily waste output in normal operating conditions). The daily production capacity is 25 000 hubcaps, with job sizes ranging from 800 to 2000 pieces, hence between 10 and 30 jobs per day. The spatial boundary considered for this product system is represented in Figure 3. Since due dates are involved, the temporal boundary for production is set as the last due date of the lots to be produced.

After having described the product system in its entirety, Substep 2 focuses on each independent subsystem to estimate their cost and impact, as well as the potential to mitigate these impacts using scheduling. As can be seen from the product system description in Figure 3, the plant is divided into three main workshops, namely the preparation, moulding, and painting/finishing. Buffer storage is present between each workshop, meaning that they can be considered as independent subsystems, as long as the buffer size and production capacity of each workshop are assumed to be sufficient.

The preparation workshop is responsible for producing the plastic used by the moulding machines. It generates few waste, namely packaging and wastewater. It has no constraints related to scheduling.

The moulding workshop manufactures the plastic pieces and raw hubcaps, and includes injection moulding machines, an assembly post as well as a quality control post. Its waste generation sources are residual plastic coming from the injection presses, and scrapped products from the quality control. From a scheduling perspective, waste production is impacted by changes in plastic compositions for the different pieces, as the machines need to be purged each time a different type of plastic is used. Once produced, the raw hubcaps and masks are sent to the painting and finishing workshop where they go through a painting line. Unicolor hubcaps only need a single coating, and go through the painting line only once before being sent to the finishing station. Bicolor hubcaps need to receive two coatings, with a mandatory 48h drying period between each coating in an intermediary storage. Painting masks are used during the second passage in the painting line and can be reused up to five times. All painted hubcaps are sent to the finishing line where a central logo is inserted and quality is controled. This workshop generates different types of wastes, namely paint sludge, scrapped products and used painting masks. Paint sludge is the result of soiled wastewater from the painting line going through an on-site flocculation process. It is considered a dangerous waste by the French environmental code (Code de l'environnement -livre V 2011) and needs to be stored in a separate building before collection for energy recovery. Paint sludge comes from two separate mechanisms: the normal functioning of the painting line, and the setup operations required when changing color. Similarly to what occurs in the moulding workshop, scheduling impacts the waste generation through the number of required setups, i.e. color changes. Thus, the moulding and the painting/finishing subsystems have been identified as opportunities for reducing waste through scheduling. In order to gather information on these subsystems and waste management, an interview was conducted with a Quality, Hygiene, Safety and Environment (QHSE) manager using another more detailed survey (Le Hesran, Corentin 2019). Missing information was extrapolated using studies from similar fields or from public sources. The yearly quantity of non-hazardous waste collected (not including scrapped products) is estimated to 54 tons. Non-hazardous waste is stored in outdoors metallic containers which were purchased by the company and do not have any renting cost. The price for waste collection and recycling was estimated at 180e per ton, based on price estimations by the French environmental agency (ADEME 2019). The price for one ton of PVC is estimated at 912e, based on recent French market prices (Ucaplast 2017), while the price for one ton of ready-to-use paint is estimated to 3000 euros. Operating prices were calculated based on the workforce of each workshop (The Boyd Company, Inc. 2019).

In the painting and finishing workshop, the company reported an average of 120 tons of paint sludge per year, with an annual cost of 38 000e for collection. This price includes neither the operation and maintenance cost of the flocculation plant nor the handling cost for packaging and transport into storage. [START_REF] Salihoglu | A review on paint sludge from automotive industries: Generation, characteristics and management[END_REF] report that costs for the flocculation station management represent around 46% of paint sludge management, which is the figure used for this study. A specific hangar is used for the paint sludge storage, further adding to the overall cost. Water is recirculated after treatment. Regarding environmental regulations, emissions levels of paint sludge are currently compliant. They are however a concern regarding the ISO 14001 certificate renewal.

Table 1 lists the information regarding the plastic and paint sludge waste flows. Treatment costs of paint sludge include management cost (152eper FU), the flocculation station operating cost (70eper FU) and waste storage cost (20eper FU). Environmental impacts were calculated using the OpenLCA 1.7.4 software and the Ecoinvent 3.1 database. LCIA method used is the eco-indicator 99, with three aggregated indicators (total ecosystems, total human health and total resources) for better clarity. As can be seen from Table 1, paint sludge has a larger environmental impact as well as a higher economic cost. Additionally, it is subject to governmental regulations, and a cause of concern regarding the ISO 14001 certification. For all these reasons, it was decided to limit this study to the painting and finishing workshop only. In order to characterize those quantity centers, the ITO method was applied to each of them. The indices corresponding to each parameter are shown in Table 2, and the detailed painting and finishing subsystem flow inventory shown in Figure 4. Paint sludge

z 1 = c r 3 × ((x 1 + x 2 ) × o w 1 + nbs 1 × s w 1 )
(1) Wastewater

z 2 = (1 -r r 3 ) × (1 -c r 3 ) × ((x 1 + x 2 ) × o w 1 + nbs 1 × s w 1 )
(2) Used masks

z 3 = x 2 × (1 -s r 1 ) × o w 6 (3) Scrapped products z 4 = x 1 × s r 1 + x 2 × s 2 r 1 (4) z 5 = x 1 × (1 -s r 1 ) × s r 4 + x 2 × (1 -s r 1 ) 2 × s r 4 (5)
Thus, the parametric representation allows us to quantify each flow circulating in the subsystem as an equation of the input flows and the operating parameters of the quantity centers. In addition, cost information is also represented, and will be used in the next flow assessment step.

Step 3: Material flow assessment After having quantified all the flows, their respective impacts and costs can be determined. More specifically, focus is given to the elementary output flows of waste, as those are the main factors to improve the environmental impact of the subsystem. Looking at the equations (1)-( 5) it can be seen that only the output flows z 1 and z 2 , respectively paint sludge and wastewater, are affected by the number of setups. Since no other parameter can be affected by scheduling, output flows z 3 , z 4 and z 5 are not considered in the rest of this analysis. While the number of setups does not appear in any other flow from the subsystem, it still affects the rest of production at the operational level in terms of lot-sizing. Increasing the number of setups tends to reduce lot-size, and conversely. This in turns affects the inventory (both for the intermediate and final storage) cost, as it depends on the number of products stored at any moment. Thus, the economic objective for this problem should include both the waste represented by flows z 1 and z 2 , as well as the inventory costs.

Based on the different activities and cost drivers described in the Ext-ABEC method (Enrico [START_REF] Cagno | Eco-efficiency for sustainable manufacturing: An extended environmental costing method[END_REF], the detailed cost equations of flows z 1 and z 2 are:

c z 1 = m c 1 × z 1 + c r 3 × o c 3 × y 1 + c r 3 × nbs 1 × set c 1 + s c 7 × y 2 + wt c 7 × y 2 (6) c z 2 = m c 2 × z 2 + (1 -c r 3 ) × o c 3 × y 1 + (1 -c r 3 ) × nbs 1 × set c 1 (7)
These costs are composed of different parts. In the case of z 1 , the meaning each term composing the equation is as follows:

• m c 1 × z 1 : material cost of z 1 , depends on the cost of flows x 3 , x 4 and x 5 ; • c r 3 × o c 3 × y 1 : cost of operating the flocculation station; • c r 3 × nbs 1 × set c 1 : setup cost for the painting line; • s c 7 × y 2 : storage cost for the paint sludge;

• wt c 7 × y 2 : waste collection and treatment cost.

Step 4: Scheduling problem identification and modeling Table 3 presents the process and outputs of the problem identification and modeling step.

β = d i , (a i , L, b i ), dependent setup-cost
Objective functions z envir : minimize waste from eq. ( 1) and ( 2) The α and β fields (workshop configuration and scheduling constraints) can be determined based on the information gathered during steps 1 and 2. The γ field (objective function) is identified during step 3 by considering all waste outputs and costs that can be influenced through scheduling. The data sets are determined using the information from step 1. Finally, decision variables are the production parameters that influence both scheduling and waste generation. They are first identified during substep 2, and then quantified during steps 2 and 3. After this step, the problem can be represented mathematically by translating the objective functions and constraints into mathematical equations using the defined data. This can be done using Mixed Integer Linear Programming (MILP). A MILP formulation of this problem can be found in Le [START_REF] Le Hesran | Reducing waste in manufacturing operations : a mixed integer linear program for bi-objective scheduling on a single-machine with coupledtasks[END_REF].

s w 1 × nbs 1 × ((c r 3 + (1 -r r 3 ) × (1 -c r 3 ));

Conclusion

This study presents an application case for a new methodology for the identification and modeling of waste-conscious scheduling problems using flow control. After a brief literature review, a four-steps methodology was proposed and applied to a case study of hubcap manufacturing. After defining the study scope, the product system was decomposed into independent subsystems. Environmental and economic impacts were estimated and the best subsystem to study chosen. Using the operational information gathered in step 1, a parametric flow inventory was conducted, providing a full description of material flows using the production parameters. An assessment of the waste flows was then conducted to identify possible improvements using scheduling. In the final step, a three field notation of the associated scheduling problem was provided and relevant data and decision variables identified. While more case studies need to be carried out and a full framework developed for this methodology, it has proven to be effective in identifying a scheduling problem with waste minimization concerns and given a basis for a full problem modeling.

Figure 2 -

 2 Figure 2 -Proposed methodology implementation steps

Figure 3 -

 3 Figure 3 -Hubcap product system description

Step 2 :

 2 Parametric flow inventoryThe quantity centers contained in the painting and finishing workshops are:

Figure 4 -

 4 Figure 4 -Painting and finishing subsystem flow inventory

Table 1 -

 1 Subsystems ranking: moulding and painting workshop wasteflows

		Impact	Scrap Plastic	Paint Sludge
		Material intensity	216 kg per FU	480 kg per FU
	Environmental	Ecosystems (PDF × m 2 × year) 0,208 per FU	4,56 per FU
		Human health (DALY)	0,27 per FU	6,62 per FU
		Resources (MJ surplus)	1,0 per FU	36,96 per FU
		Materials cost	197 euros per FU	1440 euros per FU
	Economic	Systemic cost	4901 euros per FU 3770 euros per FU
		Treatment cost	39 euros per FU	242 euros per FU

Table 2 -

 2 Parameters and flow indices Waste treatment cost nbs Number of setups set w Setup waste QC Quantity center set c Setup cost o w Operating waste

	Cost parameters	Operational parameters	Waste parameters		Flows
	m c Material cost	p r Production rate	s r Scrap rate	x Initial input flow
	o c Operating cost	cap Capacity	c r Conversion ratio	y	Intermediary flow
	s c Storage cost	s t Setup time	r r Recirculation ratio z	Final output flow
	wt c

Table 3 -

 3 Problem identification and modeling process

		Identification process	Resulting notation
		The painting line is the only relevant process to
	Workshop configuration	consider, the mask pose and finishing station have sufficient capacities and can be ignored in the	α = 1
		scheduling problem → single machine problem
		Due dates d i ; sequence-dependent setup cost;
	Constraints	Coupled tasks constraint (a i , L, b i ) (Blazewicz
		et al. 2012)

  z eco : minimize waste and inventory costs nbs 1 × set c 1 + inventory cost (intermediary, final) γ = min(z envir , z eco )

	Data sets	I: set of batches to be scheduled; J : set of operations composing a batch	I, J
		s Main decision
			variable: s ij ;
	Decision
	variables		Secondary decision
			variables: y

ij : starting time of operation j of batch i; y ijkl : 1 if operation j of batch i takes place just before operation l of batch k, 0 otherwise (operations order); t ij : drying time after operation j of batch i (intermediary inventory cost); e i : earliness of batch i, i.e. the time between the completion date and the due date of batch i (final inventory cost); nbs: nb of setups (environmental impact/cost) ijkl , t ij , e i , nbs

* This project is supported by the Auvergne Rhône-Alpes region.


M.Sc Kelly T Bocanegra-González

Mrs Kelly

Tatiana Bocanegra González

Kelly Tatiana

M.Sc Bocanegra González
email: ktbocanegrag@gmail.com

Ph.D Marie-Laure Guillemin
email: marielaure.guillemin@gmail.com

Guidelines for the restoration of the tropical timber tree Anacardium excelsum: first input from genetic

Keywords: National Strategy for Plant Conservation of Colombia, translocations, seasonally dry tropical 48 forest, molecular markers, ISSR, cpDNA

come

Guidelines for the restoration of the tropical timber tree

Anacardium excelsum: first input from genetic data

Introduction

51

During the last decades, increased efforts to restore habitats degraded by human activities throughout the 52 world have been made [START_REF] Meli | A global review of past land use, climate, and active vs. 526 passive restoration effects on forest recovery[END_REF]). In the case of tropical forests, most of these efforts have focused on 53 restoring species diversity and forest structure post-logging, while only a limited number have focused on 54 controlling the erosion of genetic diversity or help restore this diversity among species of ecological and 55 70 struggling. Translocation has been key in establishing new populations or re-establishing extirpated 71 populations and to increase numbers of individuals in critically small populations [START_REF] Menges | Restoration demography and genetics of plants: when is a translocation successful?[END_REF]Weeks et al. 72 2011). On the one hand, as in the case of the palm tree Pseudophoenix sargentii in Florida (Maschinski and

73

Duquesnel 2007), this technique has proven to be a very successful mean to counteract genetic erosion,

74

allowing to preserve the evolutionary potential of this endangered species. On the other hand, the 75 unsuccessful restoration experiments through translocations of seeds from Pinus pinaster in Western Europe 76 highlight the importance of acquisition of previous knowledge about the genetic structure of the species of 77 interest [START_REF] Benito-Garzón | Habitat restoration and climate change: 453 dealing with climate variability, incomplete data, and management decisions with tree translocations[END_REF]. Indeed, transplanted populations of P. pinaster have been unable to settle 78 in sites characterized by climatic conditions different from the ones encountered in the seedlings' area of 79 origin; a result that has been related to the existence of local adaptation across the species' natural range 80 [START_REF] Benito-Garzón | Habitat restoration and climate change: 453 dealing with climate variability, incomplete data, and management decisions with tree translocations[END_REF]). Unfortunately, it seems that the case of P. pinaster is not unique and that more 81 negative than positive experiences are reported in the process of restoration through translocations (Keller et 82 al. 2000;[START_REF] Menges | Adaptive introductions: How multiple experiments and 530 comparisons to wild populations provide insights into requirements for long-term introduction 531 success of an endangered shrub[END_REF]. It is possible that a failure at integrating knowledge not only on the species 4 taxonomy, reproductive biology, demography, and ecology but also on the distribution of genetic diversity 84 and functional adaptive traits between populations before choosing "whom to translocate where", account for 85 the majority of unsuccessful restoration programs [START_REF] Lesica | Ecological genetics and the restoration of plant communities: mix or match? 517[END_REF][START_REF] Menges | Restoration demography and genetics of plants: when is a translocation successful?[END_REF]; Aerts and 86 Honnay 2011).

87

Moreover, restoring ecosystem functions and biological diversity often requires the use of species specific to 88 that habitat [START_REF] Lesica | Ecological genetics and the restoration of plant communities: mix or match? 517[END_REF][START_REF] Mckay | How local is local?" -a review of practical and conceptual 522 issues in the genetics of restoration[END_REF][START_REF] Benito-Garzón | Habitat restoration and climate change: 453 dealing with climate variability, incomplete data, and management decisions with tree translocations[END_REF], and the need to use native 89 species in ecological restoration processes has been much emphasized lately. Although most processes of 90 reclamation of degraded land have been accomplished using exotic species [START_REF] Lamb | Large scale ecological restoration of degraded tropical forest lands: the potential role of 515 timber plantations[END_REF][START_REF] Meyerson | Exotic plant species as problems and solutions in ecological restoration: a 470 synthesis[END_REF] 91 [START_REF] Meyerson | Exotic plant species as problems and solutions in ecological restoration: a 470 synthesis[END_REF], in recent years, projects that include species native to the ecosystems to be recovered have 92 increased. For example, in Colombia, there have been recent initiatives to restore the seasonally dry tropical 93 forests (SDTF) using native species [START_REF] Thomas | Genetic diversity of Enterolobium cyclocarpum in 570 Colombian seasonally dry tropical forest: implications for conservation and restoration[END_REF][START_REF] Galindo-Rodríguez | Seed desiccation tolerance and dispersal in tropical dry forests 486 in Colombia: Implications for ecological restoration[END_REF].

94

Anacardium excelsum is a native tree of the tropics of America, and can be found in riparian forests of the dry 95 tropical and premontane areas from Guatemala to the northern part of South America, including the Guianas 96 (Bernal et al. 2015). This species is used extensively in forestry and, although it has been described as having 97 medium to slow growth rates, a recent study showed that it can reach the desired dimensions and 98 characteristics for logging in approximately 20 years [START_REF] Lozano | Estimación del crecimiento diamétrico, de Anacardium excelsum 519 (Kunth) Skeels, por medio de modelos no lineales, en bosques naturales del departamento del Tolima. 520[END_REF]. This feature, together with its 99 capacity to grow in a wide range of soils and climatic conditions, makes A. excelsum an ideal native candidate 100 for ecological restoration programs of riparian and lowland forests of Central and South America (Lamb 101 1998;[START_REF] Garen | The tree planting and protecting culture of cattle ranchers 488 and small-scale agriculturalists in rural Panama: opportunities for reforestation and land restoration[END_REF]. In spite of these interesting features, habitat loss and overexploitation are now 102 threatening A. excelsum populations, especially in Panamá and Colombia [START_REF] Condit | Demography and harvest potential of Latin American timber species: 466 data from a large, permanent plot in Panama[END_REF]

111

In the American tropics, floristic differences [START_REF] Banda-R K | Plant diversity patterns in neotropical dry forests and 442 their conservation implications[END_REF]) and the genetic structure of species of trees 112 commonly encountered in SDTF remnants [START_REF] Caetano | The history of Seasonally Dry Tropical Forests in eastern 459 South America: inferences from the genetic structure of the tree Astronium urundeuva 460 (Anacardiaceae)[END_REF][START_REF] Thomas | Genetic diversity of Enterolobium cyclocarpum in 570 Colombian seasonally dry tropical forest: implications for conservation and restoration[END_REF] have been shown to be 113 determined by historical (i.e. Pleistocene and Holocene) as well as more current changes affecting the 114 ecosystem. Indeed, the Pleistocene Arc theory proposes that species characteristic to SDTF ecosystems had 115 largely contiguous distributions during the late Pleistocene's cool and dry glacial period, and that contractions 116 in present-day SDTF remnants correspond to a more recent change due to climate warming during the 117 Holocene [START_REF] Prado | Patterns of species distributions in the dry seasonal forests of South America[END_REF][START_REF] Pennington | Neotropical seasonally dry forests and Quaternary vegetation 550 changes[END_REF]. Anacardium excelsum populations could thus have 118 been much larger and highly connected in the past. Moreover, from the 1500's onwards, anthropogenic 119 disturbance linked to human settlement after European colonization have severely degraded SDTF remnants 120 [START_REF] Valencia-Duarte | Dinámica de la vegetación en un enclave 576 semiárido del río Chicamocha, Colombia[END_REF]Pizano and García 2014). This could have led to further restriction in gene flow 121 and loss of genetic diversity through bottlenecks in A. excelsum.

122

In order to define the level of genetic structure between regions, and the level of genetic diversity within 123 regions and within remnants of A. excelsum in Colombia, we genotyped 106 trees collected from the 124 department of Bolívar in the north, down to the department of the Valle del Cauca in the south using nuclear 125

Inter-Simple Sequence Repeats (ISSRs), and also sequenced two non-coding chloroplast (cpDNA) loci. ISSR 126 and cpDNA loci represent molecular markers with highly distinct evolutionary rates. ISSRs are highly 127 polymorphic since their primers anneal directly on the SSR repeats distributed in the genome, which exhibit a 128 high rate of mutation classically linked to additions and losses of repeat units due to DNA slippage (Tuisima 129 et al. 2016). Conversely, the low rate of nucleotide substitutions, characteristic to the chloroplast genome in 130 plants, has been used to develop molecular markers well suited for phylogenetic and phylogeographic studies 131 [START_REF] Hamilton | Four primer pairs for the amplification of chloroplast intergenic regions with 495 intraspecific variation[END_REF]. Thus, we hypothesize that cpDNA loci will provide information about ancient gene flow (as 132 their sequences are rather conserved over time), while our ISSR results will reflect more recent patterns of 133 connectivity between A. excelsum remnants in Colombia. We further investigate population genetic footprints 134 of past climatic changes and current anthropogenic disturbances in a tree species characteristic to the

157

Total genomic DNA was extracted using the protocol described by [START_REF] Doyle | A rapid DNA isolation procedure for small quantities of fresh leaf tissue[END_REF].

158

Ten ISSRs primers were used for genotyping (Table S1). Amplification cycles consisted of a 5 minutes 159 initiation step at 94°C, followed by 40 amplification cycles consisting of a denaturation step at 94°C for 1 160 minute, annealing at 46°C-57°C for 45 seconds (Table S1) and a 72°C extension step for 2 minutes. Each

161

PCR was set up in a volume of 20µl using a PCR mix of: 3µl of DNA (corresponding to more or less 10ng); 162 2µl of PCR 10X buffer; 2μl of 2.5mM dNTP; 1.8μl of 25mM MgCl2; 0.9μl of 2.5mg/mL BSA; 0.4μl of 5U/μl

170

Two chloroplast intergenic spacers, trnL-trnF [START_REF] Taberlet | Universal primers for amplification of three non-coding regions of 568 chloroplast DNA[END_REF]) and rpl20-rps12 [START_REF] Hamilton | Four primer pairs for the amplification of chloroplast intergenic regions with 495 intraspecific variation[END_REF] were 171 sequenced. For trnL-trnF we used the primers B49317 5'-CGAAATCGGTAGACGCTACG-3' and A50272 172 5'-ATTTGAACTGGTGACACGAG-3' described in [START_REF] Taberlet | Universal primers for amplification of three non-coding regions of 568 chloroplast DNA[END_REF]. For rpl20-rps12, the primers 173 rpl20 5'-TTTGTTCTACGTCTCCGAGC-3' and rps12 5'-GTCGAGGAACATGTACTAGG-3' were used 174 [START_REF] Hamilton | Four primer pairs for the amplification of chloroplast intergenic regions with 495 intraspecific variation[END_REF]. Fragments were amplified by PCR in a 30µl volume using a mix similar to the one 175 described above for the ISSRs but with a final concentration of 0.2mM MgCl2, 0.4mM of each primer, and 176 0.2mM of each dNTP. A volume of 0.2μl of 5U/μl Taq DNA polymerase (Invitrogen, Carlsbad, USA) was 177 used for each PCR. For both cpDNA loci, the amplification was done following the PCR cycling parameters 178 proposed by [START_REF] Taberlet | Universal primers for amplification of three non-coding regions of 568 chloroplast DNA[END_REF]. Products were purified using MoBio UltraClean® DNA (Solana Beach, 179 USA) and sequenced using both PCR primers, at Macrogen (https://dna.macrogen.com). The distinct 180 haplotypes detected were deposited in GENBANK under accession numbers MG309720-MG309723 and 181 MG309724-MG309726 for rpl20-rps12 and trnL-trnF, respectively. 182 183

Statistical analysis

184

ISSR amplified fragments were scored manually as absent (0) or present (1) to build a binary matrix, and 185 amplification failure of a sample was treated as missing data (-1). All ISSRs were amplified and scored twice 186 in order to limit the number of missing data. Each band was assumed to represent the phenotype at a single 187 biallellic locus [START_REF] Crema | High genetic diversity detected in the endemic Primula 468 apennina Widmer (Primulaceae) using ISSR fingerprinting[END_REF]. Classical parameters of genetic diversity suitable for 188 dominant molecular markers, calculated using GenAlEx 6.5 [START_REF] Peakall | GenAlEx 6.5: genetic analysis in Excel. Population genetic software for 548 teaching and researchan update[END_REF], were: average number 189 of effective alleles (Ne), number of private bands (Pa), Shannon index (I), and expected heterozygosity (He). 190 GenAlEx 6.5 was further used to calculate Nei's genetic distance [START_REF] Nei | Estimation of average heterozygosity and genetic distance from a small number of individuals[END_REF] in order to conduct a Principal

191

Coordinates Analysis (PCoA), and to construct a dendrogram using Sequential Agglomerative, Hierarchical,

192

and Nested Clustering (SAHN, in NTSYS-pc;[START_REF] Rohlf | NTSYS-pc: Numerical taxonomy and multivariate analysis system version 1.80. Exeter 562 Software[END_REF]. A population genetic structure analysis was 193 carried out based on the model of Bayesian grouping implemented in STRUCTURE [START_REF] Pritchard | Inference of population structure using multilocus genotype 560 data[END_REF].

194

The LOCPRIOR option was implemented to detect weak signals of population structure (i.e. includes a priori 195 sampling locations as prior information in the clustering algorithm; [START_REF] Hubisz | Inferring weak population structure with the assistance of 499 sample group information[END_REF]). The admixture model 196 was run, assuming correlated allele frequencies [START_REF] Falush | Inference of population structure using multilocus genotype data: 484 linked loci and correlated allele frequencies[END_REF], using an initial burn-in of 500,000 197 followed by 1,000,000 Markov chain Monte Carlo (MCMC) iterations. Seven independent runs were 198 performed with the tested number of clusters (K) varying from 1 to 13 (i.e., the number of sites with ISSR 199 data + 1). Values obtained for the mean posterior probability of the data (Ln Pr(X|K)) and calculated for each 200 value of K [START_REF] Janes | The K = 2 conundrum[END_REF] were plotted in STRUCTURE HARVESTER

201

(http://taylor0.biology.ucla.edu/structureHarvester/). The clear maximum of the Ln Pr(X|K) values was used 202 to define the optimal K. CLUMPP [START_REF] Jakobsson | CLUMPP: a cluster matching and permutation program for dealing with 505 label switching and multimodality in analysis of population structure[END_REF] was used to combine the results of the 203 seven independent runs, applying the greedy algorithm with 100,000 random input orders. To estimate the 204 partitioning of genetic variability among regions and among individuals within regions [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances 479 among DNA haplotypes: application to human mitochondrial DNA restriction sites[END_REF],

205

an analysis of molecular variance (AMOVA) was performed with GenAlEx 6.5 by resampling 9,999 times.

206

Pairwise values of ST among regions were calculated with GenAlEx 6.5 and their significance was tested 207 using 9,999 permutations. For both the AMOVA and ST calculations, regions were defined as in Table 1 (i.e.

208

Magdalena river valley, Cauca river valley, Caribbean coast and Chicamocha canyon) and all individuals 209 from distinct sampling sites within each region were pooled. A Mantel test was performed with GENETIX 210 4.05 (Belkhir et al. 1996(Belkhir et al. -2004) ) to examine the relationship between Nei's genetic distance and geographic 211 distance (estimated as ln of the geographical distance; Rousset 1997), and its significance was tested using 212 10,000 permutations. ISSR data from BOL (i.e., the botanical garden) were not included when performing the 213 AMOVA among regions, the calculation of ST among regions, and the test for isolation by distance. 1). For both chloroplast loci, a median-joining network was constructed using NETWORK v 223 4.5 [START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF]. A rooted UPGMA tree was constructed in MEGA7 [START_REF] Kumar | MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for 513 bigger datasets[END_REF]) using the 224 concatenated trnL-trnF and rpl20-rps12 sequences with 1,000 bootstrap replicates. Only one of each distinct 225 sequences was used and gaps were coded by hand as a 5 th character. Anacardium occidentale (Genbank 226 accession number: KY635877) was used as an outgroup. To estimate the partitioning of genetic variability 227 among regions and among individuals within region [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances 479 among DNA haplotypes: application to human mitochondrial DNA restriction sites[END_REF], analyses of molecular variance

228

(AMOVA) were performed with ARLEQUIN v. 3.5 [START_REF] Excoffier | Arlequin ver. 3.0: an integrated software package for population 482 genetics data analysis[END_REF]) by resampling 10,000 times.

229

Pairwise values of ST among regions were calculated with ARLEQUIN v. 3.5 and their significance was 230 tested using 10,000 permutations. For both, the AMOVA and ST calculations, regions were defined as in 231

235

Of the ten ISSR primers used for genotyping, three showed a high level of missing data (>5%) and were 236 removed before analyses. Three individuals of the SUC population (Caribbean coast region; Table 1) also 237 showed a high level of missing data (i.e. more than half the ISSR primers did not amplify), and these

244

The seven ISSR primers selected generated 57 fragments. Total number of fragments scored per primer 245 ranged from 6 to 10 (Table S1) with an average of 69.9% polymorphic loci within locality (Table 1). Shannon 1). I and He 249 decreased then from the Magdalena river valley to the Cauca river valley, to the Caribbean coast and, at last, 250 to the Chicamocha canyon (Table 1). The lowest values were encountered in SAN (I = 0.23, He = 0.16; Table 251 10 1). Among the ISSR markers, only one private fragment was observed for the primer UBC-809 (Table S1) in

252

VEN. When estimated by region, He was 0.35 in the Magdalena river valley, 0.32 in the Cauca river valley, 253 0.25 in the Caribbean coast and 0.18 in the Chicamocha canyon.

254

Three haplotypes were observed for trnL-trnF with five polymorphic sites that included a gap of 4bp (Table 255 2). For rpl20-rps12, four haplotypes and three polymorphic sites were encountered (Table 2). Haplotypic 256 diversity within locality was generally very low and only one haplotype was observed for both loci in most 257 localities (i.e. ARM, IBG, CHA, CLI, QUI and SAN; Table 2). The highest values of diversity were observed 258 in the Caribbean coast (BOL: Hd = 1.00 and 0.67, = 4.57 10 -3 and 2.60 10 -3 for trnL-trnF and rpl20-rps12, 259 respectively; CES: Hd = 1.00, = 2.60 10 -3 for rpl20-rps12; Table 2). Intermediate values were observed in 260 MED and MAR for both loci and in VEN for rpl20-rps12 (Table 2). When estimated by region, haplotypic 261 richness (Re) was 2.0 and 2.0 in the Caribbean coast, 0.9 and 0.8 in the Cauca river valley, 0.1 and 0.

273

High genetic distances were observed between this group and COR and CES from the Caribbean coast, SAN 274 from the Chicamocha canyon and MED located at the far north of the Cauca river valley (Fig. 2). Of all these 275 four last localities, only SAN and MED clustered together. The Bayesian analysis carried out in

276

STRUCTURE for all localities shows the existence of three predominant genetic clusters (Fig. 3A). Most 277 individuals displayed varying degrees of mixed ancestry, and two of the three clusters present a rather 278 restricted geographic distribution (Fig. 3B andC). Individuals assigned, with a membership of 0.60 or higher,

279

to the red genetic cluster were only observed in the Magdalena river valley (Fig. 3C) while individuals 280 assigned to the yellow genetic cluster were observed in the Magdalena river valley, the northern part of the 281 Cauca river valley and Chicamocha canyon (Fig. 3C). Interestingly, all trees sampled in SAN were assigned 282 to the yellow genetic cluster and presented a very low level of admixture (Fig. 3B). Individuals assigned to 283 the blue genetic cluster with a membership of 0.60 or higher were observed in all the sampled localities 284 excepting SAN (Fig. 3C). Geographic isolation strongly influenced genetic differentiation among localities, as 285 demonstrated by the significant correlation observed between genetic distance and ln of geographic distance 286 (r = 0.749, p < 0.01) based on the Mantel test (Fig. S1). Significant genetic differences among regions were 287 detected by AMOVA, however, only 6% of the total of genetic variation was attributed to differences among 288 regions (Table 3). Significant genetic differentiation ( ST) was observed between all pairs of regions, except 289 between the Cauca river valley and the Caribbean coast (Table S2).

290

Due to the low level of haplotypic diversity, the structuration pattern observed with the two chloroplastic loci 291 was not clear-cut (Table 2, Fig. 2B). For trnL-trnF, haplotype H2 was observed in all regions and H1 in all

296

of trnL-trnF and rpl20-rps12 revealed two branches supported by bootstrap values >75% (Fig. 4). With the 297 exception of two individuals (one from BOL and one from MAR), the tree retrieved a geographical pattern 298 with the first group corresponding to the Caribbean coast and the Chicamocha canyon and the second group to 299 the Magdalena and the Cauca river valleys (Fig. 4). For both loci, AMOVAs attributed a majority of the 300 variance to partition among geographic regions (75% and 57% for trnL-trnF and rpl20-rps12, respectively; 301 Table 3). As for ISSR data, significant genetic differentiation ( ST) was also observed between all pairs of 302 regions, except between the Cauca river valley and the Caribbean coast, for the two chloroplastic loci (Table

314

Some of our results pointed to the existence of a pattern of spatial genetic structure in A. excelsum. Indeed, the 315 dendrogram obtained using the ISSR data grouped most of the localities of the inter-Andean valleys (except 316 MED; Fig. 2). In the same way, the concatenated sequences of the two chloroplastic loci characterizing 317 samples from the Magdalena and the Cauca river valleys form a well-supported branch in the UPGMA tree 318 (Fig. 4). Moreover, PCoA based on the ISSR dataset isolated SAN (i.e. Chicamocha canyon) and CES and 319 COR (i.e. Caribbean coast) from the rest of the populations sampled (Fig. 1). Taken together, our results 320 uncovered a slight genetic divergence between inter-Andean, Caribbean and Chicamocha canyon A. excelsum 321 remnants. Interestingly, for the two chloroplastic loci and the ISSR data, significant genetic differentiation 322 was observed between all region pairs except between the Cauca river valley and the Caribbean coast. We

323

propose that the existence of a past and/or current corridor of dispersal between these two regions may 324 explain our results. Higher genetic differences among regions were detected by AMOVAs based on the two 325 choloroplastic loci than when calculated with the ISSR data. The different evolutionary dynamics (i.e.

326

differences in strength of the genetic drift and migration and mutation rates) of nuclear and chloroplastic loci 327 could explain these contrasting results. Indeed, chloroplast 328 DNA (cpDNA) is uniparentally inherited and haploid and effective population size of chloroplast loci is four 329 times lower than the one of nuclear loci. Difference in type of inheritance could also lead to contrasting 330 migration rates between loci since cpDNA is only dispersed via seeds while nuclear loci are transmitted and 331 dispersed by both pollen and seeds. Moreover, choloroplastic loci are characterized by a low rate of 332 nucleotide substitution [START_REF] Hamilton | Four primer pairs for the amplification of chloroplast intergenic regions with 495 intraspecific variation[END_REF], while ISSR data present a high mutation rate [START_REF] Tuisima | Genetic diversity in Guazuma 573 crinita from eleven provenances in the Peruvian Amazon revealed by ISSR markers[END_REF].

333

Because of these differences among marker types, chloroplastic loci have been suggested to better reveal 334 historical changes, especially the population's response to Pleistocene glacial climatic variations, than nuclear 335 loci [START_REF] Ennos | Using organelle markers to elucidate the history, ecology and evolution of plant populations[END_REF]. Contrasting with values of ST between regions reaching up to 0.969 for the 336 choloroplastic loci, our ISSR data showed a high level of admixture between the three genetic groups detected 337 by Bayesian clustering (except in SAN, COR and CES; Fig. 3) and an increasing genetic differentiation 338 among populations with increasing geographic distance (i.e., a pattern of isolation by distance; Fig. S1). We

339

propose that these slight differences in patterns of genetic divergence detected by our two types of molecular 340 markers could be related, in part, to the complex scenario of a combination of historical isolation followed by 341 more substantial gene flow after secondary contact.

342

Historical separation between inter-Andean and Caribbean SDTF remnants has been previously supported by 343 the existence of floristic differences (Pizano and García 2014) and by a strong genetic structure within species 344 distributed in the Caribbean coast and the Cauca river valley (e.g. the tree Enterolobium cyclocarpum; 345 [START_REF] Thomas | Genetic diversity of Enterolobium cyclocarpum in 570 Colombian seasonally dry tropical forest: implications for conservation and restoration[END_REF]. The range and distribution of SDTF fragments, and concomitantly of the A. excelsum 346 populations, has been molded by the history of the region. Recent studies suggest that genetic isolation 347 between SDTF tree populations predate the late Pleistocene [START_REF] Caetano | The history of Seasonally Dry Tropical Forests in eastern 459 South America: inferences from the genetic structure of the tree Astronium urundeuva 460 (Anacardiaceae)[END_REF][START_REF] Thomas | Genetic diversity of Enterolobium cyclocarpum in 570 Colombian seasonally dry tropical forest: implications for conservation and restoration[END_REF], with 348 genetic differentiation building up in SDTF remnants isolated during the Neogene (Burnham and Carranco 349 2004; Côrtes et al., 2015). Secondary contact between differentiated genetic groups could have occurred 350 during the Pleistocene, a period during which the SDTF range is proposed to have been much larger 351 [START_REF] Pennington | Woody plant diversity, evolution, and ecology in the 552 tropics: perspectives from seasonally dry tropical forests[END_REF]. Indeed, the Pleistocene Arc theory suggests that the present-day SDTF patches represent 352 only relics of a much more contiguous formation extending all around the Amazon basin from the northeast of 353 Brazil to Paraguay, including the inter-Andean valleys of Bolivia, Peru, Ecuador and Colombia, up to the 354 Caribbean coast of Colombia and Venezuela [START_REF] Prado | Patterns of species distributions in the dry seasonal forests of South America[END_REF][START_REF] Pennington | Neotropical seasonally dry forests and Quaternary vegetation 550 changes[END_REF]. Even if the 355 Pleistocene Arc theory has been questioned [START_REF] Werneck | Revisiting the historical distribution of Seasonally Dry 580 Tropical Forests: new insights based on palaeodistribution modelling and palynological evidence[END_REF], SDTFs are still considered to have 356 occupied a much larger area in the past especially in the north of South America, the Andean region from 357 Bolivia to Peru and along a northern Amazonian block including areas of Colombia, Venezuela and the 358 Guiana highlands. In Colombia, the historically larger areas of SDTF were likely connecting our four regions 359 under study [START_REF] Pennington | Woody plant diversity, evolution, and ecology in the 552 tropics: perspectives from seasonally dry tropical forests[END_REF][START_REF] Thomas | Genetic diversity of Enterolobium cyclocarpum in 570 Colombian seasonally dry tropical forest: implications for conservation and restoration[END_REF]. The maximum extension of the SDTFs is suggested to 360 concur with the last Glacial-Holocene transition [START_REF] Werneck | Revisiting the historical distribution of Seasonally Dry 580 Tropical Forests: new insights based on palaeodistribution modelling and palynological evidence[END_REF]. Additionally, a recent palynological

363

Colombian Western Cordillera could be as recent as 4,000 yr before present (Munoz-Uribe 2012). In general, 364 patterns of isolation by distance equilibrate relatively quickly and rely mostly on recent demographic events 365 [START_REF] Barton | Inference in two dimensions: allele frequencies versus 448 lengths of shared sequence blocks[END_REF]. However, for species with long generation times such as some woody trees (trees are 366 generally considered long-lived perennials with lengthy juvenile period; e.g. in Anacardiaceae in Hormaza

367

and Wünsch 2007 and [START_REF] Dinesh | Wild crop relatives: 472 genomic and breeding resources[END_REF], it is possible that the fragmentation of SDTFs is, in our case, yet 368 too recent (< 4,000 yr) to detect its genetic footprint and explains why a pattern of isolation by distance can 369 still be retrieved in A. excelsum. Human activities during the last century, as for example the usage of tree 370 seedlings or transplantation from different localities (e.g. one tree from the botanical garden of BOL located 371 in the Caribbean coast had chloroplastic haplotypes typically encountered in the southern part of the country; 372 Fig. 4) or forestry activity could also have artificially increased gene flow between some isolated A. excelsum 373 remnants and led to a recent admixture between differentiated genetic groups.

374

For both kinds of markers (ISSRs and chloroplast sequences), the locality of SAN, within the Chicamocha

389

excelsum remnants in Costa Rica, a result possibly related to inbreeding depression in this partially self-390 incompatible tree [START_REF] Ghazoul | Reproductive ecology of tropical forest trees in logged and fragmented habitats 491 in Thailand and Costa Rica[END_REF]. Unfortunately, the markers used in our study are dominant 391 (i.e. ISSRs) or haploid (i.e. chloroplastic loci) and do not allow to test for possible traces of inbreeding in A.

392

excelsum populations. Thus, to facilitate the establishment of strategies that could mitigate risks of local 393 extinction in small isolated remnants of A. excelsum, complementary studies focusing on co-dominant 394 molecular markers such as microsatellites or SNPs (Single Nucleotide Polymorphisms) will have to be 395 undertaken.

396

Our study is the first attempt to provide genetic information at the country level for A. excelsum, a native 397 species representing a good candidate for local and national restoration programs of SDTFs. A. excelsum is 398 generally propagated by seeds since other methods, such as in vitro propagation or cuttings, have been proven 399 remarkably difficult [START_REF] Barreto | Efecto de microorganismos rizosféricos sobre germinación y 446 crecimiento temprano de Anacardium excelsum[END_REF]. No specific biological knowledge is required to grow A. excelsum 400 from seeds since the percentage of germination without any further treatment lies around 50% or higher 401 [START_REF] Barreto | Efecto de microorganismos rizosféricos sobre germinación y 446 crecimiento temprano de Anacardium excelsum[END_REF]

617

Individuals are arranged into populations from which they were sampled. C) Distribution of individuals with distinct level of admixture between genetic groups 618 across sampled localities of A. excelsum. The pie charts' color-code is: red = individuals assigned to the red genetic group with cluster membership of 0.60 or 619 higher, yellow = individuals assigned to the yellow genetic group with cluster membership of 0.60 or higher, blue = individuals assigned to the blue genetic group 620 with cluster membership of 0.60 or higher and white with black stripes = admixed individuals. Circles sizes are proportional to number of individuals. The code

621

of each locality corresponds to those noted in Table 1.

625

localities where encountered -directly followed by the number of sampled individuals presenting this sequence between parentheses -are given. Correspondence

626

with haplotypes presented in Table 2 for the trnL-trnF and rpl20-rps12 is given on the right for each unique concatenated sequence. Abbreviations of locality

627

names are noted in Table 1.

  ; Cárdenas and 103 Salinas 2007). Recent interest in the restoration of A. excelsum populations through translocation has emerged, 104 and the species has been included in the National Strategy for Plant Conservation of Colombia (NSPC) 105 (Pizano and García 2014). Indeed, A. excelsum is part of the red list of threatened Colombian timber trees 106 (Cárdenas and Salinas 2007) and is moreover characteristic to the SDTF, a habitat prioritized in the NSPC 107 (Pizano and García 2014). The high fragmentation of SDTFs and the isolation of A. excelsum populations 108 located in these remnants suggest that genetic differentiation could exist within the country. However, no 109 information about the genetic structure and diversity of A. excelsum is currently available, and no specific 110 molecular tools have been developed to work on these trees.

163

  Taq DNA polymerase (Invitrogen, Carlsbad, USA), 0.8μl of 2μM primer, and the remaining volume being 164 completed by pure water. The amplification products were visualized on a 2% agarose gel (w/v) including 2μl 165 of GelRed™ (Biotium, Fremont, USA) and separated by horizontal electrophoresis for 1h20. To ensure 166 consistent scoring of the bands, all 106 individuals were scored on the same gel for each ISSR.

167

  GeneRuler 100 bp (Thermo Scientific, Massachusetts, USA) was used as molecular weight marker. After 168 electrophoresis, agarose gels were photographed by Gel Doc NEQ-11 (Neon Creative Labware, Santiago, 169 Chile) under ultraviolet light.

214

  Chloroplast sequences were edited in CHROMAS (Conor McCarthy, Griffith University, Australia, 215 http://www.technelysium.com.au/chromas.html) and aligned in MEGA7 (Kumar et al. 2016). The number of 216 haplotypes (nH), the number of polymorphic sites (S), gene diversity (Hd, Nei 1978) and nucleotide diversity 217 ( , Nei and Li 1979) were computed using DnaSP (Rozas et al. 2003). Haplotypic richness (Re) was 218 estimated for each region using CONTRIB and the rarefaction procedure of Petit et al. (1998). Haplotypic 219 richness estimates reduce the effect of variation in sample size of A. excelsum among regions, with rarefaction 220 standardizing it to the smallest sample size. A rarefaction size of N = 6 was used, which corresponded to the 221 number of samples of the region with the lowest number of cpDNA sequences obtained (the Chicamocha 222 canyon; Table

  238 individuals were removed before analyses. For six individuals showing extremely poor DNA quality, no PCR 239 amplifications were possible. Thus, our ISSR analyses were performed on a total of 97 individuals, using 240 seven primers. Additionally, 89 sequences of trnL-trnF of 725-729bp and 92 sequences of rpl20-rps12 of 241 769bp were obtained.

  246 index (I) and Nei's gene diversity (He) values ranged between 0.23-0.52 and 0.16-0.35, respectively, and the 247 two estimators showed the same pattern within Colombia. The highest values of I and He were observed in 248 the localities of the Magdalena river valley (i.e. ARM, MAR, IBG, VEN and CHA; Table

  4 in the 262 Magdalena river valley, and 0.0 and 0.0 in Chicamocha canyon (values given for trnL-trnF and rpl20-rps12, 263 respectively).

  axes of the PCoA explained 42.97% and 23.62% of the genetic variability observed with the 267 ISSR data set for A. excelsum in Colombia (Fig.1). All localities from the Magdalena and Cauca river valleys 268 and the botanical garden sampled within the region of the Caribbean coast (BOL), clustered together. COR 269 and CES (Caribbean coast), and SAN (Chicamocha canyon) were isolated from the cluster of Magdalena, 270 Cauca and BOL (i.e. along axis 1 for the first two localities and along axis 2 for SAN). The closest localities 271 to SAN in the PCoA were CES and MED. The dendrogram based on Nei's genetic distance grouped the 272 localities from the south/center of the country (ARM, MAR, IBG, VEN, CLI, QUI, CHA) and BOL (Fig. 2C).

292

  regions except in the Chicamocha canyon. Haplotype H3 was less frequent and encountered only in the 293 Caribbean coast. For rpl20-rps12, H1 was encountered in all regions except the Chicamocha canyon and H2 294 in all regions except in the Caribbean coast. Haplotypes H3 and H4 were less frequent and only observed in 295 the Magdalena river valley and the Caribbean coast. The UPGMA tree based on the concatenated sequences

  303

For

  ISSR markers, our study revealed the existence of a gradient of genetic diversity within Anacardium 307 excelsum in Colombia, with the most diverse remnants encountered in the south (i.e. the region of the 308 Magdalena river valley and of the Cauca river valley), while the localities sampled in the Caribbean coast and 309 in the Chicamocha canyon were less diverse. Except for the Chicamocha canyon where the lowest haplotypic 310 richness was also observed for both trnL-trnF and rpl20-rps12, this pattern was not retrieved for the two 311 chloroplastic loci. Both chloroplast and nuclear markers supported the existence of genetic divergence 312 between distinct regions of Colombia. Caribbean coast exhibited unique haplotype diversity for trnL-trnF 313 chloroplastic locus while two ISSR genetic clusters seem mostly restricted to the southeast of the country.

  361and geochemical study of a site close to Mount Paramillo has estimated that vegetation change linked to 362 wetter and cooler climate, and then a possible recess of the SDTFs, in the northern termination of the

  375 canyon, presented the lowest level of genetic diversity for A. excelsum when compared to the other remnants 376 of Colombian SDTFs. Our analyses also provide evidence for possible genetic divergence of SAN from the 377 rest of the Colombian SDTF remnants of A. excelsum. During the last decades, this area has been shown to 378 suffer from high and persistent anthropogenic disturbance (Valencia-Duarte et al. 2012), a phenomenon 379 classically related to the existence of recurrent bottlenecks and low genetic diversity in plants (Young et al. 380 1996). We propose that anthropogenic disturbance could have increased the geographic isolations of relics 381 present in the Chicamocha canyon and that this may have ultimately limited gene flow between the natural 382 population of SAN and the rest of the Colombian SDTF remnants. The fragmentation of SDTFs has also led 383 to genetic isolation and inbreeding depression in other tree species in Brazil, such as Tabebuia ochracea 384 (Moreira et al. 2009), T. roseoalba, Handroanthus chrysotrichus, H. impetiginosus and H. serratifolius 385 (Collevatti et al. 2014). Inbreeding depression in plants could lead to a decrease in fecundity, in the number of 386 viable seeds as well as lower survival or lower resistance to stress. Identifying populations potentially affected 387 by inbreeding can be essential for the development of long-term conservation strategies (Keller and Waller 388 2002). An effect of forest fragmentation on pollination success and seed production was documented in A.

  . Restoration of SDTF is included in the National Strategy for Plant Conservation of 402 Colombia (NSPC), a recent environmental legislation aiming at the rehabilitation, recovery and conservation 403 of disturbed areas in Colombia (Pizano and García 2014; MiniAmbiente 2015). Herein, Regional 404 Autonomous Corporations (hereafter noted RAC) are in charge of the restoration and reforestation programs, 405 and have begun with the implementation of such plans for A. excelsum, as for example the RAC of 406 Cundinamarca (CAR) and the RAC of Tolima (CORTOLIMA) in the Magdalena river valley, and the RAC of 407 the Cauca Valley (CVC) and the RAC of Quindió (CRQ) in the Cauca river valley (Morales 2016, KT 408 Bocanegra-González pers. obs.). To meet the goals proposed in the NSPC, A. excelsum seedlings from 409 various nurseries and from different bioclimatic life zones (as defined by the Holdridge life zones system; 410 Holdridge 1947) are generally mixed up when reestablishing SDTFs, without considering the potential 411 introduction of new genetic variants and genetic pollution in the newly restored region, or the usage of seeds 412 that are poorly adapted to the new local conditions (MiniAmbiente 2015). Seedlings from mixed origins could 413 generate highly genetically diverse A. excelsum populations and potentially increase short-term resilience of 414 the reforested patches (Jump et al. 2009). However, this artificial gene flow can contribute to long-term 415 population decline via outbreeding depression and genetic swamping (see Hufford & Mazer 2003 for a 416 review). Because we do not yet know enough to predict which crossing scenarios have a risk for outbreeding 417 depression in A. excelsum, we advise against mixing seeds or plants from regions showing genetic divergence, 418 or using trees from off-regional sources in transplantation programs using A. excelsum. Finally, because of the 419 apparent genetic isolation and very low genetic diversity of A. excelsum in the Chicamocha canyon in the 420 state of Santander, we propose to prioritize this area within future conservation programs. 421 422

608 613 Figure 3 .

 6133 Figure 3. Bayesian assignment of individuals of Anacardium excelsum into genetic clusters using the program STRUCTURE based on seven nuclear inter-simple

Figure 4 .

 4 Figure 4. UPGMA tree obtained for the concatenated chloroplastic data set (i.e., joining trnL-trnF and rpl20-rps12) of Anacardium excelsum. Anacardium

  Figure 1 Click here to download Figure new Figure 1.pdf

  Figure 3 Click here to download Figure new Figure 3.pdf

Figure 4

 4 Figure 4 Click here to download Figure new Figure 4.pdf

  Table 1 and sequences from BOL were not included. As for ISSRs, all individuals from distinct sampling sites

	232
	233	within each region were pooled.
	234
		Results

Table 1 .

 1 Sampling localities of Anacardium excelsum and associated genetic diversity for seven inter-simple sequence repeat (ISSR) primers. For each locality, 585 the abbreviation (code), geographic coordinates, number of samples/number of samples genotyped with the seven ISSRs (Ns/Nsp), number of effective alleles 586 (Ne), percentage of polymorphic loci (%P), Shannon's index (I), expected heterozygosity (He) are displayed.

	587
	Region	Locality	Code	Geographic coordinates	Ns/Nsp	Ne	%P	I	He
	Magdalena river valley	Armero	ARM	5º00'42.59"N 74º54'39.35"W	15/14	1.50	89.9	0.43	0.28
		Mariquita	MAR	5°14'49.05"N 74°53'20.60"W	15/13	1.53	80.7	0.43	0.29
		Ibagué	IBG	4°25'35.47"N 75°12'46.66"W	16/16	1.52	91.2	0.46	0.30
		Venadillo	VEN	4°40'31.14"N 74°49'20.06"W	15/15	1.62	94.7	0.52	0.35
		Chaparral	CHA	3°45'10.53"N 75°32'26.29"W	6/4	1.45	71.9	0.39	0.26
	Cauca river valley	Medellín	MED	6°19'35.09"N 75°30'44.34"W	5/5	1.36	56.1	0.31	0.20
		Cali	CLI	3°26'52.31"N 76°33'30.44"W	10/10	1.50	85.9	0.45	0.30
		Armenia	QUI	4°30'43.23"N 75°39'03.07"W	5/5	1.43	75.4	0.39	0.25
	Caribbean coast	Pueblo nuevo	COR	8°28'22.14"N 75°30'36.71"W	3/3	1.42	57.8	0.34	0.23
		Turbaco	BOL	10°33'44.54"N 75°22'30.71"W	3/3	1.39	59.6	0.33	0.22
		Manaure	CES	10°13'56.89"N 73°02'03.08"W	3/3	1.30	40.3	0.24	0.16
		Corozal	SUC	9°17'16.62"N 75°13'44.68"W	4/0	-	-	-	-
	Chicamocha canyon	Bucaramanga	SAN	7°07'14.40"N 73°08'39.89"W	6/6	1.30	40.3	0.23	0.16
	Mean						69.9	0.38	0.23
	588

Table 2 .

 2 Genetic diversity of Anacardium excelsum for two chloroplastic loci (trnL-trnF and rpl20-rps12). For each locality, the abbreviation (code, as in Table1), In the haplotype list, the haplotype's name is directly followed by the number of sampled individuals carrying that haplotype, in parentheses.

	589
	590	number of sequences (N); number of haplotypes (nH); gene diversity (Hd); nucleotide diversity (π); number of polymorphic sites (S) and a list of haplotypes
	(Hap) are given. 591
							trnL-trnF						rpl20-rps12
		Region	Code	N	nH	Hd	π (.10 -3 )	S	Hap.	N	nH	Hd	π (.10 -3 ) S	Hap.
		Magdalena river	ARM	15	1	0	0	0	H1(15)	14	1	0	0	0	H1(14)
		valley
			MAR	14	2	0.14	0.98	5	H1(13),	14	2	0.14	0.37	2	H1(13),
									H2(1)						H2(1)
			IBG	16	1	0	0	0	H1(16)	14	1	0	0	0	H1(14)
			VEN	12	1	0	0	0	H1(12)	14	4	0.40	1.00	3	H1(11),
															H2(1), H3(1),
															H4(1)
			CHA	4	1	0	0	0	H1(4)	5	1	0	0	0	H1(5)
		Cauca river valley	MED	5	2	0.40	2.74	5	H1(1), H2(4)	5	2	0.40	1.04	2	H1(1), H2(4)
			CLI	7	1	0	0	0	H1(7)	8	1	0	0	0	H1(8)
			QUI	3	1	0	0	0	H1(3)	5	1	0	0	0	H1(5)
		Caribbean coast	COR	1	1	_	_	_	H2(1)	0	_	_	_	_	_
			BOL	3	3	1.00	4.57	5	H1(1), H2(1),	3	2	0.67	2.60	3	H1(2), H4(1)

Table 3 .

 3 Analysis of molecular variance (AMOVA) for Anacardium excelsum distributed across four Colombian regions for: A) seven nuclear inter-simple 596 sequence repeat (ISSR) primers, B) trnL-trnF chloroplastic sequences and C) rpl20-rps12 chloroplastic sequences. Regions were defined as in Table1:597Magdalena river valley, Cauca river valley, Caribbean coast, and Chicamocha canyon.Shades of grey correspond to the four sampled regions (Magdalena river valley: white, Cauca river valley: black, Caribbean coast: dark grey and Chicamocha 603 canyon: light grey; see Table1for more information on the localities and the abbreviation codes). Genetic differentiation among populations of Anacardium excelsum. A map including all the sampled localities is given on the left panel (A) and shades 606 of grey correspond to the four sampled regions (Magdalena river valley: white, Cauca river valley: black, Caribbean coast: dark grey and Chicamocha canyon: 607 light grey). The codes of each locality correspond to those noted in Table1. B) Median-joining networks for the trnL-trnF and rpl20-rps12 chloroplastic loci. In

	598
		Source of variation d.f. Sum of squares Variance components % total variance
		ISSR
		Among regions	3	62.417	0.633	6%
		Within regions	90	959.370	10.660	94%
		ST	0.056 (p = 0.012)
		trnL-trnF
		Among regions	3	36.415	0.868	75%
		Within regions	83	23.585	0.284	25%
		ST	0.753 (p < 0.0001)
		rpl20-rps12
		Among regions	3	11.460	0.254	57%
		Within regions	85	16.046	0.189	43%
		ST	0.573 (p < 0.0001)
	599 600	Levels of significance are based on 10,000 permutations Note that the locality of BOL, corresponding to a botanical garden, was not included in the analysis.

Acknowledgments

423

We would like to thank Fernando Fernández Méndez, Omar Aurelio Melo, Luis C. Acosta Cadenas, German 424 Urueña Serrano, Carlos A. Morales Carranza, Jeferson D. Galvis Jiménez, Sonia C. Camargo Roa, Jhon J.

425

Borda Velasquez, Ivan D. Vergara Terreros, Vannesa A. Montoya Sánchez, and Raul Rico Molina for their 426 assistance in the field and for providing samples of Anacardium excelsum from Colombia. We would like to 427 thank two anonymous reviewers and TGGE associated editor for their insightful and constructive comments.

428

This research was supported by the program "Talento Humano" N°043-16 from the Universidad del Tolima

432

The distinct haplotypes detected were deposited in GENBANK under numbers MG309720-MG309723 and 433 MG309724-MG309726 for rpl20-rps12 and trnL-trnF, respectively. This information is also available in the 434 materials and methods' section.