Mariano Alvarez

Julie Ferreira De Carvalho

Armel Salmon

Malika L Ainouche

Armand Cavé-Radet

Abdelhak El Amrani

Tammy E Foster

Sydney Moyer

Christina L Richards

C L R Orcid

Transcriptome response of the foundation plant Spartina alterniflora to the Deepwater Horizon oil spill

Keywords: Deepwater Horizon, Ecological genomics, Foundation species, Hydrocarbon response, Spartina alterniflora, Transcriptomics

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

Introduction

Human-induced environmental change is now a dominant evolutionary force, altering populations and ecological communities through direct and indirect effects of development and commerce [START_REF] Palumbi | Evolution -Humans as the world's greatest evolutionary force[END_REF][START_REF] Halpern | A global map of human impact on marine ecosystems[END_REF][START_REF] Gedan | Centuries of Human-Driven Change in Salt Marsh Ecosystems[END_REF]. Coastal ecosystems, which have great commercial value and provide valuable ecosystem services [START_REF] Pennings | Salt marsh communities[END_REF], are among the most vulnerable environments to human-mediated impacts (Halpern

Accepted Article

. 2008; [START_REF] Gedan | Centuries of Human-Driven Change in Salt Marsh Ecosystems[END_REF]). Within coastal systems, salt marshes, in particular, have long been the target of exploitation and development, and pollution arising from economic pressures threatens these productive communities worldwide [START_REF] Gedan | Centuries of Human-Driven Change in Salt Marsh Ecosystems[END_REF]). Pollutant releases, like oil spills, pose an immediate threat to organisms through chemical effects, such as polycyclic aromatic hydrocarbon (PAH) toxicity, and physical effects, such as coating [START_REF] Pezeshki | The effects of oil spill and clean-up on dominant US Gulf coast marsh macrophytes: a review[END_REF]. Over longer periods of time, pollution toxicity degrades habitat and potentially alters evolutionary trajectories in natural populations via selection and genotoxicity [START_REF] Lin | Impacts and Recovery of the Deepwater Horizon Oil Spill on Vegetation Structure and Function of Coastal Salt Marshes in the Northern Gulf of Mexico[END_REF]; [START_REF] Silliman | Degradation and resilience in Louisiana salt marshes after the BP-Deepwater Horizon oil spill[END_REF][START_REF] Reid | The genomic landscape of rapid repeated evolutionary adaptation to toxic pollution in wild fish[END_REF].

In 2010, the Deepwater Horizon (DWH) oil spill released an estimated 4.9 million barrels of oil into the Gulf of Mexico (National Commission on the BP Deepwater Horizon oil spill 2011). This oil eventually made landfall along 1,773 kilometers of the shorelines of Louisiana, Mississippi, and Alabama [START_REF] Mendelssohn | Oil Impacts on Coastal Wetlands: Implications for the Mississippi River Delta Ecosystem after the Deepwater Horizon Oil Spill[END_REF], nearly half of which was salt marsh habitat dominated by the grass Spartina alterniflora (Michel et al. 2012). As a foundation species, S. alterniflora provides crucial ecosystem functions by serving as refuge for invertebrates, as nurseries for birds and fish, and as a buffer from storm and wave action [START_REF] Day | Restoration of the Mississippi Delta: lessons from hurricanes Katrina and Rita[END_REF]; [START_REF] Mendelssohn | Oil Impacts on Coastal Wetlands: Implications for the Mississippi River Delta Ecosystem after the Deepwater Horizon Oil Spill[END_REF]. Spartina alterniflora also shows remarkable resilience to a variety of stressors [START_REF] Baisakh | Primary responses to salt stress in a halophyte, smooth cordgrass (Spartina alterniflora Loisel.)[END_REF]

Accepted Article

. However, the particular molecular mechanisms that regulate the remarkable resilience of S. alterniflora remain understudied. Understanding the molecular underpinnings of oil stress response in natura [START_REF] Shimizu | Plant sexual reproduction during climate change: gene function in natura studied by ecological and evolutionary systems biology[END_REF]; Kudoh 2016) may provide not only valuable information for the conservation and management of these threatened ecosystems, but also a novel understanding of the mechanisms of extreme stress response to an anthropogenic stressor.

Molecular biologists have developed the xenome concept to understand the genetic machinery that underlies the detection, transport and detoxification of toxic compounds, or 'xenobiotics' (Edwards et al. 2010). Although genetic resources are scarce in non-model species such as S. alterniflora, in the model plant Arabidopsis thaliana, accumulating evidence supports the involvement of six multigenic families in the activation, metabolizing, storage or excretion of xenobiotics: cytochrome P450s (CYPs), alpha/beta hydrolases, glycosyltransferases (GTs), glutathione transferases (GSTs), malonyltransferases (MTs) and ATP-binding cassette (ABC) transporters (Edwards et al. 2010; Fig. 1a). In the cell, xenobiotics are first activated by cytochrome P450s (CYPs) or alpha/beta hydrolases. Activated xenobiotics may be moved out of the cell, or conjugated with low molecular weight molecules such as glucose (via GTs) or the tripeptide glutathione (via GSTs).

Glycosylated xenobiotics are excreted or conjugated with malonate (via MTs)

Accepted Article

.

Methods

Study species and population

In August 2010, four months after the DWH oil spill, we collected plant material from four intertidal locations near Grand Isle, Louisiana and two locations in Bay St. Louis, Mississippi. These sites were naturally highly variable in conditions, but all sites supported monocultures of S. alterniflora and were moderately protected, as S. alterniflora grows predominantly in areas that are not exposed to full energy coastal conditions. In each of three contaminated and three uncontaminated populations of S. alterniflora, we collected leaf tissue from nine individuals spaced 10 meters apart (Fig. 2a). Contamination was confirmed by the visual presence of oil on the sediment in populations with substantial above ground dieback of S. alterniflora on the leading edge of the marsh. Nearby uncontaminated populations did not have any visual signs of the presence of oil or noticeable dieback of the above ground portions of S. alterniflora. From each plant, we collected the 3 rd fully expanded leaf to standardize age and minimize developmental bias in sampling. Leaf samples were immediately frozen in liquid nitrogen to prevent RNA degradation, and kept frozen during transport to the University of South Florida for processing and analysis.

RNA extraction and microarray hybridization

We extracted total RNA from each of nine plants separately per population from homogenized leaf tissue using RNeasy Plant Mini Kits (QIAGEN). The Interdisciplinary Center for Biotechnology Research at the University of Florida standardized RNA concentrations and created three pools of three randomly chosen individuals from within each population: a total of 54 samples were combined into 18 population-specific RNA pools.

Pooling is a common strategy in ecological research that sacrifices measures of individual-

Accepted Article

. level variation to increase sample size and capture population-level response [START_REF] Alvarez | Ten years of transcriptomics in wild populations: what have we learned about their ecology and evolution[END_REF]. We reverse transcribed twenty RNA pools (18 sample pools and 2 technical replicates) into cDNA and hybridized them to custom 4x44k Agilent DNA microarrays containing study were assembled with a minimum of 90% identity to account for the redundancy due to recent genome duplication events in these species, and were functionally annotated with sequence and protein homology searches using complementary datasets from Poaceae sequenced species, and A. thaliana (Ferreira de Carvalho et al. 2013). The samples used to create the transcriptome did not include plants exposed to oil, and therefore genes that are specifically up-or down-regulated during this stress response may be under-represented on the array. Probes were designed using the e-array web-portal (https://earray.chem.agilent.com/earray/) following the manufacturer requirements and printed on 4*44k chips with 2 sets of probes per 44k array. The microarray was designed to maximize the capture of overall expression of a transcript regardless of how many functional copies there are in the genome, which may have resulted in the collapse of homeologs into single probes.

Accepted Article

.

Microarray data analysis

We imported expression data into the statistical program JMP/Genomics (Version 6 for Windows; SAS Institute, Cary, NC, USA) for analysis after initial image processing by Agilent. We log-transformed the data and filtered out intensity values less than two and subsequent probes that appeared in less than half of the pools (<10). Of the total 17,049 probes, 15,950 passed our filtering protocol. Despite the potential for environmental variation in field studies, as well as possible variation introduced during the pooling process, our technical replicate pools were highly correlated with each other (r=0.995). We median normalized the raw data prior to downstream analysis, and visualized differences between individuals and populations using multidimensional scaling, implemented in the metaMDS function of the Vegan package [START_REF] Oksanen | vegan: Community Ecology Package. R package version 2[END_REF]).

We visualized the relative contribution of oil exposure, state (Louisiana or Mississippi), population, and slide using a principal variance components analysis (PVCA) on the probe-level data. This strategy uses a principal components analysis (PCA) to reduce the dimensionality of the data before calculating variance components via a linear model for each principal component [START_REF] Richards | Genome-Wide Patterns of Arabidopsis Gene Expression in Nature[END_REF]). To understand the effects of oil exposure, population and state on gene expression, we fit a mixed-model analysis of variance (ANOVA) to these data using a model that incorporated oil exposure, state, population, and microarray slide, with population nested within oil exposure and designated as a random effect (expression = oil + state + population-within-oil + slide). Populations were nested within oil because, due to the nature of the oiling event, we were not able to collect individuals affected by oil stress and individuals unaffected by oil stress from the same populations. State was included to model the effects of the comparatively large geographic distance between the four Louisiana populations and the two Mississippi populations, which also separated out on a visualization of the PCA of the expression data (Fig. S1). We

