Characterization of Laguncularia racemosa transcriptome and molecular response to oil pollution

Accepted Manuscript Title: Characterization of Laguncularia racemosa transcriptome and molecular response to oil pollution Authors: Fernanda Alves de Freitas Guedes, Priscilla de Barros Rossetto, F´abia Guimar˜aes, Maur´ıcio Wolf Wilwerth, Jorge Eduardo Santos Paes, Marisa Fabiana Nicol´as, Fernanda Reinert, Raquel Soares Peixoto, M´arcio Alves-Ferreira PII: DOI: Reference:

S0166-445X(18)30483-1 https://doi.org/10.1016/j.aquatox.2018.09.001 AQTOX 5017

To appear in:

Aquatic Toxicology

Received date: Revised date: Accepted date:

25-5-2018 6-9-2018 6-9-2018

Please cite this article as: de Freitas Guedes FA, de Barros Rossetto P, Guimar˜aes F, Wolf Wilwerth M, Santos Paes JE, Nicol´as MF, Reinert F, Soares Peixoto R, Alves-Ferreira M, Characterization of Laguncularia racemosa transcriptome and molecular response to oil pollution, Aquatic Toxicology (2018), https://doi.org/10.1016/j.aquatox.2018.09.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Characterization of Laguncularia racemosa transcriptome and molecular response to oil pollution

Fernanda Alves de Freitas Guedesa, Priscilla de Barros Rossettoa, Fábia Guimarãesa, Maurício Wolf Wilwertha, Jorge Eduardo Santos Paesb, Marisa Fabiana Nicolásc, Fernanda

IP T

Reinertd, Raquel Soares Peixotoe, Márcio Alves-Ferreiraa*

a

SC R

Laboratório de Genética Molecular e Biotecnologia Vegetal, UFRJ - Av. Prof. Rodolpho Paulo

Rocco, s/n – Bloco A – CCS Cidade Universitária – 21941-617 – Rio de Janeiro, RJ, Brasil. b

Centro

de

Pesquisa

e

Desenvolvimento

Leopoldo

Américo

Miguez

de

Mello,

U

PETROBRAS/CENPES - Av. Horácio de Macedo, nº 950 - Cidade Universitária - 21941-915 - Rio

N

de Janeiro, RJ, Brasil.

Laboratório Nacional de Computação Científica – Av. Getúlio Vargas, no333 – Quitandinha - 25651-

A

c

075 – Petrópolis, RJ, Brasil.

M

Laboratório de Ecofisiologia Vegetal, UFRJ – Av. Prof. Rodolpho Paulo Rocco, s/n – Bloco A – CCS

d

ED

Cidade Universitária – 21941-617 – Rio de Janeiro, RJ, Brasil. Laboratório de Ecologia Microbiana Molecular, UFRJ – Av. Prof. Rodolpho Paulo Rocco, s/n – Bloco

e

CC E

PT

K – CCS Cidade Universitária – 21941-617 – Rio de Janeiro, RJ, Brasil.

*corresponding author

Corresponding author telephone and fax: Telephone: +55 21 999358699 Fax: +55 21 25626380

A

Email addresses:

FAFG: [email protected] PBR: [email protected] FG: [email protected] MWW: [email protected] JESP: [email protected]

1

MFN: [email protected] FR: [email protected] RSP: [email protected] MAF: [email protected]



N

U



A quality Laguncularia racemosa transcriptome was assembled and annotated After 12 h, the oil pollution changed the expression profile of L. racemosa leaves and roots Genes related to photosynthesis and protein metabolism were modulated by oil exposure in leaves Genes related to photosynthesis, hypoxia response and ethylene pathway were modulated by oil contamination in roots A suitable pair of reference genes for qPCR data normalization of oil-exposed mangrove is provided

SC R

  

IP T

Highlights

M

A

Abstract:

ED

Mangroves are ecosystems of economic and ecological importance. Laguncularia racemosa (Combretaceae), popularly known as white mangrove, is a species that greatly contributes to the community structure of neotropical and West African mangrove forests. Despite the significance of

PT

these ecosystems, they have been destroyed by oil spills that can cause yellowing of leaves, increased sensitivity to other stresses and death of trees. However, the molecular response of plants

CC E

to oil stress is poorly known. In this work, Illumina reads were de novo assembled into 46,944 transcripts of L. racemosa roots and leaves, including putative isoform variants. In addition to improving the genomic information available for mangroves, the L. racemosa assembled transcriptome allowed us to identify reference genes to normalize quantitative real-time PCR (qPCR)

A

expression data from oil-stressed mangrove plants, which were used in RNASeq validation. The analysis of expression changes induced by the oil exposure revealed 310 and 286 responsive transcripts of leaves and roots, respectively, mainly up-regulated. Enriched GO categories related to chloroplasts and photosynthesis were found among both leaf and root oil-responsive transcripts, while “response to heat” and “response to hypoxia” were exclusively enriched in leaves and roots, respectively. The comparison of L. racemosa 12-h-oil-stressed leaf expression profile to previous Arabidopsis heat-stress studies and co-expression evidence also pointed to similarities between the heat and oil responses, in which the HSP-coding genes seem to play a key role. A subset of the L.

2

racemosa oil-responsive root genes exhibited similar up-regulation profiles to their Arabidopsis homologs involved in hypoxia responses, including the HRA1 and LBD41 TF-coding genes. Genes linked to the ethylene pathway such as those coding for ERF TFs were also modulated during the L. racemosa root response to oil stress. Taken together, these results show that oil contamination affects photosynthesis, protein metabolism, hypoxia response and the ethylene pathway in L. racemosa 12-h-oil-exposed leaves and roots.

IP T

Abbreviations:

SC R

PAH: Polycyclic Aromatic Hydrocarbon

Keywords: hypoxia; Laguncularia racemosa; mangrove; marine oil pollution; protein homeostasis; photosynthesis.PHE: Phenanthrene

U

WSF: Water Soluble Fraction

N

1. Introduction

A

Mangroves comprise tree and shrub species (Tomlinson, 1986) growing in tidal wetlands of tropical and subtropical coastlines (Hensel et al., 2014) that have evolved adaptations to high salinity

M

(halophytes) and flooding (Franke and Schreiber, 2007; Krishnamurthy et al., 2014, 2009; Reef and Lovelock, 2015; Srikanth et al., 2015). Together with species of seven other strict mangroves genera

ED