Accepted Article

. corrected for false discovery using q-values [START_REF] Storey | qvalue: Q-value estimation for false discovery rate control[END_REF], and used values of less than or equal to 0.05 to identify significantly differentially expressed transcripts.

Technical confirmation of microarray

Using the same RNA as we used in the microarray hybridization, we reverse transcribed total RNA using RetroScript kits (Ambion) to create cDNA for qPCR confirmation of the microarray. We generated primers with the primer3 program [START_REF] Rozen | Primer3 on the WWW for general users and for biologist programmers[END_REF] for three target genes as well as for a-tubulin, which has been validated as an endogenous control [START_REF] Baisakh | Primary responses to salt stress in a halophyte, smooth cordgrass (Spartina alterniflora Loisel.)[END_REF][START_REF] Ramanarao | Differential expression of the transcripts of Spartina alterniflora Loisel (Smooth Cordgrass) induced in response to petroleum hydrocarbon[END_REF]. We ran the qPCR reactions in duplicate using template from three individuals from each population (for a total subsample of 18). We quantified differential expression using the delta-delta Ct method with corrections for primer efficiency [START_REF] Schmittgen | Analyzing real-time PCR data by the comparative C-T method[END_REF], and confirmed the expression results of our DNA microarray for all three genes, which were up-or down-regulated in the same direction as they were in our microarray (Fig. S2).

Enrichment, stress annotation, and xenome comparison

To explore the functionality of genes differentially expressed in response to oil (Q<0.05), we used five different techniques based on both previous annotations and previously identified hydrocarbon-responsive genes. First, we performed a Gene Ontology (GO) enrichment test in JMP/G based on annotation data from the model species A. thaliana (TAIR 10), as measured by the Fisher exact test, to identify overrepresented GO terms within the oil-responsive genes (Q<0.05). Second, we searched for enriched Pfam annotations [START_REF] Finn | The Pfam protein families database: towards a more sustainable future[END_REF]) via a Fisher's exact test in R (Q<0.05) using HMMER [START_REF] Finn | HMMER web server: interactive sequence similarity searching[END_REF]) and custom R scripts (github.com/AlvarezMF/DWHoilspill_transcriptome).

Accepted Article

. Third, we compared transcriptome-wide response to crude oil in S. alterniflora to transcriptome-wide response to the PAH phenanthrene in two previous studies of A. thaliana. 24 hours identifying 467 genes that responded to phenanthrene in at least one of these time points. Combined, these studies captured gene expression response to phenanthrene that was affected in either the short term or the long-term in 1426 genes. While these two studies differed somewhat in their approach, plants were grown under similar abiotic conditions, and both assessed the response to phenanthrene of the same genotype of A. thaliana (Col-0). We looked for enrichment of gene categories with GO analysis on the 1426 differentially expressed genes in these two A. thaliana studies using the BioMaps function in Virtual Plant with a Fisher's exact test (FDR<0.05; [START_REF] Katari | VirtualPlant: A Software Platform to Support Systems Biology Research[END_REF]; File S1). We also looked for overlap between genes that were differentially expressed in S. alterniflora and homologs that were study. To facilitate the comparison with S. alterniflora, we only used genes that appeared on the custom S. alterniflora microarray and on either the ATH1 or CATMA chips. In total, we used 7,566 out of the 15,950 Spartina probes, all of which were homologous to A. thaliana loci and were represented on either the ATH1 or CATMA chips (Fig. S3). S1). We used χ² and Fisher exact tests (R Core Team 2015) to compare differentially expressed genes in A. thaliana and S. alterniflora by gene families.

Finally, to understand the relationship between oil-responsive genes in S. alterniflora and stress response genes in general, we compared oil-responsive genes to a previously generated list of environmental stress annotations for A. thaliana [START_REF] Richards | Genome-Wide Patterns of Arabidopsis Gene Expression in Nature[END_REF]). The stress annotation identified whether genes were previously shown to respond to abiotic or biotic stresses, including high light [START_REF] Rossel | Global changes in gene expression in response to high light in Arabidopsis[END_REF]

Accepted Article

. herbivores [START_REF] De Vos | Signal signature and transcriptome changes of Arabidopsis during pathogen and insect attack[END_REF]. Across these data sets, more than half of the genes on the ATH1 array had been associated with at least one stress (13,153 of 22,800 genes). For this comparison, we only used genes that appeared on both the custom S. alterniflora microarray and the ATH1 array on which the stress annotation was based (6770 genes).

Gene interaction networks

We used two forms of gene networks to identify specific genes that may be involved in crude oil response in S. alterniflora. First, we used a weighted gene co-expression network analysis (WGCNA; Langfelder & Horvath 2008) to assess relationships between genes and identify co-expressed clusters of genes that respond to crude oil exposure. To do this, we used the same filtered and normalized expression values that we used for our gene-by-gene linear mixed modeling (n = 15,950) and further filtered the gene set using the function goodSamplesGenes to remove 83 additional genes missing in exactly 50% of samples, which ensured that an adjacency matrix could be calculated for subsequent coexpression analysis (n = 15,867). We used the function pickSoftThreshold to quantify and plot values of scale independence and mean connectivity across a vector of possible soft thresholds in order to identify the optimal value of the soft threshold that minimized spurious noise within our network (Fig. S4a). Our scale-free topology index remained below 0.8 for all assessed values and we therefore used a soft threshold of 6 as recommended by the WGCNA creators. Next, we created a weighted adjacency matrix with the blockwiseModules function using biweight midcorrelation to represent the relationships between genes (Fig. S4a). The blockwiseModule function clustered the adjacency values into a dendrogram (Fig. S4b), which was separated using a dynamic cut algorithm to yield co-expressed clusters (File S2). For each cluster, the blockwiseModules function performed a PCA on the gene expression values for each cluster

Accepted Article

. and defined the first principal component of each cluster as the cluster eigengene, which represented the average expression of that cluster. The relationship of clusters to crude oil response was assessed using the generalized linear mixed model (eigengene = oil + state + population nested within oil) and corrected for multiple testing using a Holm correction. For this model, we removed the term for slide because it explained very little of the variance. We assessed Gene Ontology enrichment for each cluster with TopGO (Alexa & Rahnenfuhrer 2010), using genes with A. thaliana homologs as a background for Kolmogorov-Smirnov (KS) tests. We then summarized the most highly enriched GO term for each module, along with its corresponding q-value for a KS test (Table S2), and examined the protein sequences of genes within any oil-responsive clusters using HMMER [START_REF] Finn | HMMER web server: interactive sequence similarity searching[END_REF]) and custom R scripts (github.com/AlvarezMF/DWHoilspill_transcriptome). Finally, we used binomial tests to assess for enrichment of differentially expressed genes in each cluster.

We created a second gene interaction network using the software Virtual Plant [START_REF] Katari | VirtualPlant: A Software Platform to Support Systems Biology Research[END_REF]). We used the 1410 differentially expressed genes (Q<0.05) with homologs in A. thaliana, after removing duplicates and genes with no known connection (File S3).

Interactions, or edges, were parameterized with data on microRNA binding populations, protein-to-protein interactions, transcriptional regulation (which includes transcription factors, enhancers, and repressors), and transport interactions. The Virtual Plant software created additional edges using data from the metabolic interaction databases Aracyc and KEGG, and using published literature interactions [START_REF] Katari | VirtualPlant: A Software Platform to Support Systems Biology Research[END_REF]. Although these edges were parameterized with data from A. thaliana, we hypothesized that interactions between highly connected genes may be conserved, allowing us to identify important genes and pathways among the oil responsive genes in S. alterniflora. We visualized the resulting network using Cytoscape [START_REF] Shannon | Cytoscape: a software environment for integrated models of biomolecular interaction networks[END_REF]) and counted the number of connections using Virtual Plant without additional ranking. We used information from the WGCNA and Virtual

Accepted Article

. Plant networks to identify four genes that might be hubs, important regulators, or otherwise key genes in orchestrating the response to crude oil.

Functional confirmation in Brachypodium distachyon

We attempted to confirm the function of four genes in controlled experiments using the model grass Brachypodium distachyon. We chose four lines that were altered in target genes whose expression was correlated to hydrocarbon exposure in S. alterniflora with the understanding that these genes may represent the function of one or more potential homeologs within the S. alterniflora genome. We selected three genes (homologs of KCS11, ATMCB1, and ATTPS21) that were highly connected within our Virtual Plant interaction network, as measured by number of connections. Two of these genes, KCS11 and ATMCB1, were also part of the co-expression cluster that was correlated with oil exposure in our coexpression analysis (File S2). For the fourth gene, we identified an epigenetic regulator (SUVH5) that was correlated with oil exposure in the gene-by-gene analysis. SUVH5 was not part of the Virtual Plant interaction network or our oil-responsive co-expression cluster, but may participate in oil response by regulating gene expression.

We obtained seeds for wild-type (B21-3) and four previously generated and sequenceconfirmed T-DNA insertion lines (Fig. 3a) from the Western Regional Resource Center [START_REF] Bragg | Generation and Characterization of the Western Regional Research Center Brachypodium T-DNA Insertional Mutant Collection[END_REF], and stored seeds in ambient conditions before stratifying them for two weeks. We then sowed one replicate each of wild-type and T-DNA insertion lines in each of five oil treated trays, and five untreated trays, which were all grown in a single growth chamber. After all seeds had bolted, each treatment tray received 500ul of 2.5% crude oil in tap water every other day, which is a sub-lethal concentration that we found induced phenotypic response in S. alterniflora and in B. distachyon. Untreated trays received 500ul of

Accepted Article

. only tap water every other day. We grew seedlings for 27 days until a majority of plants had flowered and senesced. On the 27 th day, we measured the height and number of leaves and stems of each plant as proxies for above-ground biomass. We also measured the number of inflorescences as a component of fitness. Finally, we used a LI-6400 portable photosynthesis system (LI-COR) to measure photosynthetic rate and stomatal conductance, and to derive instantaneous water use efficiency (WUE).

Because treatment is applied at the level of trays, this design was a split plot and we modeled the effects of genotype and treatment on gene expression with a linear mixed model

Results

Oil-contaminated populations were differentiated from uncontaminated populations

Despite substantial differences among populations, oil exposure explained 25% of the variance in overall gene expression (Fig. 2b). In our principal components analysis, PC1 separated samples from the two uncontaminated sites in Louisiana from the samples collected in the uncontaminated Mississippi site and all samples exposed to oil (Fig. S1). The association of the uncontaminated Mississippi site with the contaminated sites in PC1 may be due to unobserved contamination, historical exposure, or underlying population structure. In fact, our principal variance components analysis indicated that differences among populations

Accepted Article

. explained the greatest proportion of variation in gene expression (32%, Fig. 2b). Using a gene-by-gene linear mixed model, we found that the expression of 3334 genes significantly differed by oil exposure, and 2287 genes significantly differed by state (869 genes overlapped the two main effects; Q< 0.05, Fig. 2c). We did not assess the significance of the random effect (population nested within oil).

Enrichment and comparison to PAH response in A. thaliana

Through gene set enrichment, we found three overrepresented categories among oilresponsive genes: "plant-type primary cell wall biogenesis", "sucrose biosynthetic process", and "coumarin biosynthetic process" (Q<0.05). Our Pfam enrichment test found three significantly overrepresented families among the oil responsive genes: "leucine-rich repeat", "Spc7 kinetochore protein", and "Cellulose synthase" (Q<0.05).

Of the common A. thaliana loci that were represented on our Spartina chip and either the ATH1 or CATMA chips (Fig. S3), 3063 (2417 unique A. thaliana genes) were responsive to the DWH oil spill in S. alterniflora and 526 were responsive to phenanthrene in A. thaliana (Fig. 4a). We found only 187 genes (224 homology-based S. alterniflora loci) that overlapped between these two species, including 12 genes annotated as involved in photosynthesis and 46 genes generally annotated as involved in biosynthetic processes (File S4). The discrepancy of the common number is due to the fact that multiple probes in S. alterniflora may be annotated with the same A. thaliana gene. Among these 187 genes, only 42 (22%) showed similar patterns of expression in S. alterniflora and A. thaliana, while most (78%) showed divergent expression patterns in S. alterniflora and A. thaliana (Fig. 4b), indicating little overlap between the two species' response. We found no significant overrepresentation of GO categories among genes expressed in the same direction, although we note that several Accepted Article . hydrolases and genes involved in oxidative stress appear among them. We also found that 2839 of the 3063 (93%) oil responsive transcripts in S. alterniflora (corresponding to 2230 of 2417 or 92% of the unique A. thaliana transcripts) were not responsive in either of the two previous studies of PAH response in A. thaliana.

When we assessed the overlap between the A. thaliana xenome and oil-responsive genes in S. alterniflora, we found six A. thaliana loci (seven homology-based S. alterniflora loci) that were differentially expressed in both A. thaliana and S. alterniflora (Fig. 1b,c), indicating that A. thaliana shares only 7.6% of the differentially expressed transcripts of the S. alterniflora xenome following PAH exposure. Again, the discrepancy of the common number is due to the fact that multiple probes in S. alterniflora may be annotated with the same A. thaliana gene. χ² comparisons by gene families revealed significant enrichment in α/β hydrolases and GT genes in the S. alterniflora xenome (P < 0.001, Table S1).