(Tomlinson, 1986), Laguncularia racemosa (L.) C.F. Gaertn greatly contributes to the community structure of mangrove forests distributed throughout the neotropics and West Africa (Ellison et al., 2010; Hensel et al., 2014; Record et al., 2013; GBIF: https://www.gbif.org). L. racemosa dominance

PT

and higher relative density are associated with areas under greater influence of fresh water, e.g., this species can constitute 28% of the flora species in the Benevente River mangrove forest estuary

CC E

(Espírito Santo State, Brazil) (Petri et al., 2011). Moreover, in Panama, L. racemosa can dominate 89.7% of mangrove areas (Benfield et al., 2005). Despite their economic and ecological importance (Field et al., 1998; Nagelkerken et al.,

2008), the mangals are highly threatened by anthropogenic disturbance (Hensel et al., 2014; Valiela

A

et al., 2001). In addition to coastal development, pollution, mainly oil spills, deeply impacts mangrove ecosystems (Bejarano and Michel, 2016; Beyer et al., 2016; Burns et al., 1993; Duke, 2016; Gabardo et al., 2003; Hester et al., 2016; Langangen et al., 2017; Lewis et al., 2011; Meniconi et al., 2002; Santos et al., 2012). Petroleum-derived oils are mainly composed of polycyclic aromatic hydrocarbons (PAHs), a family of hydrophobic and toxic compounds (Dabestani and Ivanov, 1999). Phytotoxicity studies on species of algae, wetland plants, mangroves and seagrasses demonstrated that oil exposure can widely affect biological processes such as photosynthesis, germination, growth

3

and protein content (Cui et al., 2016; Lewis and Pryor, 2013; Naidoo and Naidoo, 2017; Pezeshki et al., 2000; Proffitt et al., 1995; Zhang et al., 2007). In addition to growth reduction and oxidative stress, studies of Bruguiera gymnorrhiza demonstrated that PAH treatment causes disorganization of root cells and chloroplast distortion (Naidoo and Naidoo, 2016). The bioaccumulation and the impact of oil exposure depend on factors such as the oil type and amount, the duration of contact and the affected part of the plant and species (Ali et al., 2009; Hensel et al., 2014; Ke et al., 2011; Meudec et al., 2006; Naidoo, 2016; Naidoo and Naidoo, 2017; Pezeshki et al., 2000). Plant molecular responses to oil or its byproducts have been poorly investigated. A recent

IP T

global expression analysis by our group showed that Arabidopsis thaliana seedlings exposed to the water-soluble fraction of the marine fuel MF380 (WSF-MF380) displayed a complex response similar to the responses to heat, hypoxia, osmotic and oxidative stresses (Nardeli et al. 2016). This 24-h-

SC R

time-course study showed that heat shock protein (HSP) genes play an important role in the initial Arabidopsis response to oil exposure, as they were strongly up-regulated after 2 h. Additionally, the Arabidopsis transcription factors (TFs) of the APETALA2/ethylene-responsive element binding

U

protein (AP2/EREBP), basic helix-loop-helix (bHLH), MYB-related and C2C2-Co-like families were modulated by WSF-MF380. Another global gene expression analysis showed that A. thaliana

N

exposed to the PAH phenanthrene (PHE) exhibited perturbations in pathways that regulate reactive

A

oxygen species (ROS), responses related to pathogen defense and ethylene signaling (Weisman et

M

al. 2010).

In the present study, a comprehensive L. racemosa transcriptome encompassing leaf and root sequences is reported. The Illumina sequencing allowed us to de novo assemble and annotate

ED

high-quality transcripts that should improve the genomic resources currently available for mangroves. Reference genes selected exclusively from this transcriptomic dataset were accurate

PT

for oil-exposed L. racemosa expression data normalization and were used to validate the RNASeq findings. Further, this study characterized the L. racemosa molecular response to oil stress. Differential expression analysis in leaves and roots revealed that the oil exposure modulates the

CC E

expression of L. racemosa genes related to protein homeostasis, photosynthesis, hypoxia and ethylene response.

2. Materials and methods

A

2.1 Experimental conditions and plant sample harvesting

Cultivation of Laguncularia racemosa (L.) C.F. Gaertn plants and the oil experiment

were conducted as described by Wilwerth et al. (2016). Briefly, L. racemosa propagules and the used sediment were collected from the unpolluted Restinga da Marambaia. Rio de Janeiro, Brazil (23°3´27´´S 43°33´58´´W). Propagules were grown in plastic pots (3 L) containing mangrove sediment and maintained in a greenhouse under 50% natural sun light and temperature ranging from 22.7 to 38.3 °C. Equally irrigated 14-month-old plants were

4

divided into two groups, treated and control plants. Oil-treated plants were exposed to Marine Fuel 380 (MF380) oil through the spillover of 185 mL of oil per plastic pots (3 L) containing one seedling (or equivalently 12 L m -2, as performed by Proffitt and coworkers (1995)). Prior to its application, ten holes were made in the sediment to ensure the oil penetration; holes were also drilled for control plants. The MF380, which is produced by PETROBRAS Distribuidora S.A. (Rio de Janeiro, Brazil), is a PAH-rich oil (49.8 mg.g-1), mainly naphthalenes and phenanthrenes. A detailed MF380 chemical profile is shown in

IP T

Gabardo et al. (2003). In 2000, a pipeline rupture caused the release of this oil into the Guanabara Bay, Rio de Janeiro, Brazil (Gabardo et al., 2003). L. racemosa leaf and root samples were collected 12 h after the treatment. Non-treated samples harvested at the

SC R

same time point served as control. The third totally expanded leaf pair and mixed fragments from three distinct parts of the root system were collected in similar amounts. After harvesting, the root fragments passed through serial washes, followed by three water

nitrogen and stored at -80 °C until RNA extraction.

U

buckets and lastly dried with paper towels. Leaf and root samples were snap-frozen in liquid

N

2.2 Total RNA isolation, cDNA library and Illumina sequencing

A

Total RNA was extracted from 100 mg of ground tissue using the InviTrap Spin Plant Mini Kit

M

(Invitek). RNA concentration and purity were determined using a NanoDropTM Spectrophotometer ND-1000 (Thermo Scientific). RNA integrity was evaluated by 1.2% agarose gel electrophoresis and by an Agilent 2100 Bioanalyzer (Agilent Technologies). For each experimental condition, equal

ED

amounts of high-quality total RNA (RIN > 7) from five different plants were pooled and used to synthesize the four libraries to be sequenced: control and treated leaves and control and treated

PT

roots. The samples were precipitated in 0.1 volume of sodium acetate 3 M (pH 5.2) and 2 volumes of absolute ethanol, in addition to 1 µl of RNase Out (Invitrogen, Life Technologies). The mRNA was purified from 1 µg of total RNA using beads poly-T (Illumina TruSeq Stranded mRNA kit). The mRNA

CC E

was then fragmented using a 0.1 M zinc chloride solution in 0.1 M Tris-HCl (pH 7.0). These shorter fragments served as templates for the first-strand complementary DNA (cDNA) synthesis using random primers, RNase H and SuperScript II Reverse Transcriptase (Invitrogen, Life Technologies). This protocol was followed by end repair, 3' adenylation, adapter ligation and cDNA size selection

A

(≥200 bp) by agarose gel electrophoresis. Finally, the cDNA libraries were sequenced on a HiSeq 2000 Illumina platform, 2x100 bp (paired-end). The mRNA purification, library construction and sequencing were carried out by Fasteris Sequencing Company (Plan-les-Ouates, Switzerland). 2.3 De novo assembly of the L. racemosa transcriptome Raw sequence data were first preprocessed to remove adapters and poly A/T tails and filtered by base quality scores (average Q ≥ 30). The resulting filtered paired-end reads were de novo assembled into contigs using a combined approach of Velvet version 1.2.7 (Zerbino and Birney,

5

2008) and Oases version 0.2.8 (Schulz et al., 2012). The preliminary assembly produced by Velvet was inputted to Oases. Additionally, several values of kmer (81 to 91 bp) were evaluated to optimize the assembly process. To check the assembly quality and estimate reads counts for differential expression analysis, the BWA version 0.5.9 (Li and Durbin, 2009) was used to align reads against the L. racemosa transcriptome assemblies. Sequencing reads used in this article were deposited in the NCBI Short Read Archive (SRA) database under the Bioproject accession number PRJNA355123, and this Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession number GGMP00000000. The version described in this

IP T

paper is the first version, GGMP01000000. 2.4 Bioinformatics analysis

SC R

The L. racemosa sequences were annotated by comparison to different databases. First, they were searched against the non-redundant (nr) database using blastx with e-value cutoff 10-5. Gene Ontology (GO) terms were assigned to the L. racemosa sequences using Blast2GO (Conesa et al., 2005), followed by functional classification with the WEGO online tool (Ye et al., 2006).

U

Blast2GO was also used to perform the GO enrichment analysis by applying a Fisher's exact test

N

with FDR correction and q-value cutoff 0.05. The GO enriched categories were further reduced to the most specific ones by removing the general terms (terms from the highest levels of GO ontology).

A

KEGG pathway analysis was made using KOBAS standalone version 2.0 (Xie et al., 2011) by

M

searching L. racemosa sequences against A. thaliana background with e-value cutoff 10-5. The L. racemosa translated protein sequences, predicted by ESTScan version 3.03 (Iseli et al., 1999), were also annotated against the Eukaryotic Orthologous Groups (KOG) and Pfam databases by WebMGA

ED

server (Wu et al., 2011) with e-value cutoff 10-5. L. racemosa genes regulated by oil exposure were identified using the Cuffdiff package from

PT

Cufflinks version 2.0.2 (Trapnell et al., 2010). Quality reads of each library were first aligned to the assembled 87-kmer transcriptome with BWA version 0.5.9 (Li and Durbin, 2009). Fragments per

CC E

kilobase of transcript per million mapped reads (FPKM) normalization and differential expression analysis were carried out by comparing the oil-treated and control libraries of each tissue. Differentially expressed transcripts were identified using default Cuffdiff parameter values and the FDR-adjusted p-value < 0.05; no additional filters were used to determine the differential expressed

A

transcripts. Despite the Cuffdiff recommendation to set “blind” option for dispersion-method when the experiment design counts on one replicate per condition, the default option (pooled) was used here to avoid the construction of a very conservative dispersion model that could result in no significant calls. The log2FC (base-2 logarithm of the fold change) represents the ratio between the estimated expression values from treated and control libraries. For the oil-responsive genes, the ATcode assignment was done by transferring the best-hit ID obtained by blastp comparison (e-value 10-5) against the TAIR10 database. The co-expression networks were constructed by transferring high-confidence (score > 0.7) evidence from the STRING 10.5 database (Szklarczyk et al., 2015) to

6

L. racemosa responsive genes through their assigned ATcodes. Assigned ATcodes were also used to compare the log2FC values of oil-responsive genes with other, available abiotic stress expression data using Genevestigator (Hruz et al., 2008). Genevestigator hierarchical clustering tool was used to cluster these expression data by gene and condition, applying Euclidean distance. 2.5 Reference genes expression stability analysis and RNASeq validation by quantitative realtime PCR Seven assembled loci (Table S1) selected as reference gene candidates had their

IP T

expression stability evaluated in 12-h-oil-exposed L. racemosa leaves. As tubulin is a commonly used internal control, one putative L. racemosa tubulin (Locus_25136, here called TUB) was selected. Three other candidates (Locus_359, Locus_12010, Locus_5702) came from the list of

SC R

Arabidopsis reference genes proposed by Czechowski and coworkers (2005). The remaining three candidates (Locus_3266, Locus_9942, Locus_1622) were selected by exclusively searching the set of assembled transcripts, using a combination of aligned read counts and function. The expression stability of these candidates was evaluated using geNorm (Vandesompele et al., 2002) and

U

NormFinder software (Andersen et al., 2004).

N

For RNASeq validation by quantitative real-time PCR (qPCR), 12 loci were selected by

A

applying the following criteria: function, |log2FC| ≥ 1.3, number of aligned reads ≥ 100 for treated and control libraries and sequence length > 500 bp. Primer pairs were designed using Primer3Plus

M

(Untergasser et al., 2012) with the following parameters: length of 18-24 nucleotides, melting temperature of 59-61 °C, product size ranging from 70 to 150 bp and GC content of 50-60%. The

ED

primer specificity was checked against the L. racemosa transcriptome using Primer-Blast tool (http://www.ncbi.nlm.nih.gov/tools/primer-blast/).

PT

RNA extraction followed the protocol previously described by Wilwerth et al. (2016) with minor modifications. Briefly, leaf RNA samples from three different plants were pooled and treated with TURBO DNase (Ambion, Thermo Fisher) to compose each biological replicate. To improve the

CC E

RNA quality, cleanup steps were added. First, 200 μl of phenol:chloroform:isoamyl alcohol (25:24:1) was added to RNA samples (100 μl). The resulting mixture was centrifuged at 18,000 g for 10 min at -4 °C. RNA was precipitated from the aqueous phase by adding 1/10 volume of sodium acetate 3 M and 1 volume of cold isopropanol, mixing and incubating at -20 °C for 1 h. After a 20 min

A

centrifugation, the supernatant solution was removed, and the RNA pellets were washed with 400 μl cold ethanol 70%. The samples were centrifuged for 5 min at 18,000 g at -4 °C, and the ethanol solution was removed. Finally, the RNA pellets were air-dried for 10 min and resuspended by adding 10 μl nuclease- and DNA-free water. The total RNA concentration and purity were evaluated by a NanoDropTM spectrophotometer ND-1000 (Thermo Scientific) and qubit RNA (Thermo Fisher Scientific). The first-strand cDNA was synthesized by adding 2 μg of clean DNase-treated total RNA (control and treated leaves), 50 μM poly-dT (oligo dT20) and 10 μM each deoxyribonucleoside 5'triphosphate (dNTP). This mixture was incubated at 65 °C for 5 min and briefly chilled on ice. Then,

7

20 mM dithiothreitol (DTT), 1X first-strand buffer and 200 U of SuperScript III Reverse Transcriptase (Invitrogen, Life Technologies) were added to the mixture. The total volume (20 μl) was incubated at 50 °C for 1 h following the manufacturer’s instructions. The reverse transcriptase was inactivated by incubating the mixture at 70 °C for 15 min. The qPCR reactions were carried out in an optical 96-well plate in a 7500 Fast Real-Time PCR detection system (Applied Biosystems) following the manufacturer’s instructions. The amplification reactions were performed in a 20 μl final volume containing 10 μl cDNA (1:50), SYBER® Green 1X, 0.4 μM of each primer (forward and reverse), 0.025 mM dNTP, PCR buffer 1X,

IP T

3 mM MgCl2, 0.25 U Platinum Taq Polymerase and 0.4 μl ROX reference dye. The reaction mixtures were incubated at 94 °C for 5 min, followed by 40 cycles of 94 °C for 15 s, 60 °C for 10 s, 72 °C for 15 s and 60 °C for 35 s, and a final denaturation curve from 30 °C to 100 °C for 1 min. The melting

SC R

curve analysis of the amplification products confirmed that the primers amplified only a single product with the expected size (data not shown). Three independent biological replicates of each experimental condition were evaluated in technical triplicate. Both the efficiency and the

U

quantification cycle (Cq) values generated by the qPCR were estimated using the Miner online tool (Zhao and Fernald, 2009). Cq values were converted using qBase software (Hellemans et al., 2007)

N

into non-normalized relative quantities. The relative expression ratio between oil-treated and control

A

samples was calculated using the Relative Expression Software Tool (REST) (Pfaffl et al. 2002). The REST statistical test applied on Cq values was the pair-wise fixed reallocation randomization

M

test.

ED

3. Results

3.1 Sequencing and assembly of transcriptome

PT

Leaf and root samples harvested from control and 12-h-oil-exposed L. racemosa plants were used to construct the cDNA libraries. Illumina paired-end sequencing yielded between 36 and 51 million quality reads for L. racemosa libraries (Table S2). These reads were de novo assembled into

CC E

contigs using a combination of Velvet and Oases assemblers and different kmer values (81, 85, 87, 89, 91) (Table S3). For each tested kmer value, the assembly produced by Velvet served as input to Oases. To evaluate the quality of the assemblies, the reads were aligned to each assembled transcriptome using BWA. Although contig number and mean length were also considered, the

A

percentage of aligned reads was the main criterion to choose the kmer value. Given that the Oases final 87-kmer assembly showed the highest mean percentage of aligned reads (84.89%) as well as for leaf libraries (Table S3), this value was the chosen value. For kmer=87, the first assembly step carried out with Velvet resulted in 59,473 transcripts with an average size of 670 bp and N50=957 bp (Table S3). Finishing the assembly with Oases resulted in a reduction of contig number accompanied by an increase in the mean length and N50 (Table S3). The final L. racemosa assembly contained 46,944 transcripts corresponding to 35,795 putative loci and a total length of 55.3 Mb (Table 1). The mean contig length was 1,179 bp, the longest transcript had 11,799 bp, and 45.92%

8

(21,556) of transcripts had at least 1 kb. A reduction in the number of sequences shorter than 1,000 bp was observed together with an increase in the number of longer transcripts in the final Oases assembly (Figure S1). A comparison of the L. racemosa transcriptome assembly with RNASeq studies recently conducted with other mangrove and semi-mangrove species (Table 2) showed that the present study had both deep sequencing and longer assembled transcripts. 3.2 Functional annotation of L. racemosa root and leaf transcripts The L. racemosa transcriptome was annotated against different biological databases to get

IP T

the broadest possible functional information. Annotation based on sequence similarity revealed that 73.88% (34,683) of L. racemosa transcripts matched with nr sequences, whose most represented top-hit species were Theobroma cacao (8,156 sequences), Vitis vinifera (7,323 sequences), Prunus

SC R

persica (3,939 sequences), Ricinus communis (3,890 sequences) and Populus trichocarpa (3,372 sequences) (Figure S2). For sequences with blast hits, 23,021 (66.38%) also had at least one assigned GO term. The transcriptome GO annotation overview obtained with the WEGO tool is shown in Figure S3. For further functional classification, L. racemosa sequences were also searched

U

against the KOG database. Almost 42% (19,628) of the assembled transcripts were assigned to KOG classes (Figure 1). The most prevalent KOG class was “signal transduction mechanisms”

N

(3,088; 15.73%), followed by “general function prediction only” (2,893; 14.74%), “posttranslational

A

modification, protein turnover, chaperones” (2,140; 10.90%) and “transcription” (1,304; 6.64%).

M

Using the KOBAS tool to perform the KEGG pathway annotation, 71.12% (33,387) of L. racemosa transcripts showed similarity to A. thaliana sequences, of which 6,021 (18.03%) were assigned to 124 different pathways. Among them, “metabolic pathways”, “biosynthesis of secondary

ED

metabolites”, “ribosome”, “plant hormone signal transduction” and “glycolysis/gluconeogenesis” were the top five most represented pathways. Moreover, the KEGG pathway enrichment analysis

spliceosome.

PT

showed three statistically significant pathways (corrected p-value < 0.01) (Table 3), including the

CC E

3.3 Identification and GO enrichment of L. racemosa oil-responsive transcripts To identify the L. racemosa transcripts whose expression levels significantly changed in

response to oil, differential expression analysis was carried out with Cuffdiff. In leaves, 310 regulated transcripts (Table S4) were found, most of which (237 transcripts, 76.45%) were up-regulated.

A

Similarly, root oil-responsive transcripts (286) were mainly up-regulated (251 transcripts, 87.76%). Of all the L. racemosa modulated transcripts (586), only 10 (0.02%) were shared between leaves and roots (Table S5). The GO analysis of oil-responsive transcripts revealed 118 enriched categories (Table S6), which could be further reduced to the most specific ones. Among these specific enriched categories (leaves – 23; roots – 30) (Figure 2), the overrepresented “metal ion binding” (57 sequences) stands out in roots. For biological process terms, functional categories related to response to abiotic factors, photosynthesis and rRNA transcription were enriched both in leaves and roots (Figure 2; Table S6). Notably, “response to heat”, “heat shock protein binding”, “unfolded

9

protein binding” and “protein folding” were exclusively enriched among the leaf modulated transcripts, while “response to hypoxia” was exclusive to roots. 3.4 Comparative analysis of the L. racemosa expression profile changes induced by MF380 oil Following the GO enrichment results, the L. racemosa differentially expressed loci (leaves – 235; roots – 238) were compared to expression data of Genevestigator stress experiments collection through their assigned ATcodes. We identified 124 and 171 putative A. thaliana homologs

IP T

corresponding to 143 (60.85%) and 187 (78.57%), respectively, of L. racemosa leaf and root oilresponsive loci. In leaves, the comparison of L. racemosa genes regulated by the MF380 oil with previous Arabidopsis heat studies (Figure 3; Figure S4; Table S8) revealed an overlapping group of

SC R

highly induced HSP-coding genes. Another group of L. racemosa oil-responsive genes whose Arabidopsis homologs were induced by heat comprised the ribosomal protein RPSL2, HSFB2A (heat shock transcription factor B2A) and the PDX1.2 (pyridoxine biosynthesis 1.2) (Figure 3). Additionally, the set of oil-responsive root genes was paralleled by hypoxia-responsive genes (Figure 4; Figure

U

S5; Table S8). The outcome of this comparative analysis revealed a group of L. racemosa up-

N

regulated genes whose Arabidopsis homologs are involved in the response to hypoxia/anoxia conditions, such as ADH1 (alcohol dehydrogenase 1), PCO2 (plant cysteine oxidase 2), PDC1

A

(pyruvate decarboxylase 1) and the TFs HRA1 (hypoxia response attenuator 1) and LBD41 (LOB

M

domain-containing protein 41) (Figure 4). Three L. racemosa ethylene-related genes whose expression was induced by the oil exposure (Table 4) had their corresponding Arabidopsis homologs

ED

also induced by hypoxia, including the TF ERF71 (ethylene response factor 71). Further, the L. racemosa oil-responsive loci were specifically compared to the Arabidopsis oil- and PAH-modulated genes identified by previous studies (Nardeli et al., 2016; Weisman et al.,

PT

2010). In leaves, 16 of 19 L. racemosa modulated genes whose putative homologs were also modulated in Arabidopsis seedlings treated with WSF-MF380 (Figure 5) belonged to the HSP class

CC E

(Table S7). On the other hand, the comparison between leaf genes regulated by MF380 and PHE as well as those regulated by the three experiments revealed more functionally diversified sequences (Figure 5; Table S7). Among the leaf genes commonly regulated by the three experiments (3), one transcript (Locus_7764) that encodes for the putative homolog of LHCA1 (photosystem I light

A

harvesting complex gene 1) was also regulated in L. racemosa roots (Table S7). Surprisingly, 22 root genes whose putative homologs were also regulated by PHE in Arabidopsis were functionally related to chloroplasts and/or photosynthesis (Table S7). Additionally, six TFs modulated by the other treatments were up-regulated by oil exposure in L. racemosa roots (Table S7), including ERF71, ERF72 and HRA1. Only one putative TF, AGAMOUS-like 20, was down-regulated by oil exposure in L. racemosa roots. 3.5 The co-expression network underlying L. racemosa response to oil

10

To obtain a more systemic view of its oil response, the co-expression information was incorporated into the analysis of L. racemosa regulated genes through the construction of networks using their assigned ATcodes and STRING high-confidence co-expression evidence. The L. racemosa oil-responsive leaf network (Figure 6A) revealed two major clusters of co-expressed genes. One cluster (Figure 6A, cluster 1) was formed mainly by the up-regulated HSP-coding genes, including those shared with the Arabidopsis WSF-MF380 experiment (Table S7). Another cluster of induced co-expressed genes coding for ribosomal proteins (RPs) (Figure 6A, cluster 2) emerged from the L. racemosa leaf network, embracing the heat-modulated RPSL2 (Figure 4) and the PHE-

IP T

modulated RPL10C (Table S7). A pair of photosynthesis-related genes (Figure 6A, cluster 4) also appeared in this co-expression network. Surprisingly, a large cluster (Figure 6B, cluster 5) of 24 oilinduced genes linked to photosynthesis were co-expressed in L. racemosa roots. Another cluster

SC R

(Figure 6B, cluster 6) of the co-expression root network contained four up-regulated genes related to response to hypoxia, including the TFs HRA1 and LBD41. All these hypoxia-related genes were also regulated in Arabidopsis by WSF-MF380 (Table S7). A third remarkable oil-responsive root cluster (Figure 6B, cluster 7) combined four induced TF-coding genes, ERF1, ERF105, ZAT10, and

U

ZAT12. It is interesting to note that as opposed to that observed for leaves, the two co-expressed

N

root HSP-coding genes were down-regulated by oil exposure (Figure 6B, cluster 8).

A

3.6 Evaluation of reference genes and RNASeq validation

M

Seven L. racemosa candidate reference genes had their expression stability evaluated by qPCR in 12-h-oil-exposed leaf samples (Table S1). The selected genes belonged to different functional classes to reduce the chances of co-regulation among them. The Cq values obtained for

ED

these candidate genes ranged from 20.77 up to 28.63 (Table S9), and Locus_25136 (TUB), Locus_9942 (Ref4) and Locus_1622 (Ref6) had the highest expression. The NormFinder and

PT

geNorm tools, which are based on distinct statistical algorithms, were combined to assess the expression stability of L. racemosa candidate reference genes. The NormFinder identified Ref1/Ref6 as the most stable pair of genes (Table S10). On the other hand, the stability rank generated by the

CC E

geNorm average stability values (M) indicated to P91/P97 as the most stable pair, followed by Ref1 and Ref6 genes. All the M values calculated by geNorm were smaller than the cutoff (M<1.5) suggested by Vandesompele et al. (2002) (Figure 7A). Furthermore, the geNorm pairwise variation (V) assessed the optimal gene number for accurate normalization. Adopting the cutoff V<0.15, the

A

observed value V2/3 (V=0.097) indicated that two reference genes were enough for L. racemosa expression data normalization (Figure 7B). Combining the stability results from NormFinder and geNorm, Ref1 and Ref6 were chosen as reference genes. Notably, instead of being classical and/or well-studied internal genes, the most suitable reference genes were selected exclusively on the basis of the L. racemosa assembled transcriptome. Ref1 and Ref6 were then used to normalize the expression of target sequences (Table S1) from qPCR data for RNASeq validation. Primer efficiencies and mean Cq values for the selected target genes are shown in Table S9. The statistical

11

significance analysis carried out with REST software confirmed the RNASeq-estimated expression changes for 10 of 12 (83.33%) L. racemosa tested genes (Figure 8). The selected up-regulated genes were all validated, including those that appeared in the heat-induced group (Figure 3) and coexpression network (Figure 6A). Although most of the oil-responsive genes were induced, three down-regulated genes were also included in the set of target genes for validation, one of which was validated by qPCR. 4. Discussion

IP T

4.1 Assembly of quality L. racemosa sequences Despite the great ecological and economic relevance of mangrove ecosystems, mangrove plant sequences are poorly represented in biological databases. Searching in the NCBI database

SC R

returned a small number of L. racemosa nucleotide (28) and protein (12) sequences. Even in the specific Mangrove Transcriptome Database (MTDB; http://mangrove.illinois.edu) (Dassanayake et al., 2010), gene and protein sequences from L. racemosa, which ranks among the dominant species in South American mangrove areas, are scarce. In the present study, the final assembled L.

U

racemosa transcriptome contained 46,944 sequences. Compared to previous RNASeq studies

N

carried out with other semi-mangrove and mangrove species (Table 2), in this study, the L. racemosa

A

transcriptome contained sequences with a higher mean length, probably representing more fulllength sequences. Factors such as sequencing depth, the transcriptome non-uniform coverage

M

(intrinsic variation of transcript expression levels) and the user-defined kmer (Surget-Groba and Montoya-Burgos, 2010) can affect the de novo transcriptome assembly. Since the De-Bruijn-graph

ED

assemblers such as Oases (Schulz et al., 2012) are very sensitive to the kmer parameter, L. racemosa transcriptome assemblies generated with different kmer values (Table S3) were first evaluated, and the intermediate value (87) was chosen, representing a compromise between

PT

sensitivity (lower values) and specificity (higher values). Therefore, the combination of the attained sequencing depth and the selected kmer value for the transcriptome assembly might have helped

CC E

us obtain quality L. racemosa sequences. The transcriptome assembled with Velvet (Zerbino and Birney, 2008) and Oases (Schulz et

al., 2012) included variants, such as putative isoforms, expressed in leaves and roots of L. racemosa. Isoforms originate from the alternative splicing process carried out by the spliceosome and provide

A

variation through the creation of multiple transcripts from a single gene (Barbazuk et al., 2008). Alternative splicing may act as an important regulatory mechanism in the plant response to environmental stimuli (Mastrangelo et al., 2012; Staiger and Brown, 2013). The spliceosome pathway was significantly enriched in the L. racemosa transcriptome (Table 3). This pathway was also among the most represented KEGG pathways for the semi-mangrove Millettia pinnata (Huang et al., 2012) and for the mangrove species Avicennia marina (Huang et al., 2014). Moreover, a comparative analysis conducted between L. racemosa and M. pinnata revealed the GO category “mRNA splicing, via spliceosome” was enriched in the set of homolog sequences (data not shown).

12

Taken together, these findings suggest that the splicing might be an active RNA editing mechanism in the semi-mangrove and mangrove species such as L. racemosa. Hence, the L. racemosa assembled transcriptome, including its putative variants, may help enrich databases and improve the current knowledge about mangrove species. 4.2 Supplying reference genes for accurate normalization of qPCR expression data from oilstressed plants and RNASeq validation Even though qPCR is a robust and useful tool for studying gene expression profiles (Gachon

IP T

et al. 2004), finding stable reference genes for accurately normalizing qPCR data is a critical step. To search for appropriate reference genes to normalize the expression data from L. racemosa oilstressed plants, we evaluated the stability of seven selected leaf loci. Despite being a frequently

SC R

used reference gene (Czechowski et al., 2005), the stability values found for the tested L. racemosa tubulin gene (Figure 7; Table S10) indicated that it was not the most appropriate gene to normalize L. racemosa expression data, as already reported for other plant species (Artico et al., 2010; Silveira et al., 2009). The Arabidopsis list provided by Czechowski and colleagues (2005) was also used to

U

select three L. racemosa candidate reference genes, two of which (P91 and P97) were considered

N

the most stable according to the geNorm algorithm (Figure 7). GeNorm analysis also indicated that two was the ideal number of reference genes to normalize the L. racemosa expression dataset

A

(Figure 7). Finally, three other candidates were selected based only on the RNASeq expression

M

estimates, two of which (Ref1 and Ref6) were picked by NormFinder as the best pair (Table S10). In contrast to the geNorm algorithm (Vandesompele et al., 2002), the NormFinder algorithm (Andersen et al., 2004) assesses the expression variation in a broader way, taking both intra- and

ED

intergroup variations into account for normalization factor calculations. Moreover, its suggested best pair was well positioned in the geNorm stability rank. Thus, the Ref1/Ref6 pair was selected for L.

PT

racemosa expression data normalization. To the best of our knowledge, this is the first study to provide reference genes for normalization of qPCR expression data from mangroves contaminated

CC E

with oil.

To assess the reliability of the L. racemosa selected reference genes and to validate the

RNASeq predicted oil-induced expression changes, the expression of 12 leaf differentially expressed loci was checked by qPCR (Table S9). Only leaf samples were chosen to perform the qPCR

A

validation. The exclusion of root samples was due to the poor quality of RNA after extraction. As demonstrated by Wilwerth and colleagues (2016), oil-impregnated L. racemosa root samples must be washed several times before RNA extraction. Additionally, PAHs are mostly retained in root cells (Naidoo and Naidoo, 2016; Tao et al., 2006; Wang et al., 2015), probably necessitating extra steps of RNA cleanup. The qPCR analysis using Ref1 and Ref6 as reference genes confirmed the significant regulation predicted by Cuffdiff for 10 of 12 (83.33%) tested genes (Figure 8), confirming the reliability of the chosen reference genes and the RNASeq differential expression analysis. 4.3 Providing a functional profile of the L. racemosa transcriptome

13

As a non-model mangrove, functional information is scarce even for L. racemosa’s closely related species. Nevertheless, 74% (34,683) of L. racemosa (Combretaceae) transcripts had similarity with nr sequences. In addition to the comprehensive and widely accepted annotations currently available for the top five blast-hit species of L. racemosa sequences (Figure S2), an overlap was found with annotation sources of the mangroves A. marina (Acanthaceae) (Huang et al., 2014), Rhizophora mangle (Rhizophoraceae) and Heritiera littoralis (Malvaceae) (Dassanayake et al., 2009) and for the semi-mangrove M. pinnata (Fabaceae) (Huang et al., 2012). Additionally, since four of the five most represented top-hit species were trees, the L. racemosa blast-hit results might

IP T

be related to its tree habitat. Taken together, these annotation results indicate that the L. racemosa transcriptome generated in the present study was confidently assembled and annotated.

The GO annotation of L. racemosa transcriptome (23,021 sequences; 49.04%) allowed us to

SC R

construct its functional profile (Figure S3). Similar to the observed in other semi-mangrove and mangrove transcriptome studies (Chen et al., 2011; Huang et al., 2014, 2012; Yang et al., 2015), “metabolic process”, “cellular process”, “biological regulation” and “response to stimuli” contained

U

the greatest number of L. racemosa transcripts among the biological process categories. L. racemosa protein sequences obtained by translating predicted coding regions (35,961) were

N

searched against KOG database (Figure 1). Like the extreme-environment-adapted species

A

Firmiana danxiaensis (Chen et al., 2015), the three KOG categories with the greatest number of L. racemosa transcripts were “signal transduction mechanisms”, “general function prediction only” and

M

“posttranslational modification, protein turnover, chaperones”. The major KOG clusters found for the A. marina mangrove (Huang et al., 2014) and for the recretohalophyte Reaumuria trigyna (Dang et

ED

al., 2013) transcriptomes also included these top-three categories. These transcriptome annotation results suggest that mechanisms of signal transduction and protein modification might be involved

PT

in L. racemosa adaptation to the mangrove environment. 4.4 Oil pollution affects photosynthesis and protein homeostasis in L. racemosa leaves

CC E

Mangrove plants are often exposed to pollution caused by oil spills (Beyer et al., 2016; Duke, 2016; Gabardo et al., 2003; Hester et al., 2016; Lewis and Pryor, 2013; Proffitt et al., 1995; Santos et al., 2012). In L. racemosa leaves, the PAH-rich MF380 oil mainly promoted the induction of genes expression 12 h after the exposure (Table S4). The GO analysis showed that the “photosynthesis,

A

light harvesting” category was enriched in the set of oil-responsive leaf transcripts (Figure 2). Among them, three L. racemosa up-regulated genes coding for light harvesting complex proteins (LHCA1, LHCB1.3 and LHCB4.2) were identified. The co-expressed LHCA1 and LHCB4.2 genes (Figure 6, cluster 4) also had putative homologs that were modulated by PHE in Arabidopsis (Table S7). Effects of PAHs on photosynthesis have been reported in several plant species (Ahammed et al., 2012a, 2012b; Jin et al., 2017; Kummerová et al., 2006; Kummerová and Váňová, 2007; Li et al., 2008; Tomar and Jajoo, 2014, 2013), as well as the modulation of photosynthesis-related genes (Weisman et al., 2010). Our group recently observed an increase in the photosynthetic performance indexes

14

PITotal and PIABS in L. racemosa leaves after 12 h of oil exposure followed by the reduction of these indexes from the 7th day onwards, revealing the negative effects of oil contamination (Reinert et al., 2016). Taken together, our results show that L. racemosa photosynthesis is affected on different levels by oil exposure. We hypothesize that the early induction of L. racemosa genes coding for LHC proteins, which are components of the antenna system of the photosynthetic apparatus (Jansson, 1999), is not a sustainable response to oil exposure to increase the efficiency of light energy absorption and/or transfer that would be required to overcome the metabolic requirements of oil stress acclimation.

IP T

L. racemosa protein homeostasis also seemed to be altered by oil exposure, as indicated by the modulation of genes coding for RPs, HSPs and chaperones. Weisman and coworkers (2010) have already reported the regulation of genes involved in protein biosynthesis by PAH treatment.

SC R

Among the L. racemosa responsive co-expressed RPs (Figure 6, cluster 2), RPSL2 also grouped with the heat-induced genes (Figure 3), and RPL10C was induced by PHE treatment (Table S7). In addition to being ribosome constituents, RPs can act off these cellular components to alert the cell

U

to stress or defects in ribosome synthesis itself (Carvalho et al., 2008; Ferreyra et al., 2013, 2010; Moin et al., 2016; Nagaraj et al., 2016; Sormani et al., 2011; Warner and McIntosh, 2009). It is still

N

noteworthy that, in contrast to RPSL2, the remaining co-expressed oil-responsive RP-coding genes

A

(Figure 6) did not show the same remarkable expression induction profile in heat experiments (Figure S4). These results suggest that RPs might play a role in the L. racemosa leaf response to oil stress

M

in a way that differs from the heat stress response.

As indicated by the enriched categories “protein folding”, “unfolded protein binding” and “heat

ED

shock protein binding” (Figure 2) and by the observed similarities with the heat expression profile (Figure 3), the protein folding mediated by HSPs might be affected by oil pollution in L. racemosa

PT

leaves. HSPs and chaperones act in protein folding, assembly, translocation, degradation, and stabilization and can assist in their refolding under stress conditions (Jacob et al., 2017; Wang et al., 2004). An A. thaliana transcriptional profiling study showed that HSP genes exhibit expression

CC E

changes in a wide range of biotic and abiotic stress types, revealing an extensive overlap between heat and non-thermal stressful conditions (Swindell et al., 2007). Nardeli and colleagues (2016) additionally demonstrated that the treatment of A. thaliana seedlings with WSF-MF380 led to the regulation of HSP-coding gene expression, whose putative homologs were also modulated and co-

A

expressed in L. racemosa leaves (Figure 6, cluster 1). These results suggest that the HSPs play a key role in the plant response to oil stress, probably acting in a coordinated way, as occurs in other stressful conditions. The identified L. racemosa HSFB2A TF (Locus_18599), whose Arabidopsis homolog exhibited a heat-induced expression profile, might act as a transcriptional regulator of L. racemosa HSPs in response to oil stress. Taken together, our results strongly suggest that PAH compounds can widely affect L. racemosa leaf protein metabolism and that, to address the negative effects of oil pollution, L.

15

racemosa plants seem to resort to molecular mechanisms that both control the protein synthesis and folding to neutralize protein damage. 4.5 Chloroplast-located processes, hypoxia response and ethylene pathway are modulated by oil in L. racemosa roots Differences between the molecular response of L. racemosa roots and leaves to oil stress were observed. In addition to the expected differences in the transcriptome of these organs, the oilexposed roots tend to accumulate higher amount of PAHs than aerial parts (Kang et al., 2010;

IP T

Naidoo and Naidoo, 2016; Tao et al., 2006; Wang et al., 2015). Additionally, high-molecular-weight (four or more rings) PAHs generally are retained in roots, while the low-molecular-weight (two or three rings) PAHs can also be translocated to the leaves (Naidoo and Naidoo, 2016). Another factor

SC R

that could explain the differences observed between leaves and roots is the oil accumulation on the sediment surface, which might promote, by exposure to sunlight, changes in the PAH components. PAHs can absorb ultraviolet light (290-400 nm), and photoinduced processes can generate photoproducts with increased toxicity (Dabestani and Ivanov, 1999; Huang et al., 1997, 1996;

U

Kummerová and Kmentová, 2004; Tomar and Jajoo, 2015). For a detailed view of the experiment,

N

see Figure 2 of Wilwerth et al. (2016).

A

Chloroplast-located processes were modulated by the PAH-rich MF380 oil in L. racemosa roots in a more extensive way than found for leaves, as revealed by the GO enrichment analysis

M

(Figure 2) and the co-expression evidence (Figure 6, cluster 5). Consistent with this, a L. racemosa gene (Locus_3066) that encodes for a TF of GATA family was induced in 12-h-oil-exposed roots

ED

(Table S4). Its putative Arabidopsis homolog (AT4G26150) is a master regulator of chloroplast biogenesis (Chiang et al., 2012). Some mangrove species count on above-ground adapted roots that can photosynthesize, supplying adequate oxygen levels to underground roots and/or fixing

PT

exogenous CO2 (Aschan and Pfanz, 2003; Dromgoole, 1988; Srikanth et al., 2015). The photosynthesis by aerial roots can contribute to the carbon economy of the entire plant (Flores et al.,

CC E

1993). In L. racemosa, the photosynthesis-related co-expression cluster with 24 loci putatively encoding proteins of photosystems I and II, representing 22 Arabidopsis homologs, were induced by oil exposure (Figure 6, cluster 5). Photosystems I and II are functional units that, together with chlorophyll and accessory pigments, absorb light energy and transfer it, ultimately resulting in ATP

A

synthesis (Taiz and Zeiger, 2002). Although the putative Arabidopsis homologs were previously shown to be PAH-responsive (Table S7) (Weisman et al., 2010), the results obtained in the present study show for the first time that oil exposure can regulate the expression of photosynthesis-related genes specifically in roots. Notably, their corresponding Arabidopsis genes are also modulated during the response to a 2-h-long hypoxia stress (Figure S5). Taken together, our results suggest that, like other mangroves, L. racemosa roots can adapt to photosynthesize, and they support the hypothesis that the activation of this process in response to oil exposure might be due to the lowoxygen condition caused by the oil root coating.

16

Consistent with the finding that the “response to hypoxia” category was exclusively enriched in roots, the comparative analysis against Arabidopsis hypoxia expression data (Figure 4) showed a group of L. racemosa up-regulated hypoxia markers (Giuntoli et al., 2017; Licausi et al., 2011). These genes included PCO2, PDC1, ADH1 and the HRA1 and LBD41 TFs, which also formed coexpression clusters (Figure 6, cluster 6 and 10). We hypothesize that the percolation of the viscous MF380 oil through the sediment followed by the root impregnation reduced the oxygen availability to L. racemosa roots, subsequently inducing the hypoxia response. McKee (1996) found that L. racemosa plants may exhibit greater vulnerability to low-oxygen tension in roots than other

IP T

mangroves, such as R. mangle. In another cluster of root co-expressed genes (Figure 6, cluster 7), two up-regulated genes encoding ERF TFs (ERF1 and ERF105) were found. Genes encoding members of the AP2/EREBP TF family and the ethylene pathway were also up-regulated in 12-h-

SC R

oil-exposed L. racemosa roots (Table 4). The study conducted by Weisman and colleagues (2010) with ethylene mutants found that this hormone appears to be involved in the Arabidopsis response to PHE. Moreover, the AP2/EREBP TF family was one of the most represented TF families found by Nardeli and colleagues (2016) in the set of genes responsive to WSF-MF380. These findings

U

corroborate the hypothesis that the ethylene pathway is activated and that AP2/EREBP family TFs

N

may act as transcriptional regulators during the plant response to oil stress. Notably, genes in these TF-containing root clusters are connected in a topology that resembles common transcriptional

A

regulatory network motifs such as bifan and feed-forward loops (FFLs) (Alon, 2007; Babu et al.,

M

2004; Burda et al., 2011; Milo et al., 2002; Wang and Purisima, 2005). The FFL is formed by three nodes, where two regulatory genes control the expression of a target, while in the four-node bifan

ED

motif the products of two regulatory genes control the expression of two targets. The results obtained here suggest possible targets for the responsive L. racemosa root TFs that might act in an intricate

5. Conclusion

PT

way during the plant response to oil pollution.

The characterization of the L. racemosa molecular response to oil stress revealed that the

CC E

expression of photosynthesis-related genes was modulated in both leaves and roots exposed to this pollutant. On the other hand, leaves and roots also showed specific gene expression responses. While the genes involved in protein metabolism were mainly regulated in leaves, the direct contact of oil with roots induced hypoxia- and ethylene-related responses in L. racemosa roots. Additionally,

A

the assembled L. racemosa transcriptome allowed us to provide suitable reference genes for normalization of qPCR expression data from an oil-exposed mangrove plant species. Finally, improving the understanding of plant molecular responses to oil is an important step to provide tools for biomonitoring and bioremediation after an oil spill.

Author’s contribution:

17

FAFG: data analysis, manuscript preparation PBR: experimental design and execution FG: qPCR experiment and analysis JESP: MF380 oil supply MWW: experimental design and execution MFN: data analysis

IP T

FR: physiological characterization of plant material RSP: physiological characterization of plant material

SC R

MAF: experimental design, manuscript preparation

U

Declarations of interest:

A

N

There are no conflicts of interest to declare.

M

Acknowledgments

ED

The authors would like to thank Petrobras for carefully providing the oil MF380 used in the experiment and to Prof. Wanke Zhang for kindly providing the M. pinnata sequences. We

PT

also acknowledge Sarah Muniz Nardeli and Carolina Farias Saad for their contributions to the experiments. The author's research was supported by Conselho Nacional de

CC E

Desenvolvimento Científico e Tecnológico (CNPq) (Grant 308832 / 2017-5 to MA-F); Fundação de Amparo à Pesquisa do Rio de Janeiro (FAPERJ) (E-26 / 203.039 / 2016 and E-26/010.001942/2014) and Coordenação de Aperfeiçoamento de Pessoal de Nível

A

Superior (CAPES). F.A.F. Guedes is grateful to CAPES for partial support.

References

18

Ahammed, G.J., Wang, M.M., Zhou, Y.H., Xia, X.J., Mao, W.H., Shi, K., Yu, J.Q., 2012a. The growth, photosynthesis and antioxidant defense responses of five vegetable crops to phenanthrene stress. Ecotoxicol. Environ. Saf. 80, 132–139. https://doi.org/10.1016/j.ecoenv.2012.02.015 Ahammed, G.J., Yuan, H.L., Ogweno, J.O., Zhou, Y.H., Xia, X.J., Mao, W.H., Shi, K., Yu, J.Q., 2012b. Brassinosteroid alleviates phenanthrene and pyrene phytotoxicity by increasing detoxification activity and photosynthesis in tomato. Chemosphere 86, 546– 555.https://doi.org/10.1016/j.chemosphere.2011.10.038 Ali, B.A.A., Ali, H.. H., Shaker, G.. A., 2009. The Impact of Ascending Levels of Crude Oil Pollution on Growth of Olive (Olea europaea Linn) Seedlings. Ibn Al-Haitham J. Pure Appl. Sci. 22.

IP T

Alon, U., 2007. Network motifs: Theory and experimental approaches. Nat. Rev. Genet. 8, 450–461. https://doi.org/10.1038/nrg2102

SC R

Andersen, C.L., Jensen, J.L., Ørntoft, T.F., 2004. Normalization of real-time quantitative reverse transcription-PCR data: A model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 64, 5245–5250. https://doi.org/10.1158/0008-5472.CAN-04-0496

N

U

Artico, S., Nardeli, S.M., Brilhante, O., Grossi-de-Sa, M.F., Alves-Ferreira, M., 2010. Identification and evaluation of new reference genes in Gossypium hirsutum for accurate normalization of real-time quantitative RT-PCR data. BMC Plant Biol. 10, 49. https://doi.org/10.1186/1471-222910-49

A

Aschan, G., Pfanz, H., 2003. Non-foliar photosynthesis - A strategy of additional carbon acquisition. Flora 198, 81–97. https://doi.org/10.1078/0367-2530-00080

ED

M

Babu, M.M., Luscombe, N.M., Aravind, L., Gerstein, M., Teichmann, S.A., 2004. Structure and evolution of transcriptional regulatory networks. Curr. Opin. Struct. Biol. 14, 283–291. https://doi.org/10.1016/j.sbi.2004.05.004

PT

Barbazuk, W.B., Fu, Y., McGinnis, K.M., 2008. Genome-wide analyses of alternative splicing in plants: Opportunities and challenges. Genome Res. 18, 1381–1392. https://doi.org/10.1101/gr.053678.106

CC E

Bejarano, A.C., Michel, J., 2016. Oil spills and their impacts on sand beach invertebrate communities: A literature review. Environ. Pollut. 218, 709–722. https://doi.org/10.1016/j.envpol.2016.07.065 Benfield, S.L., Guzman, H.M., Mair, J.M., 2005. Temporal mangrove dynamics in relation to coastal development in Pacific Panama. J. Environ. Manage. 76, 263–276. https://doi.org/10.1016/j.jenvman.2005.02.004

A

Beyer, J., Trannum, H.C., Bakke, T., Hodson, P. V., Collier, T.K., 2016. Environmental effects of the Deepwater Horizon oil spill: A review. Mar. Pollut. Bull. 110, 28–51. https://doi.org/10.1016/j.marpolbul.2016.06.027 Burda, Z., Krzywicki, A., Martin, O.C., Zagorski, M., 2011. Motifs emerge from function in model gene regulatory networks. Proc. Natl. Acad. Sci. 108, 17263–17268. https://doi.org/10.1073/pnas.1109435108 Burns, K.A., Garrity, S.D., Levings, S.C., 1993. How many years until mangrove ecosystems recover from catastrophic oil spills? Mar. Pollut. Bull. 26, 239–248. https://doi.org/10.1016/0025326X(93)90062-O

19

Carvalho, C.M., Santos, A.A., Pires, S.R., Rocha, C.S., Saraiva, D.I., Machado, J.P.B., Mattos, E.C., Fietto, L.G., Fontes, E.P.B., 2008. Regulated nuclear trafficking of rpL10A mediated by NIK1 represents a defense strategy of plant cells against virus. PLoS Pathog. 4. https://doi.org/10.1371/journal.ppat.1000247 Chen, S., Zhou, R., Huang, Y., Zhang, M., Yang, G., Zhong, C., Shi, S., 2011. Transcriptome sequencing of a highly salt tolerant mangrove species Sonneratia alba using Illumina platform. Mar. Genomics 4, 129–136. https://doi.org/10.1016/j.margen.2011.03.005

IP T

Chen, S.F., Li, M.W., Jing, H.J., Zhou, R.C., Yang, G.L., Wu, W., Fan, Q., Liao, W.B., 2015. De novo transcriptome assembly in Firmiana danxiaensis, a tree species endemic to the danxia landform. PLoS One 10, 1–13. https://doi.org/10.1371/journal.pone.0139373

SC R

Chiang, Y.-H., Zubo, Y.O., Tapken, W., Kim, H.J., Lavanway, A.M., Howard, L., Pilon, M., Kieber, J.J., Schaller, G.E., 2012. Functional Characterization of the GATA Transcription Factors GNC and CGA1 Reveals Their Key Role in Chloroplast Development, Growth, and Division in Arabidopsis. Plant Physiol. 160, 332–348. https://doi.org/10.1104/pp.112.198705

U

Conesa, A., Götz, S., García-Gómez, J.M., Terol, J., Talón, M., Robles, M., 2005. Blast2GO: A universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. https://doi.org/10.1093/bioinformatics/bti610

N

Cui, B., Zhang, X., Han, G., Li, K., 2016. Antioxidant Defense Response and Growth Reaction of Amorpha fruticosa Seedlings in Petroleum-Contaminated Soil. Water. Air. Soil Pollut. 227. https://doi.org/10.1007/s11270-016-2821-3

M

A

Czechowski, T., Stitt, M., Altmann, T., Udvardi, M.K., 2005. Genome-Wide Identification and Testing of Superior Reference Genes for Transcript Normalization 139, 5–17. https://doi.org/10.1104/pp.105.063743.1

ED

Dabestani, R., Ivanov, I.N., 1999. A compilation of physical, spectroscopic and photophysical properties of polycyclic aromatic hydrocarbons. Photochem. Photobiol. 70, 10–34. https://doi.org/10.1111/j.1751-1097.1999.tb01945.x

PT

Dang, Z., Zheng, L., Wang, J., Gao, Z., Wu, S., Qi, Z., Wang, Y., 2013. Transcriptomic profiling of the salt-stress response in the wild recretohalophyte Reaumuria trigyna. BMC Genomics 14, 29. https://doi.org/10.1186/1471-2164-14-29

CC E

Dassanayake, M., Haas, J.S., Bohnert, H.J., Cheeseman, J.M., 2010. Comparative transcriptomics for mangrove species: An expanding resource. Funct. Integr. Genomics 10, 523–532. https://doi.org/10.1007/s10142-009-0156-5

A

Dassanayake, M., Haas, J.S., Bohnert, H.J., Cheeseman, J.M., 2009. Shedding light on an extremophile lifestyle through transcriptomics. New Phytol. 183, 764–775. https://doi.org/10.1111/j.1469-8137.2009.02913.x Dromgoole, F.I., 1988. Carbon dioxide fixation in aerial roots of the new zealand mangrove Avicennia marina var. Resinifera (note). New Zeal. J. Mar. Freshw. Res. 22, 617–619. https://doi.org/10.1080/00288330.1988.9516331 Duke, N.C., 2016. Oil spill impacts on mangroves: Recommendations for operational planning and action based on a global review. Mar. Pollut. Bull. 109, 700–715. https://doi.org/10.1016/j.marpolbul.2016.06.082

20

Ellison, A., Farnsworth, E., Moore, G., 2010. Laguncularia racemosa., The IUCN Red List of Threatened Species 2010: e.T178798A7609219. http://dx.doi.org/10.2305/IUCN.UK.20102.RLTS.T178798A7609219.en Ferreyra, M.L.F., Casadevall, R., Luciani, M.D., Pezza, A., Casati, P., 2013. New evidence for differential roles of l10 ribosomal proteins from Arabidopsis. Plant Physiol. 163, 378–91. https://doi.org/10.1104/pp.113.223222 Ferreyra, M.L.F., Pezza, A., Biarc, J., Burlingame, A.L., Casati, P., 2010. Plant L10 Ribosomal Proteins Have Different Roles during Development and Translation under Ultraviolet-B Stress. Plant Physiol. 153, 1878–1894. https://doi.org/10.1104/pp.110.157057

IP T

Field, C.B., Osborn, J.G., Hoffman, L.L., Polsenberg, J.F., Ackerly, D.D., Berry, J.A., Björkman, O., Held, A., Matson, P.A., Mooney, H.A., 1998. Mangrove biodiversity and ecosystem function. Glob. Ecol. Biogeogr. 7, 3–14. https://doi.org/10.2307/2997693

SC R

Flores, H.E., Dai, Y., Cuello, J.L., Maldonado-Mendoza, I.E., Loyola-Vargas, V.M., 1993. Green Roots: Photosynthesis and Photoautotrophy in an Underground Plant Organ. Plant Physiol. 101, 363–371. https://doi.org/101/2/363 [pii]

U

Franke, R., Schreiber, L., 2007. Suberin - a biopolyester forming apoplastic plant interfaces. Curr. Opin. Plant Biol. 10, 252–259. https://doi.org/10.1016/j.pbi.2007.04.004

A

N

Gabardo, I.T., Meniconi, M.F.G., Carneiro, M.E.R., Barbanti, S.M., Falcão, L. V, Platte, E.B., 2003. Oil Spills in a Tropical Country – Brazilian Case Studies. Int. Oil Spill Conf. 3, 1–12. https://doi.org/10.7901/2169-3358-2003-1-1039

M

Gachon, C., Mingam, A., Charrier, B., 2004. Real-time PCR: What relevance to plant studies? J. Exp. Bot. 55, 1445–1454. https://doi.org/10.1093/jxb/erh181

ED

Giuntoli, B., Licausi, F., van Veen, H., Perata, P., 2017. Functional Balancing of the Hypoxia Regulators RAP2.12 and HRA1 Takes Place in vivo in Arabidopsis thaliana Plants. Front. Plant Sci. 8, 1–10. https://doi.org/10.3389/fpls.2017.00591

PT

Hellemans, J., Mortier, G., De Paepe, A., Speleman, F., Vandesompele, J., 2007. qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 8, R19. https://doi.org/10.1186/gb-2007-8-2-r19

CC E

Hensel, P., Proffitt, E.C., Delgado, P., Shigenaka, G., Yender, R., Hoff, R., Mearns, A.J., 2014. Oil spills in mangroves: Planning & response considerations. Hester, M.W., Willis, J.M., Rouhani, S., Steinhoff, M.A., Baker, M.C., 2016. Impacts of the Deepwater Horizon oil spill on the salt marsh vegetation of Louisiana. Environ. Pollut. 216, 361–370. https://doi.org/10.1016/j.envpol.2016.05.065

A

Hruz, T., Laule, O., Szabo, G., Wessendorp, F., Bleuler, S., Oertle, L., Widmayer, P., Gruissem, W., Zimmermann, P., 2008. Genevestigator v3: a reference expression database for the metaanalysis of transcriptomes. Adv. Bioinformatics 2008, 420747. https://doi.org/10.1155/2008/420747 Huang, J., Lu, X., Yan, H., Chen, S., Zhang, W., Huang, R., Zheng, Y., 2012. Transcriptome characterization and sequencing-based identification of salt-responsive genes in Millettia pinnata, a semi-mangrove plant. DNA Res. 19, 195–207. https://doi.org/10.1093/dnares/dss004

21

Huang, J., Lu, X., Zhang, W., Huang, R., Chen, S., Zheng, Y., 2014. Transcriptome sequencing and analysis of leaf tissue of Avicennia marina using the Illumina platform. PLoS One 9, e108785. https://doi.org/10.1371/journal.pone.0108785 Huang, X.D., Mcconkey, B.J., Babu, T.S., Greenberg, B.M., 1997. Mechanisms of photoinduced toxicity of photomodified anthracene to plants: Inhibition of photosynthesis in the aquatic higher plant Lemna gibba (duckweed). Environ. Toxicol. Chem. 16, 1707–1715. https://doi.org/10.1002/etc.5620160819

IP T

Huang, X.D., Zeiler, L.F., Dixon, D.G., Greenberg, B.M., 1996. Photoinduced toxicity of PAHs to the foliar regions of Brassica napus (canola) and Cucumbis sativus (cucumber) in simulated solar radiation. Ecotoxicol. Environ. Saf. 35, 190–197. https://doi.org/10.1006/eesa.1996.0099

SC R

Iseli, C., Jongeneel, C. V, Bucher, P., 1999. ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences. Proc. Int. Conf. Intell. Syst. Mol. Biol. 138–148. Jacob, P., Hirt, H., Bendahmane, A., 2017. The heat-shock protein/chaperone network and multiple stress resistance. Plant Biotechnol. J. 15, 405–414. https://doi.org/10.1111/pbi.12659

U

Jansson, S., 1999. A guide to the Lhc genes and their relatives in Arabidopsis. Trends Plant Sci. 4, 236–240. https://doi.org/10.1016/S1360-1385(99)01419-3

A

N

Jin, L., Che, X., Zhang, Z., Li, Y., Gao, H., Zhao, S., 2017. The mechanisms by which phenanthrene affects the photosynthetic apparatus of cucumber leaves. Chemosphere 168, 1498–1505. https://doi.org/10.1016/j.chemosphere.2016.12.002

M

Kang, F., Chen, D., Gao, Y., Zhang, Y., 2010. Distribution of polycyclic aromatic hydrocarbons in subcellular root tissues of ryegrass (Lolium multiflorum Lam.). BMC Plant Biol. 10, 1–7. https://doi.org/10.1186/1471-2229-10-210

ED

Ke, L., Zhang, C., Wong, Y.S., Tam, N.F.Y., 2011. Dose and accumulative effects of spent lubricating oil on four common mangrove plants in South China. Ecotoxicol. Environ. Saf. 74, 55–66. https://doi.org/10.1016/j.ecoenv.2010.09.011

PT

Krishnamurthy, P., Ranathunge, K., Franke, R., Prakash, H.S., Schreiber, L., Mathew, M.K., 2009. The role of root apoplastic transport barriers in salt tolerance of rice (Oryza sativa L.). Planta 230, 119–134. https://doi.org/10.1007/s00425-009-0930-6

CC E

Krishnamurthy, P., Tan, X.F., Lim, T.K., Lim, T.M., Kumar, P.P., Loh, C.S., Lin, Q., 2014. Proteomic analysis of plasma membrane and tonoplast from the leaves of mangrove plant Avicennia officinalis. Proteomics 14, 2545–2557. https://doi.org/10.1016/j.dib.2015.10.016

A

Kummerová, M., Kmentová, E., 2004. Photoinduced toxicity of fluoranthene on germination and early development of plant seedling. Chemosphere 56, 387–393. https://doi.org/10.1016/j.chemosphere.2004.01.007 Kummerová, M., Krulová, J., Zezulka, Š., Tříska, J., 2006. Evaluation of fluoranthene phytotoxicity in pea plants by Hill reaction and chlorophyll fluorescence. Chemosphere 65, 489–496. https://doi.org/10.1016/j.chemosphere.2006.01.052 Kummerová, M., Váňová, L., 2007. Chlorophyll fluorescence as an indicator of fluoranthene phototoxicity. Plant, Soil Environ. 53, 430–436.

22

Langangen, Olsen, E., Stige, L.C., Ohlberger, J., Yaragina, N.A., Vikebø, F.B., Bogstad, B., Stenseth, N.C., Hjermann, D., 2017. The effects of oil spills on marine fish: Implications of spatial variation in natural mortality. Mar. Pollut. Bull. 119, 102–109. https://doi.org/10.1016/j.marpolbul.2017.03.037 Lewis, M., Pryor, R., 2013. Toxicities of oils, dispersants and dispersed oils to algae and aquatic plants: Review and database value to resource sustainability. Environ. Pollut. 180, 345–367. https://doi.org/10.1016/j.envpol.2013.05.001

IP T

Lewis, M., Pryor, R., Wilking, L., 2011. Fate and effects of anthropogenic chemicals in mangrove ecosystems: A review. Environ. Pollut. 159, 2328–2346. https://doi.org/10.1016/j.envpol.2011.04.027 Li, H., Durbin, R., 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760. https://doi.org/10.1093/bioinformatics/btp324

SC R

Li, J.H., Gao, Y., Wu, S.C., Cheung, K.C., Wang, X.R., Wong, M.H., 2008. Physiological and Biochemical Responses of Rice (Oryza Sativa L.) to Phenanthrene and Pyrene. Int. J. Phytoremediation 10, 106–118. https://doi.org/10.1080/16226510490454803

N

U

Licausi, F., Weits, D.A., Pant, B.D., Scheible, W.R., Geigenberger, P., van Dongen, J.T., 2011. Hypoxia responsive gene expression is mediated by various subsets of transcription factors and miRNAs that are determined by the actual oxygen availability. New Phytol. 190, 442–456. https://doi.org/10.1111/j.1469-8137.2010.03451.x

M

A

Mastrangelo, A.M., Marone, D., Laidò, G., De Leonardis, A.M., De Vita, P., 2012. Alternative splicing: Enhancing ability to cope with stress via transcriptome plasticity. Plant Sci. 185–186, 40–49. https://doi.org/10.1016/j.plantsci.2011.09.006

ED

McKee, K.L., 1996. Growth and physiological responses of neotropical mangrove seedlings to root zone hypoxia. Tree Physiol. 16, 883–889. https://doi.org/10.1093/treephys/16.11-12.883

PT

Meniconi, M. de F.G., Gabardo, I.T., Carneiro, M.E.R., Barbanti, S.M., Silva, G.C. da, Massone, C.G., 2002. Brazilian oil spills chemical characterization - Case studies. Environ. Forensics 3, 303–321. https://doi.org/10.1006/enfo.2002.0101

CC E

Meudec, A., Dussauze, J., Deslandes, E., Poupart, N., 2006. Evidence for bioaccumulation of PAHs within internal shoot tissues by a halophytic plant artificially exposed to petroleum-polluted sediments. Chemosphere 65, 474–481. https://doi.org/10.1016/j.chemosphere.2006.01.058 Milo, R., Shen-Orr, S., Itzkovitz, S., Kashtan, N., Chklovskii, D., Alon, U., 2002. Network motifs: Simple building blocks of complex networks. Science (80-. ). 298, 824–827. https://doi.org/10.1126/science.298.5594.824

A

Moin, M., Bakshi, A., Saha, A., Dutta, M., Madhav, S.M., Kirti, P.B., 2016. Rice Ribosomal Protein Large Subunit Genes and Their Spatio-temporal and Stress Regulation. Front. Plant Sci. 7, 1– 20. https://doi.org/10.3389/fpls.2016.01284 Nagaraj, S., Senthil-Kumar, M., Ramu, V.S., Wang, K., Mysore, K.S., 2016. Plant Ribosomal Proteins, RPL12 and RPL19, Play a Role in Nonhost Disease Resistance against Bacterial Pathogens. Front. Plant Sci. 6, 1–10. https://doi.org/10.3389/fpls.2015.01192 Nagelkerken, I., Blaber, S.J.M., Bouillon, S., Green, P., Haywood, M., Kirton, L.G., Meynecke, J.O., Pawlik, J., Penrose, H.M., Sasekumar, A., Somerfield, P.J., 2008. The habitat function of

23

mangroves for terrestrial and marine fauna: A review. Aquat. Bot. 89, 155–185. https://doi.org/10.1016/j.aquabot.2007.12.007 Naidoo, G., 2016. Mangrove propagule size and oil contamination effects: Does size matter? Mar. Pollut. Bull. 110, 362–370. https://doi.org/10.1016/j.marpolbul.2016.06.040 Naidoo, G., Naidoo, K., 2017. Are pioneer mangroves more vulnerable to oil pollution than later successional species? Mar. Pollut. Bull. 121, 135–142. https://doi.org/10.1016/j.marpolbul.2017.05.067

IP T

Naidoo, G., Naidoo, K., 2016. Uptake of polycyclic aromatic hydrocarbons and their cellular effects in the mangrove Bruguiera gymnorrhiza. Mar. Pollut. Bull. 113, 193–199. https://doi.org/10.1016/j.marpolbul.2016.09.012

SC R

Nardeli, S.M., Saad, C.F., Rossetto, P. de B., Caetano, V.S., Ribeiro-Alves, M., Paes, J.E.S., Danielowski, R., Maia, L.C. da, Oliveira, A.C. de, Peixoto, R.S., Reinert, F., Alves-Ferreira, M., 2016. Transcriptional responses of Arabidopsis thaliana to oil contamination. Environ. Exp. Bot. 127, 63–72. https://doi.org/10.1016/j.envexpbot.2016.03.007

U

Petri, D.J.C., Bernini, E., Souza, L.M. de, Rezende, C.E., 2011. Species distribution and structure of mangrove of the Benevente River, Anchieta, ES. Biota Neotrop. 11, 107–116. https://doi.org/10.1590/S1676-06032011000300009

A

N

Pezeshki, S.R., Hester, M.W., Lin, Q., Nyman, J.A., 2000. The effects of oil spill and clean-up on dominant US Gulf coast marsh macrophytes: A review. Environ. Pollut. 108, 129–139. https://doi.org/10.1016/S0269-7491(99)00244-4

M

Pfaffl, M.W., Horgan, G.W., Dempfle, L., 2002. Relative expression software tool (REST©) for groupwise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 30.

ED

Proffitt, C.E., Devlin, D.J., Lindseyt, M., 1995. Effects of Oil on Mangrove Seedlings Grown Under Different Environmental Conditions. Mar. Pollut. Bull. 30, 788–793. https://doi.org/10.1016/0025-326X(95)00070-4

PT

Record, S., Charney, N.D., Zakaria, R.M., Ellison, A.M., 2013. Projecting global mangrove species and community distributions under climate change. Ecosphere 4, 1–23. https://doi.org/10.1890/ES12-00296.1

CC E

Reef, R., Lovelock, C.E., 2015. Regulation of water balance in Mangroves. Ann. Bot. 115, 385–395. https://doi.org/10.1093/aob/mcu174

A

Reinert, F., Pinho, C.F. De, Ferreira, M.A., 2016. Diagnosing the level of stress on a mangrove species (Laguncularia racemosa) contaminated with oil - a necessary step for monitoring mangrove ecosystems. Mar. Pollut. Bull. 6. https://doi.org/10.1016/j.marpolbul.2016.08.070 Santos, L.C.M., Cunha-Lignon, M., Schaeffer-Novelli, Y., Cintrón-Molero, G., 2012. Long-term effects of oil pollution in mangrove forests (Baixada Santista, Southeast Brazil) detected using a GIS-based multitemporal analysis of aerial photographs. Brazilian J. Oceanogr. 60, 159–170. https://doi.org/10.1590/S1679-87592012000200006 Schulz, M.H., Zerbino, D.R., Vingron, M., Birney, E., 2012. Oases: Robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics 28, 1086–1092. https://doi.org/10.1093/bioinformatics/bts094

24

Silveira, E.D., Alves-Ferreira, M., Guimarães, L.A., da Silva, F.R., Carneiro, V.T.D.C., 2009. Selection of reference genes for quantitative real-time PCR expression studies in the apomictic and sexual grass Brachiaria brizantha. BMC Plant Biol. 9, 84. https://doi.org/10.1186/14712229-9-84 Sormani, R., Masclaux-Daubresse, C., Daniele-Vedele, F., Chardon, F., 2011. Transcriptional Regulation of Ribosome Components Are Determined by Stress According to Cellular Compartments in Arabidopsis thaliana. PLoS One 6. https://doi.org/10.1371/journal.pone.0028070

IP T

Srikanth, S., Lum, S.K.Y., Chen, Z., 2015. Mangrove root: adaptations and ecological importance. Trees - Struct. Funct. 30, 451–465. https://doi.org/10.1007/s00468-015-1233-0

SC R

Staiger, D., Brown, J.W.S., 2013. Alternative Splicing at the Intersection of Biological Timing, Development, and Stress Responses. Plant Cell 25, 3640–3656. https://doi.org/10.1105/tpc.113.113803 Surget-Groba, Y., Montoya-Burgos, J., 2010. Optimization of de novo transcriptome assembly from next-generation sequencing data. Genome Res. 20, 1432–1440. https://doi.org/10.1101/gr.103846.109.1432

N

U

Swindell, W.R., Huebner, M., Weber, A.P., 2007. Transcriptional profiling of Arabidopsis heat shock proteins and transcription factors reveals extensive overlap between heat and non-heat stress response pathways. BMC Genomics 8, 1–15. https://doi.org/10.1186/1471-2164-8-125

M

A

Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., Simonovic, M., Roth, A., Santos, A., Tsafou, K.P., Kuhn, M., Bork, P., Jensen, L.J., Von Mering, C., 2015. STRING v10: Protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–D452. https://doi.org/10.1093/nar/gku1003

ED

Taiz, L., Zeiger, E., 2002. Plant physiology, third. ed. Sinauer asssociates, Sunderland.

PT

Tao, S., Jiao, X.C., Chen, S.H., Liu, W.X., Coveney, R.M., Zhu, L.Z., Luo, Y.M., 2006. Accumulation and distribution of polycyclic aromatic hydrocarbons in rice (Oryza sativa). Environ. Pollut. 140, 406–415. https://doi.org/10.1016/j.envpol.2005.08.004

CC E

Tomar, R.S., Jajoo, A., 2015. Photomodified fluoranthene exerts more harmful effects as compared to intact fluoranthene by inhibiting growth and photosynthetic processes in wheat. Ecotoxicol. Environ. Saf. 122, 31–36. https://doi.org/10.1016/j.ecoenv.2015.07.002 Tomar, R.S., Jajoo, A., 2014. Fluoranthene, a polycyclic aromatic hydrocarbon, inhibits light as well as dark reactions of photosynthesis in wheat (Triticum aestivum). Ecotoxicol. Environ. Saf. 109, 110–115. https://doi.org/10.1016/j.ecoenv.2014.08.009

A

Tomar, R.S., Jajoo, A., 2013. A quick investigation of the detrimental effects of environmental pollutant polycyclic aromatic hydrocarbon fluoranthene on the photosynthetic efficiency of wheat (Triticum aestivum). Ecotoxicology 22, 1313–1318. https://doi.org/10.1007/s10646-0131118-1 Tomlinson, P.B., 1986. The botany of mangroves, second. ed. Cambridge University Press, New York. Trapnell, C., Williams, B., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M.J., Salzberg, S.L., Wold, B.J., Pachter, L., 2010. Transcript assembly and quantification by RNA-Seq reveals

25

unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. https://doi.org/10.1038/nbt.1621 Untergasser, A., Cutcutache, I., Koressaar, T., Ye, J., Faircloth, B.C., Remm, M., Rozen, S.G., 2012. Primer3-new capabilities and interfaces. Nucleic Acids Res. 40, 1–12. https://doi.org/10.1093/nar/gks596 Valiela, I., Bowen, J.L., York, J.K., 2001. Mangrove Forests: One of the World’s Threatened Major Tropical Environments. Bioscience 51, 807. https://doi.org/10.1641/00063568(2001)051[0807:MFOOTW]2.0.CO;2

IP T

Vandesompele, J., De Preter, K., Pattyn, F., Poppe, B., Van Roy, N., De Paepe, A., Speleman, F., 2002. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 3, RESEARCH0034. https://doi.org/10.1186/gb2002-3-7-research0034

SC R

Wang, E., Purisima, E.O., 2005. Network motifs are enriched with transcription factors whose transcripts have short half-lives. Trends Genet. 21, 492–495. https://doi.org/10.1016/j.tig.2005.07.004

U

Wang, W., Vinocur, B., Shoseyov, O., Altman, A., 2004. Role of plant heat-shock proteins and molecular chaperones in the abiotic stress response. Trends Plant Sci. 9, 244–252. https://doi.org/10.1016/j.tplants.2004.03.006

M

A

N

Wang, Y., Wang, S., Luo, C., Xu, Y., Pan, S., Li, J., Ming, L., Zhang, G., Li, X., 2015. Influence of rice growth on the fate of polycyclic aromatic hydrocarbons in a subtropical paddy field: A life cycle study. Chemosphere 119, 1233–1239. https://doi.org/10.1016/j.chemosphere.2014.09.104

ED

Warner, J.R., McIntosh, K.B., 2009. How Common are Extra-ribosomal Functions of Ribosomal Proteins? Mollecular Cell 34, 3–11. https://doi.org/10.1016/j.molcel.2009.03.006

PT

Weisman, D., Alkio, M., Colón-Carmona, A., 2010. Transcriptional responses to polycyclic aromatic hydrocarbon-induced stress in Arabidopsis thaliana reveal the involvement of hormone and defense signaling pathways. BMC Plant Biol. 10. https://doi.org/https://doi.org/10.1186/14712229-10-59

CC E

Wilwerth, M.W., Rossetto, P. de B., Reinert, F., Peixoto, R.S., Alves-Ferreira, M., 2016. Efficient rna extraction protocol for the wood mangrove species Laguncularia racemosa suited for nextgeneration RNA sequencing. Pakistan J. Bot. 48, 661–672.

A

Wu, S., Zhu, Z., Fu, L., Niu, B., Li, W., 2011. WebMGA: a customizable web server for fast metagenomic sequence analysis. BMC Genomics 12, 444. https://doi.org/10.1186/1471-216412-444 Xie, C., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., Kong, L., Gao, G., Li, C.Y., Wei, L., 2011. KOBAS 2.0: A web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, 316–322. https://doi.org/10.1093/nar/gkr483 Yang, Y., Yang, S., Li, J., Deng, Y., Zhang, Z., Xu, S., Guo, W., Zhong, C., Zhou, R., Shi, S., 2015. Transcriptome analysis of the Holly mangrove Acanthus ilicifolius and its terrestrial relative, Acanthus leucostachyus, provides insights into adaptation to intertidal zones. BMC Genomics 16, 1–12. https://doi.org/10.1186/s12864-015-1813-9

26

Ye, J., Fang, L., Zheng, H., Zhang, Y., Chen, J., Zhang, Z., Wang, J., Li, S., Li, R., Bolund, L., Wang, J., 2006. WEGO: A web tool for plotting GO annotations. Nucleic Acids Res. 34, 293–297. https://doi.org/10.1093/nar/gkl031 Zerbino, D.R., Birney, E., 2008. Velvet: Algorithms for de novo short read assembly using de Bruijn graphs. Genome Res. 18, 821–829. https://doi.org/10.1101/gr.074492.107 Zhang, C.G., Leung, K.K., Wong, Y.S., Tam, N.F.Y., 2007. Germination, growth and physiological responses of mangrove plant (Bruguiera gymnorrhiza) to lubricating oil pollution. Environ. Exp. Bot. 60, 127–136. https://doi.org/10.1016/j.envexpbot.2006.09.002

A

CC E

PT

ED

M

A

N

U

SC R

IP T

Zhao, S., Fernald, R.D., 2009. Comprehensive Algorithm for Quantitative Real-Time Polymerase Chain Reaction. J. Comput. Biol. 12, 1047–1064. https://doi.org/10.1089/cmb.2005.12.1047

27

Tables Table 1. General statistics for the L. racemosa reference transcriptome. This transcriptome assembly was produced by an approach combining Velvet and Oases software with kmer=87. Statistics

46,944

IP T

Contigs number

Total length (bp)

55,349,684

1,179

SC R

Mean length (bp)

N50 (bp)

Counts until N50

1,789

10,433

N

U

Minimum length (bp)

A

Maximum length (bp)

Number of contigs with nr annotation

M

Number of contigs ≥1kb

163

11,799

21,556

34,683

23,021

Number of contigs with KOG annotation

19,628

PT

ED

Number of contigs with GO annotation

6,021

A

CC E

Number of contigs with KEGG Pathway annotation

28

Table 2. Comparison of mangrove and semi-mangrove assembled transcriptomes.

Laguncularia racemosa

Illumina HiSeq 2000 Platform

Velvet and Oases

454/Roche GSFLX

Newbler and Phrap

Rhizophora mangle

Reference

169.84

46,944

1,179

This study

0.23

25,386

433

0.31

31,548

360

Dassanayake et al., 2009

15.51

30,628

661

Chen et al., 2011

Millettia pinnata

Illumina Genome SOAPdenovo Analyzer IIx

192

108,598

606

Huang et al., 2012

Avicennia marina

Illumina Genome SOAPdenovo Analyzer IIx

91,125

N

463

Huang et al., 2014

34.13

73,039

580

38.53

69,580

913

Acanthus ilicifolius

40

Trinity

Yang et al., 2015

A

CC E

PT

ED

Acanthus leucostachyus

Illumina HiSeq 2000 Platform

U

Illumina Genome Abyss, Velvet Analyzer and Edena

M

Sonneratia alba

A

Heritiera littoralis

#Contigs

IP T

Assemblers

Species

Mean length (bp)

#Read (million)

SC R

Sequencing technology

29

Table 3. KEGG pathways annotated for L. racemosa transcriptome. Asterisk indicates statistically significant pathways.

Rank

KEGG Pathway

#L. racemosa transcripts

Corrected p-

(%)

value 0.21296071

Metabolic pathways (ath01100)

2701 (43.56)

2

Biosynthesis of secondary metabolites

1335 (21.53)

3

Ribosome (ath03010)

315 (5.08)

IP T

1

4

Plant hormone signal transduction (ath04075)

267 (4.31)

0.99999998

5

Glycolysis / Gluconeogenesis (ath00010)

262 (4.23)

0.42188958

6

Plant-pathogen interaction (ath04626)

249 (4.02)

0.99999998

7

Purine metabolism (ath00230)

8

Spliceosome* (ath03040)

9

0.99999998

U

SC R

(ath01110)

0.99999998

N

234 (3.77)

0.09360004 0.00000068

Starch and sucrose metabolism (ath00500)

207 (3.34)

0.85951972

10

RNA transport (ath03013)

206 (3.32)

0.29212087

12

Pyrimidine metabolism* (ath00240)

185 (2.98)

0.00772265

44

Propanoate metabolism* (ath00640)

87 (1.40)

0.00825595

A

CC E

PT

ED

M

A

215 (3.47)

30

Table 4. L. racemosa oil-responsive root genes functionally related to the ethylene pathway, with their putative Arabidopsis homologs, log2FC values, and FPKM values estimated from RNASeq. Asterisk indicates that values shown for the gene represent the average of the transcripts from the same locus.

FPKM

L. racemosa

Ath gene name/description

locus ID

Log2FC

Locus_5692

ERF1B (AT3G23240); ethylene responsive TF

Locus_11039

ERF106 (AT5G07580); ethylene responsive TF

Locus_22841

EIN4 (AT3G04580); ethylene receptor

Locus_33481

3.429

2.505

SC R

ERF17 (AT1G19210); ethylene responsive TF

26.975

6.488

27.911

1.982

2.316

9.146

1.871

5.456

19.963

RAP2.1 (AT1G46768); ethylene responsive TF

1.804

8.111

28.317

Locus_3084

ERF105 (AT5G51190); ethylene responsive TF

1.752

6.691

22.540

Locus_2424

ERF1A (AT4G17500); ethylene responsive TF

1.742

14.823

49.588

Locus_35375

ACO1 (AT2G19590); ethylene biosynthesis

1.613

107.547

329.042

Locus_32033*

ERF71 (AT2G47520); ethylene responsive TF

1.526

39.859

114.890

Locus_13549

ERF114 (AT5G61890); ethylene responsive TF

1.462

6.867

18.917

Locus_2661

ETR2 (AT3G23150); ethylene signaling

1.438

13.256

35.915

ERF72 (AT3G16770); ethylene responsive TF

1.390

205.843

539.299

N

A

M

ED

CC E

A

Locus_2001

U

2.105

PT

Locus_30244

IP T

Control Oil Treated

31

Figures

Figure 1. KOG functional classification of L. racemosa transcriptome. The x-axis represents KOG classes (25), and the y-axis represents the number of sequences annotated in each class. (2-

A

CC E

PT

ED

M

A

N

U

SC R

IP T

column fitting; color online only)

32

Figure 2. Enriched GO categories found for oil-responsive transcripts of L. racemosa leaves and roots. The x-axis shows the number of sequences annotated into MF (molecular function), CC

A

CC E

PT

ED

M

A

N

U

SC R

IP T

(cellular component) and BP (biological process) categories. (2-column fitting; color online only)

33

Figure 3. Comparison of L. racemosa leaf oil-responsive genes with heat stress expression data. The log2FC values calculated for L. racemosa loci were compared to those of Arabidopsis genes in heat stress studies available in the Genevestigator collection through their assigned ATcodes. The log2FC values are shown in a color-coded scale. The description of the heat experiments used to generate this heatmap can be found in Table S8, and the complete heatmap is in Figure S4. (2-

A

CC E

PT

ED

M

A

N

U

SC R

IP T

column fitting; color in print)

34

Figure 4. Comparison of L. racemosa root oil-responsive genes with hypoxia stress expression data. The log2FC values calculated for L. racemosa loci were compared to those of Arabidopsis genes in hypoxia stress studies available in the Genevestigator collection through their assigned ATcodes. The log2FC values are shown in a color-coded scale. The description of the hypoxia experiments used to generate this heatmap can be found in Table S8, and the complete heatmap is in Figure S5.

A

CC E

PT

ED

M

A

N

U

SC R

IP T

(2-column fitting; color in print)

35

Figure 5. Comparative analysis of L. racemosa genes regulated by MF380 oil. Venn diagrams show L. racemosa oil-responsive leaf (A) and root genes (B) whose putative homologs were also regulated in A. thaliana by phenanthrene (PHE) (Weisman et al., 2010) and by the water-soluble fraction of MF380 (WSF-MF380) (Nardeli et al., 2016). The detected L. racemosa genes were compared to A. thaliana findings through their assigned ATcodes. The number of shared ATcodes is shown together with the corresponding number of L. racemosa loci in parenthesis. (single-column fitting; color

A

CC E

PT

ED

M

A

N

U

SC R

IP T

online only)

36

Figure 6. L. racemosa oil-responsive leaf (A) and root co-expression networks (B). These networks were constructed by searching the ATcodes assigned to L. racemosa modulated genes against the STRING database. Only high-confidence co-expression evidence (score ≥ 0.7) is shown in these networks. The highest absolute log2FC value represented in these networks was 6.647. (2-column

A

CC E

PT

ED

M

A

N

U

SC R

IP T

fitting; color in print)

37

Figure 7. Expression stability values M (A) and pairwise variation V of the tested candidate reference genes (B). These values were calculated by geNorm using L. racemosa control and oil-exposed leaf expression data. The lower the M value is, the more stable the candidate gene expression. V denotes the optimal number of reference genes necessary for an accurate normalization. Asterisk indicates the optimal number of genes (2) used for normalization of L. racemosa expression data. (single-

A

CC E

PT

ED

M

A

N

U

SC R

IP T

column fitting)

38

Figure 8. RNASeq validation by qPCR analysis. The Ref1 and Ref6 genes were used as internal controls for REST analysis. The qPCR results are the means of three biological replicates, each one

IP T

with three technical replicates. RNASeq oil-responsive genes were identified by Cuffdiff software. Asterisk indicates that the differential expression predicted by RNASeq was confirmed by qPCR.

A

CC E

PT

ED

M

A

N

U

SC R

(1.5-column fitting)

39