In a comparison of oil-responsive S. alterniflora genes to a previously generated "stress annotation" in A. thaliana, we found that our custom microarray had homologs for 3907 genes that were previously annotated for response to stress in A. thaliana, including 356 "super responsive" genes that responded to more than four stressors in the A. thaliana studies (File S5). Of these "super responsive" A. thaliana genes, 85 gene homologs (including six from the xenome) in S. alterniflora were also responsive to the DWH crude oil, highlighting the potential importance of this small subset of genes in a diverse array of responses that may be functionally conserved across taxa. Additionally, 149 "stress annotation" genes responded in at least one of the previous A. thaliana phenanthrene studies, including 18 "super responders", all of which also responded in S. alterniflora (File S5, Weisman et

Accepted Article

.

Gene interaction networks identify targets for confirmation

Using weighted gene co-expression network analysis (WGCNA; Langfelder & Horvath 2008) with the expression of 15,867 genes across the 18 replicate pools, we identified 39 co-expressed clusters of genes, ranging in size from 36 to 3974 genes (Table S2; File S2). Of the 39 co-expressed clusters, the value of the eigengene from one cluster (#3 in File S2) was significantly predicted by oil exposure after a Holm multiple testing correction. This cluster contains 2054 genes, which were significantly enriched for differentially expressed genes (665/2054, P<0.001). The only enriched GO category in this gene cluster after Q-value correction is "biological process". However, the cluster also contains genes with GO annotations that include RNA splicing, response to osmotic stress, response to salt stress, and plant-type cell wall organization (Fig. S5). This gene cluster is also enriched for a single protein family: RNA polymerase Rpb1 C-terminal repeat (PF05001, Q<0.05).

In Virtual Plant [START_REF] Katari | VirtualPlant: A Software Platform to Support Systems Biology Research[END_REF]), 6 of the top 20 most highly connected genes were part of the co-expressed cluster (#3) that was significantly associated with response to oil in our WGCNA analysis. However, many highly connected genes from our Virtual Plant interaction network did not fall within this cluster. For example, a gene homologous to ATMCB1 had 116 connections in our interaction network, but was part of a WGCNA module that was not significantly correlated with oil exposure.

T-DNA insertion genotypes differ in phenotypes from wild-type plants, and in response to oil

After exposing replicates of each of the T-DNA insertion lines to crude oil, we found overall treatment effects for height, number of leaves, number of stems, number of inflorescences, photosynthetic rate and water use efficiency (WUE), but not stomatal Accepted Article . conductance. We also found a significant effect of T-DNA insertion line for number of leaves, number of inflorescences, photosynthetic rate and WUE (Table 1). Using post-hoc assessments of significance for individual lines, we found a significant effect of a SUVH5 homolog (Bradi3g35330) overexpression on photosynthetic rate, and a KCS11 homolog (Bradi3g07730) overexpression on stomatal conductance and WUE (Table 2).

Finally, we compared response to oil in each line to the wild-type line and found significant T-DNA insertion line-by-treatment effects for two genes. We found that Bradi3g35330, which overexpresses a homolog of SUVH5, resulted in increased number of inflorescences in response to oil compared to the wild type (Fig. 3b; Table 2). We also found that Bradi1g62540, which overexpresses a homolog of ATTPS21, increased the number of leaves produced under crude oil exposure compared to wild-type lines in B. distachyon (Fig. 3c, Table 2).

Discussion

Large-scale anthropogenic impacts, such as crude oil spills, can be leveraged as natural "treatment and control" designs to understand molecular function in ecologically relevant settings [START_REF] Whitehead | Genomic and physiological footprint of the Deepwater Horizon oil spill[END_REF][START_REF] Reid | The genomic landscape of rapid repeated evolutionary adaptation to toxic pollution in wild fish[END_REF]. We found 3334 genes that responded to crude oil exposure in S. alterniflora, including 2230 genes that were not previously responsive in the model plant A. thaliana exposed to the PAH phenanthrene. We found that most (78%) of the shared, differentially expressed transcripts showed divergent expression patterns between S. alterniflora and A. thaliana. While the two species may share a broadly defined set of stress-responsive genes (including but not limited to the more narrowly-defined xenome), the specific oil-responsive transcriptional pathways in S. alterniflora, and the xenome in particular, diverge from those in A. thaliana.

Accepted Article

.

Transcriptional hydrocarbon response in S. alterniflora may be species specific

Plants may adapt to cope with toxic molecules [START_REF] Antonovics | Evolution in closely adjacent plant populations VIII. Clinal patterns at a mine boundary[END_REF], which may be reflected in the fine-tuning of metabolic detoxification pathways that are transcriptionally activated under challenging toxic conditions (El Amrani et al. 2015; [START_REF] Dumas | Unraveling the early molecular and physiological mechanisms Accepted Article . involved in response to phenanthrene exposure[END_REF]). In the Gulf of Mexico, natural oil seeps create the possibility of repeated hydrocarbon exposure in coastal species like S. alterniflora over evolutionary time [START_REF] Macdonald | Remote Sensing Inventory of Active Oil Seeps and Chemosynthetic Communities in the Northern Gulf of Mexico[END_REF]. We used the framework of the xenome to compare the regulation of metabolic pathways involved in organic xenobiotic detoxification in S. alterniflora and A. thaliana. We hypothesized that comparison of the xenome of S. alterniflora, which shows hyper tolerance in the presence of phenanthrene (Cavé-Radet & El Amrani, unpublished) and crude oil exposure in laboratory conditions and in the field [START_REF] Alvarez | Molecular response of Spartina alterniflora to the Deepwater Horizon oil spill[END_REF], to that of the PAH-sensitive model plant species A. thaliana exposed to hydrocarbons, would provide insight into the conserved mechanisms that underlie tolerance in plants, as well as nonoverlapping divergent components of tolerance. We found remarkably little overlap in genome-wide as well as xenome response between S. alterniflora and A. thaliana. Only 42 (22%) of the shared, differentially expressed transcripts showed similar direction of expression in S. alterniflora and A. thaliana, while most (78%) showed divergent expression patterns between S. alterniflora and A. thaliana. Less than 100 of the 3334 differentially expressed genes in S. alterniflora were among the six multigenic families that make up the A. thaliana xenome as strictly defined. Gene set enrichment on oil-responsive genes in S. alterniflora identified functional categories involved in metabolism (plant-type primary cell wall biogenesis, sucrose biosynthetic process) and defense (coumarin biosynthetic process) that were not enriched in previous A. thaliana studies, but may participate in the detoxification process more broadly (Edwards et al. 2010[START_REF] El Amrani | Omics" Insights into PAH Degradation toward Improved Green Remediation Biotechnologies[END_REF].

Accepted Article

. Highly connected genes are expected to evolve more slowly when they exhibit pleiotropic effects on phenotype [START_REF] Cork | The evolution of molecular genetic pathways and networks[END_REF]. We expected that highly connected stress-responsive genes would be functionally conserved and potentially play important roles in the response to crude oil. Therefore, we hypothesized that genes that are highly connected in A. thaliana, would appear within gene clusters that are most responsive to crude oil exposure, allowing us to identify genes that are more likely to modulate phenotype during crude oil response. However, we found little concordance between highly connected genes in our Virtual Plant network and oil-responsive genes identified in our coexpression analysis. In combination with the lack of overlap in our xenome survey, these results suggest that S. alterniflora has a substantially divergent genomic architecture underlying hydrocarbon response relative to A. thaliana.

The differences between the two species could be partly because of differences between the response to complex crude oil and response to phenanthrene. In the DWH oil spill, Napthalene was the dominant PAH (64% of the total PAHs) while phenanthrene made up only 17% [START_REF] Liu | The weathering of oil after the Deepwater Horizon oil spill: insights from the chemical composition of the oil from the sea surface, salt marshes and sediments[END_REF]. Additionally, the DWH contamination included other toxic compounds, such as methylnaphthalene and dispersants [START_REF] Liu | The weathering of oil after the Deepwater Horizon oil spill: insights from the chemical composition of the oil from the sea surface, salt marshes and sediments[END_REF], which may induce different transcriptional responses. However, although there were differences between the field and laboratory conditions, including xenobiotic composition, the hyper-tolerant phenotype of S. alterniflora may also be in part due to the fact that S. alterniflora detoxification involved more genes for detoxification than were identified in A. thaliana. The enrichment of S. alterniflora xenome response with more alpha/beta hydrolases and GTs may be related to the fact that, as a hexaploid, S. alterniflora has potentially many more functional or subfunctionalized copies of xenome genes [START_REF] Fortune | Evolutionary dynamics of Waxy and the origin of hexaploid Spartina species (Poaceae)[END_REF][START_REF] Boutte | Reference Transcriptomes and Detection of Duplicated Copies in Hexaploid and Allododecaploid Spartina Species (Poaceae)[END_REF]). The increased copy number of functional genes may also partly explain the overrepresentation of genes annotated as being responsive to cell wall biogenesis among oil-responsive genes in S.

Accepted Article

. alterniflora. However, enrichment of functional categories in S. alterniflora, particularly those involved in cell wall biogenesis in the GO analysis and cellulose synthase in the Pfam analysis, may be also related to a separate, species-specific detoxification process, which includes incorporation of the xenobiotic into cell wall polymers such as lignin (El Amrani et al. 2015) as well as the involvement of regulatory elements such as the Rpb1 protein family and CHH methylation, which may modulate gene expression (as discussed below). This behavior, to our knowledge, has not been documented outside of this study, and highlights both limitations in previously published annotations like Gene Ontology as well as the potential for novel behavior in non-model species or in natural settings [START_REF] Colbourne | The ecoresponsive genome of Daphnia pulex[END_REF][START_REF] Whitehead | Genomic and physiological footprint of the Deepwater Horizon oil spill[END_REF]). We emphasize, however, that our study compares the response of only two species, one of which is quite tolerant and the other sensitive to hydrocarbons, and which also differ in both physiology and evolutionary history. For example, monocots like S. alterniflora differ substantially from dicots, such as A. thaliana, in the composition of both primary and secondary cell walls [START_REF] Vogel | Unique aspects of the grass cell wall[END_REF], which may contribute to the observed enrichment of these genes in response to crude oil exposure. Ideally, experiments that include a gradient of hydrocarbon-tolerant and hydrocarbon-sensitive species across a broad phylogeny would assess the conservation of hydrocarbon response mechanisms across plant species.

Two candidate genes modulate conserved response to hydrocarbon stress

In our candidate gene screen using knockout lines of B. distachyon, we found that overexpression of a homolog of SUVH5 (Bradi3g35330), which was down-regulated in oil exposed S. alterniflora, resulted in increased number of inflorescences in response to oil compared to the wild type. SUVH5 is a methyltransferase that broadly contributes to non-CG methylation patterning in A. thaliana (At2G35160, [START_REF] Stroud | Non-CG methylation patterns shape the epigenetic landscape in Arabidopsis[END_REF]. Non-CG methylation, in particular CHH methylation, may be dynamically regulated in response to specific

Accepted Article

. Although the putative S. alterniflora SUVH5 was downregulated in response to crude oil compared to uncontaminated conditions in the field, our study suggests that individuals that can maintain relatively higher SUVH5 expression may display increased reproduction, and perhaps fitness, under oil stress compared to lines that do not. Alternatively, overexpression of SUVH5 may increase number of inflorescences, but reduce a more downstream component of fitness, such as germination success. Further tests with additional mutants will be useful in resolving this ambiguity, particularly with concurrent assessment of TE activity, and increased understanding of the function of DNA methylation variation. For example, manipulation of methylation variation, either through methylation inhibitors or the selective silencing of other methyltransferases, may help resolve the relationship between DNA methylation and fitness during crude oil stress [START_REF] Richards | Ecological plant epigenetics: Evidence from model and non-model species, and the way forward[END_REF]).

We also found that overexpression of a homolog of ATTPS21 (Bradi1g62540), which was up-regulated in oil exposed S. alterniflora, increased the number of leaves produced under crude oil exposure compared to wild-type lines in B. distachyon. ATTPS21 is a terpene synthase gene involved in the production of volatiles in flower petals in A. thaliana [START_REF] Liu | A Conserved Cytochrome P450 Evolved in Seed Plants Regulates Flower Maturation[END_REF]. Plant volatiles are often activated during stress, and may play a role in relieving oxidative stress [START_REF] Holopainen | Multiple stress factors and the emission of plant VOCs[END_REF]. Although the mechanism for this relationship remains unclear, the overrepresentation of the coumarin synthesis pathway among differentially expressed genes provides additional evidence for an interaction between Although we characterize gene function in this manuscript, further dissection of the genomic mechanisms of oil resilience within S. alterniflora will be required for understanding the evolutionary fate of oil-exposed populations of S. alterniflora. Our field sites possess an unknown ecological history, and latent, unaccounted-for effects may confound and bias our gene expression assays, making controlled follow-up studies important to accurately characterize the functional genomic response to crude oil. We found genetic differentiation among these same populations [START_REF] Robertson | Genetic and epigenetic variation in Spartina alterniflora following the Deepwater Horizon oil spill[END_REF], which may explain some of the observed response to oil in gene expression patterns observed in our study. Additionally, the observed variation in gene expression between oil-exposed and unexposed populations may be modulated by epigenetic regulators like SUVH5, which can exert effects in concert with or in addition to the effects of genetic variation [START_REF] Robertson | Non-genetic inheritance in evolutionary theory -the importance of plant studies[END_REF][START_REF] Richards | Ecological plant epigenetics: Evidence from model and non-model species, and the way forward[END_REF]. Finally, we note that our sample pooling strategy may miss a number of subtle or cell-specific transcriptional behaviors involved in hydrocarbon response. Future efforts to catalog gene expression variation using RNAseq and additional cell types (including rhizomes and roots) will enhance our capacity to understand the functional response to an increasingly common anthropogenic stressor in this foundation marsh species.

Considering that the resilience of S. alterniflora may attenuate many types of anthropogenic damage in salt marshes, understanding the molecular underpinnings of oil stress response may provide valuable information for the conservation and management of these ecosystems, as well as a novel understanding of the mechanisms of phenotypic plasticity.

Accepted Article

. on resident marsh fishes. Proceedings of the National Academy of Sciences of the United States of America, 109 [START_REF]Report to the President[END_REF]

Accepted Article

. Population was modeled as a random effect and not estimable. (a)

Accepted Article

  ; Baisakh & Subudhi 2009; Pennings & Bertness 2001; Silliman et al. 2012). Oil-affected populations showed up to 100% recovery within seven months of the DWH spill, despite the immediate effects of reduced carbon fixation and transpiration (Lin & Mendelssohn 2012; RamanaRao et al. 2012; Silliman et al. 2012; Lin et al. 2016), and evidence for genetic divergence of oil exposed populations from nearby uncontaminated populations (Robertson et al. 2017). Prior studies in S. alterniflora under controlled conditions have examined expression of candidate genes in response to heat, salt or oil stressors (Baisakh et al. 2008; Baisakh & Subudhi 2009; RamanaRao et al. 2012), and Bedre et al. (2016) recently characterized the full transcriptome response to salinity.

  Finally, xenobiotics are stored in the vacuole via ABC-transporters (El Amrani et al. 2015). Although studies have investigated the xenome response of A. thaliana exposed to one component of crude oil, phenanthrene (Weisman et al. 2010; Dumas et al. 2016), the xenome response to complex crude oil has not been quantified either in controlled conditions or in situ. Further, we know of no investigation of the xenome in non-model plants. However, targeted gene expression assays in S. alterniflora have identified several oil responsive candidate genes, including cell wall proteins and regulatory genes, which may participate in the detoxification Accepted Article . process more broadly (RamanaRao et al. 2012). The degree to which detoxification and tolerance pathways are conserved across evolutionary divergence remains understudied. Species-specific xenome features may reflect the differential ability to respond to challenging conditions under the selective pressure of toxic molecules (El Amrani et al. 2015). Therefore, measuring xenome and stress-response activity in non-model species under complex, natural conditions is critical to understanding the evolution of the xenome. To understand the mechanisms of response and resilience to oil stress in S. alterniflora, we assessed (1) differential gene expression in individuals from natural S. alterniflora populations exposed to the DWH oil spill, (2) the xenome and stress response of the oil-resilient S. alterniflora compared to that of PAH-sensitive A. thaliana, and (3) the phenotypic effects of four candidate oil response genes. Because native S. alterniflora populations harbor high levels of genetic diversity (Richards et al. 2004; Hughes & Lotterhos 2014; Foust et al. 2016; Robertson et al. 2017), we expected expression differentiation among populations. However, we also anticipated that large differences in expression due to oil exposure, particularly in the xenome and in genes involved in stress response pathways, would allow us to identify novel candidate genes since complex environments have exposed novel transcript behaviors in other studies (Colbourne et al. 2011; Whitehead et al. 2012), and S. alterniflora has a complex hexaploid genome (Ainouche et al. 2012) with potential for diversification of duplicated gene copies (Fortune et al. 2007; Roulin et al. 2013; Boutte et al. 2016).

Weisman

  et al. (2010) germinated seeds in phenanthrene and harvested the plants at 21 days, which identified 1074 genes that responded to phenanthrene. Dumas et al. (2016) ran a time series experiment exposing 15 day old seedlings to phenanthrene for 30 minutes, 2, 4, 8 and

  differentially expressed in either the Weisman et al. (2010) study or the Dumas et al. (2016)

Fourth, we compared

  the A. thaliana xenome with the putative S. alterniflora xenome by examining loci found in the six xenome gene families in Spartina (our study) and A. thaliana (Weisman et al. 2010; Dumas et al. 2016; i.e. represented on either ATH1 or CATMA chips). We defined a list of putative A. thaliana xenome genes, which was composed of 945 A. thaliana candidate loci (191 α/β hydrolases, 53 GST, 319 GT, 130 ABC transporters, 245 CYP450 and 7 related-malonyltransferases), based on published studies Accepted Article . (Edwards et al. 2010; Weisman et al. 2010; Skipsey et al. 2011; Dumas et al. 2016). Large families (GST, GT, ABC transporters and CYP450) were downloaded from the TAIR10 database (http://www.arabidopsis.org/browse/genefamily), whereas loci coding for α/β hydrolases protein superfamily were directly retrieved from the TAIR10 A. thaliana genome annotation file, and related-malonyltransferases from Taguchi et al. (2010). The in silico survey of the whole putative genes set of the xenome revealed that 239 candidate loci were represented in A. thaliana arrays (including 61 α/β hydrolases, 13 GST, 87 GT, 49 ABC transporters, 29 CYP450 and 0 related-malonyltransferases) related to 437 candidate loci in our S. alterniflora array (99 α/β hydrolases, 18 GST, 193 GT, 92 ABC transporters, 35 CYP450 and 0 related-malonyltransferases). We identified the phenanthrene responsive A. thaliana xenome by combining differentially expressed xenome genes in the Weisman et al. (2010) and Dumas et al. (2016) studies. We then compared only the loci whose homologs were present on the Spartina array and were also represented on either the ATH1 or CATMA. Hence, 22 loci met these criteria in A. thaliana (3 α/β hydrolases, 4 GSTs, 7 GTs, 4 ABC transporters, 4 CYP450 and 0 related-malonyltransferases) and 99 in S. alterniflora (26 α/β hydrolases, 6 GST, 49 GT, 14 ABC transporters, 3 CYP450 and 0 relatedmalonyltransferases, Table

  implemented in R (R Core Team 2015): (Response ~ Treatment + Genotype + (Treatment * Genotype) + Error (Block*Treatment)), with block as a random effect (Richards et al. 2008).We assessed effects using a Type II Wald Chi-squared test implemented in the Car package[START_REF] Fox | An R Companion to Applied Regression, Second Edition[END_REF]) in R. We interpreted significant genotype by treatment interactions to indicate that manipulation of the gene of interest affected the response to crude oil.

  stressors and may contribute to the regulation of transcription under stress or novel environments (Dowen et al. 2012; Dubin et al. 2015). CHH methylation is enriched near transposable elements (TE), and the maintenance of TE silencing may prevent deleterious insertions in protein-coding regions of the genome (Matzke and Mosher 2014). Alternatively, CHH methylation within TEs can also regulate nearby genic regions (Matzke and Mosher 2014), which would allow SUVH5 to function as a global regulator of stress response.

  the production of volatiles and fitness during crude oil exposure in both B. distachyon and S. alterniflora. Conclusions While these findings show novel patterns of response to hydrocarbon stress, and our experimental results support a role for SUVH5 and ATTPS21 in regulating the response to crude oil, the divergence between S. alterniflora, B. distachyon, and A. thaliana limits the conclusions that we can draw. The Spartina genus (Chloridoideae subfamily) diverged from Brachypodium (Pooideae subfamily) more than 40 million years ago (Rousseau-Gueutin et al. 2015), and from Arabidopsis approximately 150 million years ago (Kumar et al. 2017). During that time, the Spartina genus has undergone several hybridization and polyploidization events (Ainouche et al. 2004) that have increased chromosome number and ploidy level from 2n=4x=40 (tetraploid) to 2n=12x=120-124 (dodecaploid). Genome duplication and copy number variation can complicate our understanding of the role of specific genes, and increase the proportion of false positive matches (Fortune et al. 2007; Ferreira de Carvalho 2013; Primmer et al. 2013; Boutte et al. 2016; Ferreira de Carvalho et al. 2017). Although we characterized gene function using B. distachyon homologs, confirmation of gene function in S. alterniflora would require detailed forward and reverse genetic screens, using both transcriptome-wide and single-mutant assays in both field and controlled conditions, to fully characterize gene function in response to crude oil exposure. However, leveraging the genetic resources of model organisms provides a useful baseline and comparison to describe gene function during oil exposure in S. alterniflora. Additionally, genomic assays in non-model species that inhabit diverse ecologies will be insightful in the effort to characterize the behavior of genes in natural settings, providing a more complete Accepted Article . "ecological transcriptome" or "ecological annotation" (Landry & Aubin-Horth 2007; Richards et al. 2009; Alvarez et al. 2015; Kudoh 2016).

:

  Fig. 1: Xenome comparison between Spartina alterniflora and Arabidopsis thaliana. (a) Simplified scheme of metabolic pathways of xenobiotic (X) detoxification in plant cells: the thick line indicates the cell membrane, xenome components are indicated in red. (b) The total number of xenome genes expressed in S. alterniflora in response to oil contamination in the field, compared to those observed for Col-0 shoot tissue in response to the PAH phenanthrene (Weisman et al. 2010; Dumas et al. 2016). (c) Graphical matrix depicting differentially expressed xenome loci in response to phenanthrene in A. thaliana (Weisman et al. 2010; Dumas et al. 2016), and in response to the DWH oil spill in S. alterniflora. Arabidopsis thaliana based gene IDs are categorized by the five xenome gene families on the left. Changes in gene expression were scored as a binary trait: no significant response to PAH are indicated in black, differentially expressed genes are indicated in red.

Fig. 2 :

 2 Fig. 2: Response to the DWH oil spill in S. alterniflora across Gulf of Mexico study sites. (a) Green markers represent sites with no visible oil and red markers represent sites with visible oil in the sediment and substantial dieback of above ground Spartina alterniflora stems. (b) The percent of the transcriptional variance explained by oil, state, population, slide and the residual. (c) The number of genes that significantly differ between populations exposed to oil or not, and location (State).Population was modeled as a random effect and not estimable.

Fig. 3 :

 3 Fig. 3: Brachypodium distachyon T-DNA insertion lines used for confirmation of candidate genes. (a) T-DNA lines are listed as the catalog numbers from the WRRC with corresponding S. alterniflora contigs from our custom microarray and the closest A. thaliana homolog. Modification type represents the effect of the T-DNA insertion, either knocking out the function of the particular gene or tagging the promoter region to induce overexpression. Expression response in S. alterniflora is listed along with the justification for choosing each line. Reaction norms for (b) number of inflorescences, and (c) number of leaves for each T-DNA knockout line, referenced by its homolog in A. thaliana, in oil treatment and control.

Fig. 4 :

 4 Fig. 4: Comparison of oil responsive genes across studies in S. alterniflora and A. thaliana. (a) From a consistent set of genes surveyed (represented either on the Spartina array and either the ATH1 or CATMA arrays), the total number of genes expressed in S. alterniflora in response to oil contamination in the field, compared to those observed for Col-0 shoot tissue in response to the PAH phenanthrene (Weisman et al. 2010; Dumas et al. 2016). (b) Direction of difference in expression in 187 genes that respond to hydrocarbon stress in both S. alterniflora and A. thaliana.

  , cold, drought, heat, osmotic stress,

oxidative stress, salt, genotoxins, UV-B exposure, wounding (Kilian et al. 2007), infection by RNA virus (Whitham et al. 2003; Babu et al. 2008), bacterial pathogens, fungi and

Table 1 :

 1 , 20298-20302. doi: 10.1073/pnas.1109545108 85. Whitham, S. A., Quan, S., Chang, H. S., Cooper, B., Estes, B., Zhu, T., … Hou, Y. M. (2003). Diverse RNA viruses elicit the expression of common sets of genes in susceptible Arabidopsis thaliana plants. Plant Journal, 33(2), 271-283. doi: 10.1046/j.1365-313X.2003.01625.x Results of linear mixed models and generalized linear mixed models for phenotypes assessed. Height, photosynthetic rate, stomatal conductance, and water use efficiency were modeled separately with a normal distribution, while number of leaves, number of stems, and number of inflorescences were modeled with a Poisson distribution. Effects were assessed through Wald type II chi squared tests. Marginal and conditional R 2 values were assessed through the piecewiseSEM R package. df indicates degrees of freedom.

	Tables
				Oil	Line	Line x Oil
				treatment		treatment
				df = 1	df = 4	df = 4
		Marg.	Cond.	Chi Sq	Chi Sq	Chi Sq
		R 2	R 2
	Final height	0.27	0.42	6.55 *	3.45	4.32 NS
	Leaf number	0.43	0.56	9.15 **	11.14 *	7.08 NS
	Stem number	0.33	0.42	8.26 **	5.51	3.59 NS
	Infl number	0.35	0.54	8.27 **	1.46	5.41 NS
	Photosynthesis 0.38	0.65	4.87 *	21.65 ***	7.13 NS
	Stom conduct	0.23	0.48	0.14	18.29 **	2.43 NS
	WUE	0.25	0.34	4.80 *	6.54	3.87 NS

* P < 0.05; ** P < 0.01; *** P < 0.001; NS P > 0.05.

Table 2 :

 2 Significant terms in pairwise comparisons of T-DNA insertion lines to wild type, and interaction of line x treatment with T values (for linear mixed effects models) and Z values (for generalized linear mixed effects models) using Satterthwaite approximations to estimate degrees of freedom.

	Effect	Phenotype	Estimate	t/z
	Line (Bradi3g07730)	Stomatal conductance	-0.24	-2.15*
		WUE	15.10	0.013*
	Line (Bradi3g35330)	Photosynthetic	5.74	2.54*
		rate
	Interaction of Line x	Leaf number	0.45	2.10*
	Treatment (Bradi1g62540)
	Interaction of Line x	Inflorescence number	0.53	1.97*
	Treatment (Bradi3g35330)

Acknowledgements

We thank Steve Pennings, Brittney DeLoach McCall, Aaron Schrey, Christina Meals, and Ashley Shayter for access to field sites and assistance with plant sampling. This work was supported by funding from the University of South Florida New Investigator Grant program, the National Science Foundation (U.S.A.) DEB-1419960 and IOS-1556820, and the Franco-American Fulbright Commission (to CLR). Sequencing analyses were supported by Genoscope (Evry, France), the Environmental Genomics and the GenOuest (Bioinformatics) Biogenouest platforms (University of Rennes 1, France) and The Partner University Funds.

Author Contribution

CLR, AS & MLA designed the research; CLR, MA & SM collected plants, performed greenhouse studies and data collection; AS, JFC, MLA designed the Spartina microarray; CLR, MA, AC-R, AEA & TEF performed data analysis. All authors contributed to data interpretation, and writing of the manuscript.

Data accessibility statement

Unprocessed and processed microarray data have been deposited at the Gene Expression Omnibus (GEO) and are available at accession number GSE114635. Other processed data files, R scripts, and site coordinates are available on Dryad at doi:10.5061/dryad.kp0rc5s and on GitHub at https://github.com/AlvarezMF/DWHoilspill_transcriptome/.

Supporting Information

Table S1: Comparison of differentially expressed xenome genes in A. thaliana and the custom S. alternflora.

Table S2: List of 39 co-expressed gene clusters identified by WGCNA, with the number of genes in each cluster and the most significant enrichment category from Gene Ontology analysis (For clusters with equally significant GO categories, all most significant categories are shown). File S5: Stress-responsive genes (n = 3907) across stresses, the A. thaliana xenome and S. alterniflora oil responsive genes.


M Al Heib

Feedback of the empirical approach to design the room and pillar mines -Application on chalk mines (France)

The goal of this study is to improve the analysis of the stability of pillars and to answer the question of the choice of the pillar strength using an empirical approach. The pillar strength depends mainly on two factors: the uniaxial compression strength and the pillar geometry. We have shown in this study the importance of geometry of pillars by comparing calculations with and without the effect of W/H ratio (slenderness). The sensitivity of the parametric study based on the empirical relationships converges to the same resistance of the pillar despite the difference between the coefficients of these empirical relationships. The methodology was applied to case studies (France). The application of these relationships can help to explain the collapse of Clamart mine (Paris). Generally, and in the absence of data, it will certainly be more prudent to use the relationships available taking into account the effect of uncertainty.

INTRODUCTION

The room and pillar mining method is conventionally used for shallow mines. The method consists in practice to create rooms separated by small parts of the deposit, deliberately not exploited (pillars). The depth is controlled by the deposit; the extraction ratio should be selected according the bearing capacity of pillars. The pillar function is to ensure local and overall stability of the surface and underground mine. The stability of mines should be checked for short-term, long-term, local (isolated pillar) and large scale.

Several methods (approaches) are used to design the pillars, from the simplest to the most sophisticated one. The design is generally based on empirical methods from lessons learned, sometimes with dramatic consequences due to large collapse. Geologists and engineers nowadays have more scientific approaches to assess the stability of these structures. They are based on analytical and numerical calculations, sometimes in 3D and use sophisticated models with numerous parameters used to evaluate the short and long term pillar behavior.

All approaches should assess the risk through a comparison between the stress of the pillar and strength. The simplest analytical approach is often preferred. It is enriched by the feedback. The design of the pillars passes by estimating the average stress on the vertical pillar and comparing it with the strength of the pillar. The safety factor corresponds to the ratio of the resistance of the pillar (Rp) with respect to the stress of the pillar (σp).

It requires great care in the choice of parameters to be taken into account to establish the design and the fair and accurate diagnosis based on uncertainties of underground conditions. Moreover, the choice of safety factor should take into account the presence of discontinuities, age etc. In case we retain a safety factor near one, the failure probability could increase significantly. The choice of an optimum safety factor value depends on the experience and extensive knowledge of the studied area.

The risk of instability faced by room and pillar mining method has different shapes and match a variety of failure mechanisms. The historical accidents in underground mines are often associated with the failure of the pillars (ex. Coalbrook collapse in South Africa, Clamart in France, etc.). The failure of the pillars is due to:

• the mechanical properties due to an insufficient intrinsic compression resistance; • the buckling of the pillars becoming too slender or such as a result of roof falls; • the geological structural and discontinuities may affect isolated pillars or whole districts of the mine; • the exogenous factors can significantly reduce the resistance, it is primarily water and time.

Sudden collapses risk increases with increasing height pillar relative to its width [START_REF] Mark | The evolution of intelligent coal pillar design 1981-2006[END_REF].

In this paper, we will discuss the empirical relationship to estimate the safety factor of room and pillar mining method.

PILLAR VERTICAL STRESS

The excavation of the chambers changes the state of initial stresses in the massive and transfer the initial stress to the pillars. The pillar stress can be calculated using an analytical method to estimate the average vertical stress on the pillars: it depends on the depth of the pillar (h), density of the rock mass (ρ) and extraction ratio (E), g gravity acceleration This approach gives the average vertical stress and it does not inform about the spatial variation of this vertical stress through the pillar section and the other components of the stress tensor. The method may be used to estimate the stresses in the pillars of operating at several mine levels, mainly when the pillars are perfectly superimposed. Indeed, if the pillars are poorly superimposed and if the thickness of the parting bed is insufficient, the approach of tributary area cannot be used to estimate the vertical pillar stress and the state of stability.

PILLAR STRENGTH

The strength and behavior of the pillar is often linked to the behavior of specimens tested in the laboratory. Pillar strength is characterized from the uniaxial compressive strength (Rc) of the rock determined by a uniaxial compression test on a sufficient number of samples. Pillar resistance is also related to the geometric shape and the slenderness of the pillar (W/H): width (W) and height (H). Based on these parameters, several empirical equations have been proposed to estimate the relationship between the issue value of uniaxial compression tests in laboratory R c and the pillar strength of R p . The most general form of the equation is a mathematical power function.

λ, α and β are numerical constants, they are determined through feedback and in-situ observations of the behavior of the stable and collapsed pillars. The first known relation was established after the disaster in the Coalbrook Coal Mine (South Africa), in which more than 900 pillars were collapsed and the loss of 437 lives. The relation takes into account the effect of scale and heterogeneity of the massif. Table 1 provides the values of these coefficients available in bibliographic references for different rock types (coal, limestone, etc.). [START_REF] Salamon | A study of the strength of coal pillars[END_REF] have proposed coal pillar design relation from more than 125 historical cases. [START_REF] Hedley | Stope and pillar design for the Elliot Lake Uranium Mines[END_REF], have shown that the equation proposed by Salmon and Munro can also be used for hard rock, checked for mines of quartzite with a uniaxial compression strength of 230 MPa. The analysis was conducted on an area of 60 pillars out of which 23 were collapsed. Esterhuizen et al. (2007) and others have established a relation for limestone mines. Their approach takes into account the conditions of the pillars, the state of the rock mass Figure 1. Evolution of pillar strength (Rp)-uniaxial compression strength (Rc) against W/H ratio according to different empirical approaches (see table 1).

fracturing. They proposed a sizing equation from thousands of stable pillars observed. [START_REF] Martin | The strength of the hardrock pillars[END_REF] used the Hoek-Brown criterion to study the stability of pillar and to take into account the effect of the pillar dimension and the width-to-height ratio.

It is noted that the influence of λ (related to uniaxial compression strength Rc) is considered linear, however the influence of the pillar geometry follows power function. [START_REF] Bieniawski | The effect of specimen size on compressive strength of coal[END_REF] noted that scale effect is a power function, the strength decreases with the increasing of the sample dimensions. Spanish iron compression strength was determined in laboratory, for different height/diameter ratio (H/D), ranged from 5 to 0.5. The strength follows a power function comparable to those proposed for estimating the strength of the pillar [START_REF] Oyanguren | La stabilité des mines de fer de la région de Bilbao (Espagne)[END_REF].

Based on the parameters of the Table 1, Figure 1 shows the evolution of pillar strength/compressive strength ratio (Rp/Rc) for different values of slenderness (H is assumed to be constant and equal to 5 m). At first glance of the coefficients (α, β and λ) are all different because they have been established for different rocks. Nevertheless, we note that the results are close when W/H ratio is lower than 0.5 or smaller width of pillars (or tall). On the other hand, the gap is widening after this value. We compare the previous results with the effect of scale with the following values (α = 1, β = 1 and λ = 0.93).

It is clear that the equations take into account in-directly the effects of confining and buckling associated with the geometry of the pillar. Note that most low slenderness (W/H>> 1), more horizontal confining stress in the pillar is important. The increase in the horizontal stress tends to increase the capacity of the pillar supporting the vertical load, which results in a reduction of buckling failure mechanism.

Rp pillar strength must be greater than the pillar stress, taking into account indirectly the geotechnical and geometrical parameters related to pillar design.

The use of a three of four equations can be conceded as acceptable in first step of the analysis (Table 1), especially for configurations and rocks with similar characteristics. One could even make the calculation for the three hypotheses (Salamon and Munro, 1968, Herdley et al. 1972[START_REF] Esterhuizen | Pillar strength and design methodology for stone mines in the United States[END_REF]) and consider the average value.

Discussion and comparison

Numerical methods are very useful tools to understand complex rock and structure behavior. [START_REF] Laouafa | Modélisations des exploitations de sel[END_REF] considered that the method of the tributary area ignores friction action between the immediate roof and pillar.

The strength of the pillar increases because of friction forces between the immediate roof and pillar. He introduced a coefficient Cp of a variable value according to the pillar geometry and the mechanical properties of the mine layer (width, height, friction angle, cohesion). Laouafa carried out a simple numerical model to determine the evolution of the reduction/increased coefficient based on the pillar slenderness. The effect of the geomechanical properties (cohesion, friction) are well below the effect of pillar geometry. For pillars with (W/H) higher than 3, the pillar strength can reach three times the uniaxial compression strength. For (W/H) very low, the pillar strength is equals to uniaxial compression strength. Figure 2 compares the numerical relationship with the empirical relation of Salamon (α = 0.46, β = 0.66 and λ = 0.93) and the linear relation (α = 1, β = 1 and λ = 0.93). This comparison calls for the following remarks: the three relationships confirm that the increasing of (W/H) ratio increases the strength of the pillar; the pillar strength can be more or less than 1 according to the empirical approaches, however for the numerical approach, the pillar strength is greater or equal to the compression strength of the rock.

The strength of the pillar continues to increase according the "linear" approach, somewhat less for numerical approach, however it stabilizes around 1.10 for Salamon approach. Van der Merwe and Madden (2002) consider the relationship of Salamon and Munro tends to underestimate the strength of the large pillars and overestimate the resistance of small pillars. In fact, the numerical modelling with numerous parameters can certainly reflect the real behavior but due to time consuming issues, it may have limited their engineering application.

Effect of pillar

The influence of the width and height of the pillar to the same slenderness or equal to 1 was studied, therefore W = H = 1 and for different values of L and H. We used assumptions equations 1 to 3 (Table 1). Figure 3 shows the Rp/Rc ratio. The large pillars are, according to these relationships, weaker than the small pillars and short. The relationship is not linear and the dimensions of the pillar regardless of its slenderness plays an important role in determining its strength.

Variability of parameters

The determination of parameters, λ, α and β was made from data analysis. We do not have the variability of these parameters. To study the variability of such influence on the estimation of the strength of pillars, we adopted three parameters (α, β and λ) with a normal distribution characterized by a mean value (m) and standard deviation (s). The standard deviation is 10% of the average value. We analyzed this analysis for different W/H ratio values (H = ct.).

Figure 4 presents three curves expressing the ratio estimating strength of the pillar relative to the uniaxial compressive strength (mean, mean + standard deviation, mean-standard deviation). The standard deviation variation of Rp/Rc is not linear, it varies between 12 and 20%. For the design of an operating room and pillar area, it is safer to consider the lower limit (mean value minus the standard deviation).

This sensitivity parametric study can also be applied to all the variables in the empirical equation.The uncertainties on the pillar dimensions (width and height) and pillar strength lead to lower values (and higher) of the ratio (Rp/Rc) and can impact the design of pillars.

Conclusion

The methodology adopted and discussed herein takes into account the effect of scale on rock strength and pillar dimensions. The average stress of the pillar is determined by the method of the tributary area. Pillar strength is based on the simple compression standard test specimens and empirical coefficients related to the geometry of the pillars and can be determined by feedback analysis. The necessary data is limited to geometrical and geotechnical data from the underground mine. The uncertainty of the parameters can be considered thanks to sensitivity parametric studies.

FEEDBACK OF THE CLAMART CHALK MINE

The Clamart chalk mine was operated until 1880 by the method of abandoned room and pillar, [START_REF] Heib | Analysis of the historical collapse of an abandoned underground chalk mine in 1961 in Clamart[END_REF]. The disastrous collapse of the chalk mine of Clamart and Issy-les-Moulineaux, two bordering suburbs of Paris, on June 1, 1961 is one catastrophic example. The collapse of the Clamart mine caused the brutal subsidence of three hectares of ground surface and the destruction of a whole urban district, including residential buildings, several roads and sports facilities. Twenty-one people died and more than fifty people were injured.

The Clamart collapse was associated with two seismic events, 7 seconds between them, one at 10:30 am, when the soil subsided and loud noises had been heard, and the second collapse was recorded a halfhour later. Actually, the existing data of the events cannot be analysed today in order to obtain more information concerning the magnitude, the energy and the localization.

Two chalk levels were extracted separated by thin layer. In the southern part, the upper level has been exploited to the south and north with a different development scheme but in both cases very regular. The junction between the two zones takes place on a row of irregular and hybrid pillars. The geometric characteristics of the different areas are shown in Figure 5, they are also given in Table 2. The ratio (W/H) is varied between 0.73 (ZNS1) and 1.4 (ZSS1).

The geometry of the pillars of the lower level is different from the upper level by its irregularity, that could result from an early cessation of mine operations.

To understand the origin of the collapse, we compare first the pillar stress (<2 MPa, see table 2) to laboratory uniaxial compression strength (2 MPa), one can consider the pillars are very close to stability. We analyzed then the pillar stability using Salamon and Munro relation to calculate the strength of pillar taking into account the (W/H) ratio, considering that coal strength is the nearest strength of chalk. According to this choice, all areas of Clamart mine had a safety factor less than 1, except the area ZSI3. The ZNS1 area (upper level) would even have a value of 0.41 and note that all the pillars of the sector were very degraded and had probably been reinforced during the mining operation (of black lines in Fig. 5).

In light of this analysis, the mine could not be stable, including the top level. The collapse took place well after the mine operation and the value of the compressive strength of the chalk remains unknown during operation (short-term). The evaluation of the stability of the Clamart mine is essentially based on average values of the width, height, the compressive strength of the chalk and the coefficients (α, β and λ). The compressive strength of the chalk was measured on ten samples; its value varies relatively little. We considered a normal distribution with a mean value of the resistance is 2 MPa and a standard deviation of 0.2 MPa. The width of pillars is relatively constant, however the height of the pillars (galleries) is variable, we have measures particularly in the upper level. We calculated, for a normal distribution, the average height and the standard deviation are 6.45 m and 0.97 m respectively. From these data we calculated the pillar strength and the maximal and minimal values of R p /R c . according to several hypotheses shown in Table 3.

The arbitrary choice of coefficients (α, β and λ) and the uncertainty or variability of geometric parameters gives an average pillar strength of the order of 1.34 MPa instead of 1.23 MPa (see Table 3) and a standard deviation of between 0.23 MPa and 0.32 MPa.

The minimum and maximum ratios of Rp/Rc, calculated basing on the statistical variation, are respectively 0.8 and 0.54. This approach is used to evaluate the stability of pillars for a homogeneous area in which the dimensions are not perfectly identical, and the rock strength is variable some pillars have very low strengths because of their dimensions and also their resistance. The collapse certainly started at pillars with very low safety factor. In addition, the presence of competent limestone bed in the upper roof had played an essential role on the stability during mining operation [START_REF] Heib | Analysis of the historical collapse of an abandoned underground chalk mine in 1961 in Clamart[END_REF]. (2, 0.2) (7, 0.1) (6.45, 0.97) (1.36, 0.32) (0,8, 0.58) (3)

(2, 0.2) (7, 0.7) (6.45, 0.97) (1.36, 0.23) 0,72, 0,62, It is true that the accumulation of variations on the coefficients (α, β and λ), the uniaxial compression strength and the dimensions of the pillars is a strong assumption and can affect the design results. Nevertheless, it encourages us for the purpose of stability assessment to determine the most accurate way these coefficients (α, β and λ) and know the geometric variability of the pillars in the area in question.

STABILITY OF CHALK MINES (MEUDON)

In the zone of Clamart, there are several underground mines, comparable to Clamart collapsed mine. The Meudon chalk mines were operated on two levels (sometimes more). The mine depth varies between 20 and 50 m. The mining method is room and pillar abandoned. The pillars are generally square in shape and variable dimensions according to the adopted extraction ratio. This extraction rate varies between 30 and 70%. Former mechanical tests were performed to determine the uniaxial compression strength of the chalk, the range of variation is between 1.8 and 2.45 MPa. In addition, 34 cases were analyzed including 15 cases of collapses in nearby areas of similar configurations (from 1843-1961). The objective is to distinguish between stable and collapsed abandoned chalk mines.

We performed an estimation of the safety factor based on the collected data in particular on the geometry of the chalk mines. We considered an average uniaxial compression strength of 2.25 MPa; while the lowest strength is 1.8 MPa, an average height of 6 m pillars (varies between 7 and 5 m). The average strength of the pillars is calculated based on the relationship of Salamon and Munro. Figure 6 shows the histogram of the safety factor for underground chalk mines, the blue columns correspond to stable mines (until now) and red for collapsed chalk mines. The average safety factor for stable mines is equals to 1.2. The average value of collapsed mines safety factor is 0.76, so half the mines regarded to be stable today. More precise analysis is recommended for the mines with an average safety factor is less than 1 and still stable.

If one analyzes collapsed over time, Figure 6 shows that the most recent collapsed mines have generally higher a safety factor than the oldest collapses. This

CONCLUSION

The goal of the study is to discuss the analysis of the stability of pillar. The method of the tributary area is dependent on the choice of parameters and the feedback. We discussed here the pillar strength where two effects should be taken into account: the scale effect (λ < 1) and the pillar geometry. Empirical relationships have been used to design the pillar. We have shown the importance of geometry and W/H ratio. The empirical relationships converge to almost the same results despite the difference between the coefficients of these empirical relationships.

The application of these relationships can help to explain the collapse of Clamart mine. The pillars of Clamart are unstable (SF<1) due to the geometry of pillars and the extraction ratio.

In conclusion, the design or back-analysis of room and pillar stability must take into account the pillar stress ratio and pillar of strength. The stress can be determined using the method of the tributary area and the pillar strength by using the empirical relationships taking into account the scale effect (λ) and the effect of the geometry (α and β). Sophisticated tools are useful to understand the rock mass and mines (pillars) behavior.

Figure 3 .

 3 Figure 3. Evolution of pillar strength to the uniaxial compression strength for different values of W = H (m) W/H = 1 and different empirical equations (see table 1).

  Figure 3. Evolution of pillar strength to the uniaxial compression strength for different values of W = H (m) W/H = 1 and different empirical equations (see table 1).

Figure 4 .

 4 Figure 4. Evolution of pillar strength to the uniaxial compression strength for different values W/H, taking into account standard deviation (s) = 10%, (m = mean).

Figure 5 .

 5 Figure 5. Clamart chalk mine, collapse zone in upper and down levels.

Figure 6 .

 6 Figure 6. Safety factor for 34 underground chalk mines -15 (red and left) are collapse chalk mines. could be explained by two reasons: the first is the reduction of the uniaxial compressive strength over time and the second reason is the limestone bed role The limestone bed has a role important on the stability. It can reduce the load up to its rupture. The role of the stiff bed has been discussed in the case of mass collapses (Al Heib et al., 2014).

Table 1 .

 1 Different parameter values to estimate the pillar strength.

	Reference	Rock type	α	β	λ
	Salamon and Munro (1967)	Coal	0,46	0,66	0,93
	Hedley et Grant (1972)	Quartzite	0,5	0,75	0,57
	Von Kimmelmann et al. (1984)	Sedimentary rock	0,46	0,66	0,69
	Esterhuizen et al. (2011)	Limestone	0,3	0,59	0,92
	W/H effect only		1	1	0,93

Figure 2. Evolution of pillar strength to the uniaxial compression strength for different values of W-to-H ratio and different empirical equations and numerical modelling.

Table 2 .

 2 Different parameter values to estimate the pillar strength.

						SF =
		D W R H W/H Rp/R C σp	Rp/σ P
	Zone	m m % m			MPa
	ZSS1	42 7 57 7	1	0.63	1.95 0.65
	ZSS2	-55 7			1.87 0.67
	ZNS1	5 68 6.8 0.73 0.55	2.63 0.41
	ZSI1	52 7 45 5	1.4	0.85	1.75 0.97
	ZSI2	51 5	1.4		1.96 0.86
	ZSI3	36 5	1.4		1.5	1.26
	D: depth, W: pillar width, H: pillar height, E: extraction ratio.

Table 3 .

 3 Different parameter values to estimate the pillar strength.

	Simul R c	W	H	R p	R p /R c
		(m, s) (m, s) (m, s)	(m, s)	(max, min)
		MPa m	m	MPa
	(1)	(2, 0) (7, 0) (6.45, 0)	(1.34, 0.25) (0.8, 0,54)
	(2)