Local adaptation of a dominant coastal tree to freshwater availability and solar radiation suggested by genomic and ecophysiological approaches

Local adaptation is often a product of environmental variations in geographical space and has implications for biodiversity conservation. We investigated the role of latitudinal heterogeneity in climate on the organization of genetic and phenotypic variation in the dominant coastal tree Avicennia schaueriana. In a common garden experiment, samples from an equatorial region, with pronounced seasonality in precipitation, accumulated less biomass, and showed lower stomatal conductance and transpiration, narrower xylem vessels, smaller leaves and higher reflectance of long wavelengths by the stem epidermis than samples from a subtropical region, with seasonality in temperature and no dry season. Transcriptomic differences identified between trees sampled under field conditions at equatorial and subtropical sites, were enriched in functional categories such as responses to temperature, solar radiation, water deficit, photosynthesis and cell wall biosynthesis. Remarkably, the diversity based on genome-wide SNPs revealed a north-south genetic structure and signatures of selection were identified for loci associated with photosynthesis, anthocyanin accumulation and the responses to osmotic and hypoxia stresses. Our results suggest the existence of divergence in key resource-use characteristics, likely driven by seasonality in water deficit and solar radiation. These findings provide a basis for conservation plans and for predicting coastal plants responses to climate change.

dispersion capabilities but also by means of natural selection, caused by varying establishment success over the broad environmental heterogeneity across latitudes 18 . For instance, distinct coastal species in the North Atlantic show north-south organization of diversity with evidence of selection for different thermal regimes [19][20][21] . Similarly, in the Southwest Atlantic, an overlapping north-south structure of genetic diversity, with reduced gene flow between northern and southern populations, has been observed in phylogenetically distant habitat-forming species [22][23][24] . In this context, one could expect northerly and southerly populations of these species to adapt differently to contrasting environments over their latitudinal ranges, especially due to their large population sizes 2,25 . However, the neutrality of the molecular markers used in previous studies has precluded inferences regarding adaptive variation. The existence of local adaptation remains virtually unknown in species that play a central role in sustaining the coastal biodiversity in the South Atlantic. This limited knowledge about the organization of non-neutral variation compromises accurate predictions and suitable conservation efforts for sustainable management.
Here, we tested the hypothesis that latitudinal heterogeneity in climate variables shapes the variation in genotypes and phenotypes involved in the optimization of resource use in widely distributed coastal species. To reduce the potential for incorrect conclusions about selection 26 , we integrated three independent but complementary approaches, making predictions as follows: (1) using a common garden experiment, individuals from contrasting latitudes would show genetically based phenotypic divergence in ecophysiological traits; (2) under contrasting latitudes, transcriptomic changes would be detected in genes involved in responses to environmental variation; and (3) signatures of selection along the species distribution would be detected in genes involved in responses to latitudinal variation in air temperature, solar radiation and freshwater availability determined by air water vapor pressure deficit (VPD), rainfall and tidal regime. We considered the new-world coastal tree species Avicennia schaueriana Stapf & Leechman ex Moldenke, found in the Lower Lesser Antilles and widespread from Venezuela to the southernmost temperate mangroves on the Atlantic coast of South America (~28 °S) 27 , as a model to test this hypothesis. The broad latitudinal range of A. schaueriana, spanning wide environmental variation (Fig. 1), and the previously detected north-south structure of neutral variation 22 , facilitate ecological predictions and the accumulation of divergence, motivating this choice. Our results indicated that fluctuations in VPD, rainfall and solar radiation may be associated with the observed phenotypic and genotypic variation as well as the regulation of gene expression in this dominant coastal tree. We discuss the implications of our results for the persistence of coastal biodiversity in the context of climate change.

Materials and Methods
propagule sampling. Mature propagules were collected from 15 Avicennia schaueriana Stapf & Leechman ex Moldenke mother trees that were at least 100 m apart from one another in the following two genetically contrasting populations 22 : (1) the southernmost range limit of American mangroves, in the subtropical region, and (2) an equatorial site in one of the world's largest macrotidal mangrove forests 28,29 , near the northernmost limit of the species range (Fig. 1). We refer to samples collected in the former site as "subtropical" and those in the latter as "equatorial" throughout this work. Sampling in more localities was impracticable due to the absence of mature propagules during fieldwork. A detailed characterization of each of these contrasting sites can be found in Table 1, in Supplementary Methods and in the Supplementary Information Fig. S1. comparative ecophysiology using a common garden experiment. Propagules were germinated as described for Avicennia germinans 30 . We grew propagules in trays with local mangrove soil for two months. After this period, 44 similar-sized seedlings from 30 distinct mother trees, 15 from equatorial and 15 from subtropical The green-shaded area represents the presence of the species; colored points represent sampling sites of propagules, used in a common garden experiment, and of plant tissues, used for both DNA-and RNA-sequencing; black points represent sampling sites of plant tissues, used for genomic DNA-sequencing only; and arrows represent the flow directions of major oceanic currents. Environmental variation across sampling sites was obtained from public databases of climate ('WorldClim' 110,111 ) and tide variables ('Environment Climate Data Sweden' 112 ).
sites -with an average height of 18 cm, most with three leaf pairs and senesced cotyledons -were transplanted to 6 L pots filled with topsoil and sand (1:1). Seedlings were cultivated for seven months under homogenous conditions in a glasshouse at the University of Campinas, São Paulo, Brazil (22°49′S 47°04′W), where automatic sensors coupled with a data logger (Onset Computer Corp.) measured the atmospheric humidity and temperature every 15 minutes. Seedlings were irrigated daily at 10 a.m. and 5 p.m. with a 3-minute freshwater spray. Twice a week, nutrients were added to the soil using 600 mL of 0.4X Hoagland solution with 15.0 g L −1 NaCl per pot. Pots were rotated weekly to reduce the effects of environmental heterogeneity. The environmental conditions in the glasshouse were different from those at each sampling site, as is shown in the Supplementary Information Fig. S2.
The light reflectance of stems was measured in ten plants from each sampling site using a USB4000 spectrophotometer (OceanOptics, Inc.) coupled to a deuterium-halogen light source (DH-2000; OceanOptics) using a light emission range from 200-900 nm. Photosynthesis, stomatal conductance and transpiration rates were measured every 2.0-2.5 hours in five six-month-old individuals from each sampling site on two different days using a Li-Cor 6400 XT (Li-Cor Corp.).
After harvest, three plants without flowers or flower buds from each sampling site were split into leaves, stems and roots, washed with distilled water, dried for 7 days at 70 °C and weighed. The individual leaf area, total leaf area and leaf lamina angle per plant were measured through photographic analyses using ImageJ 31 . The specific leaf area (SLA, cm 2 leaf area kg −1 leaf dry biomass) was also calculated for these samples. Stems were fixed in FAA (Formaldehyde Alcohol Acetic acid), stored in 70% alcohol for wood anatomy analysis and cut into 30-μm thick transverse sections. Sections were stained with a mixture of 1% Astra Blue and 50% alcohol (1:1) followed by 1% Safranin O. Micrographs were taken using an Olympus BX51 microscope coupled to an Olympus DP71 camera (Olympus Corp.). The following wood traits were quantified using ImageJ and R v.4.0.0: vessel lumen area (A), vessel density in xylem (number of vessels/xylem area), vessel grouping index (mean number of vessels per vessel grouping) and vessel lumen area in sapwood (vessel lumen area/sapwood area). The vessel arithmetic diameter (D) was estimated according to Scholz et al. 32 .
Shapiro-Wilk's normality tests and Fisher's F-tests of equality of variances were performed in R v.4.0.0 to determine the suitability of group comparison statistical tests. In comparisons between equatorial and subtropical samples, we used unpaired Student's T-tests or, alternatively, Mann-Whitney-Wilcoxon tests. Multiple-group comparisons were conducted for the analysis of environmental conditions in the field vs. in the glasshouse using one-way analysis of variance (ANOVA) with Tukey's post hoc honest significant difference (HSD) tests. A significance level of 0.05 was used as the alpha for all tests. The Bonferroni adjustment of P-values was conducted for multiple comparisons.
Plant material for RNA extraction and RNA-sequencing. Plant material used for RNA-sequencing (RNA-Seq) was collected from the equatorial and subtropical sites described in the "Propagule sampling" section and immediately stored in RNAlater (Ambion, Inc.). Because precise species identification requires the analysis of both vegetative branches and flower morphology, opportunities for sampling were limited by the reproductive phenology at each visited site. Leaves, stems and flowers from three adult trees at least 100 m apart were collected from July-August of 2014, which corresponds to the end of winter at the subtropical site and the beginning of the dry season at the equatorial site. To minimize the detection of differential transcript expression due to circadian changes, sampling was conducted during low tide and from 10:30 AM to 4:00 PM. A detailed description of the environmental conditions at the time of sampling is available in Supplementary Information Table S1. We extracted RNA according to Oliveira et al. 33 and evaluated its integrity and purity using agarose gel electrophoresis and a NanoVue spectrophotometer (GE Healthcare Life Sciences). Illumina TruSeq RNA Sample Preparation kits (Illumina, Inc.) were used to construct cDNA libraries. cDNA quality was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies) and concentrations were quantified by qPCR using the Sequencing Library qPCR Quantification kit (Illumina, Inc.). Sequencing was performed using two 36-cycle TruSeq SBS paired-end kits (Illumina, Inc.) and the Genome Analyzer IIx platform (Illumina, Inc.).
Assembly and characterization of the A. schaueriana transcriptome. Adapter sequences were trimmed, and 72-bp paired-end reads were filtered by quality (Phred score ≥ 20 for at least 70% of the read length) using the NGS QC Toolkit v.2.3 34 . High-quality reads were used for transcriptome assembly in the CLC Genomics Workbench (https://www.qiagenbioinformatics.com/). We used the default settings, except for the distance between pairs (300-500 bp) and k-mer size (45 bp).
Reads were mapped to transcripts using bowtie1 35 in single-read mode using the default parameters. Transcripts without read-mapping support were removed. Functional annotation was performed using blastx v.2.2.31 36 with an e-value < 1 −10 . NCBI RefSeq. 37 , The Arabidopsis Information Resource (TAIR) 38 and NCBI non-redundant (nr) databases were used as reference databases. We excluded transcripts that were exclusively similar to non-plant sequences. Protein family domains were identified using HMMER3 39 , which iteratively searched transcripts against the Pfam database. To assign Gene Ontology (GO) terms to transcripts, we used the Arabidopsis thaliana gene association file from the GO Consortium 40 and retrieved the information for transcripts showing similar coding sequences in the A. thaliana genome. Redundant transcripts were clustered with Cd-hit-est v.4.6.1 41 using local alignment mode with 95% identity and 70% coverage of the shortest sequence thresholds. Open reading frames (ORFs) were identified using Transdecoder (http://transdecoder.sf.net). We reduced the redundancy of transcripts in the final assembly by retaining for each CD-HIT-EST cluster either the sequence with the longest ORF or, in the absence of sequences containing ORF, the longest sequence.
The completeness of the final transcriptome was assessed using BUSCO 42 . Additionally, a reciprocal blastn alignment using an e-value threshold of 10 −10 and a minimum alignment length of 100 nucleotides with at least 70% identity was used to compare the A. schaueriana transcriptome with other publicly available transcriptomes of congeneric species.
Comparative transcriptomics using RNA-sequencing. Tissue-specific count data were obtained from the number of reads uniquely mapped to each transcript of the non-redundant transcriptome using bowtie1 35 and normalized using edgeR 43 . Differentially expressed transcripts (DETs) between equatorial and subtropical tree tissue-specific samples were detected using the exact test for negative binomial distributions with an FDR < 0.05 (Benjamini-Hochberg). GO term enrichment analyses of the DETs were performed using GOseq 44 Table S2). We isolated DNA using the DNeasy Plant Mini Kit (QIAGEN) and NucleoSpin Plant II (Macherey Nagel) following the manufacturers' instructions. DNA quality and quantity were assessed using 1% agarose electrophoresis and the QuantiFluor dsDNA System with a Quantus fluorometer (Promega). Nextera-tagmented reductively-amplified DNA (nextRAD) libraries 45 were prepared and sequenced by SNPsaurus (SNPsaurus) on a HiSeq. 2500 (Illumina, Inc.) with 100-bp single-end chemistry. Genomic DNA fragmentation and short adapter ligation were performed with the Nextera reagent (Illumina, Inc.) followed by amplification with one of the primers matching the adapter and extending nine arbitrary nucleotides into the genomic DNA. Assembly, mapping and single nucleotide polymorphic loci (SNP) identification were performed using custom scripts (SNPsaurus), which created a reference catalogue of abundant reads across the combined sample set and mapped reads to this reference, allowing two mismatches and retaining biallelic loci present in at least 10% of the samples. We further filtered markers by allowing no more than 65% missing data, a Phred score > 30, 8x minimum coverage, only one SNP per locus and a minor allele frequency ≥0.05 using vcftools v.0.1.12b 46 . To reduce paralogy or low-quality genotype calls, we used a maximum read coverage of 56 (the average read depth times 1.5 standard deviation).
After excluding plants morphologically identified as A. schaueriana with genomic signs of hybridization with A. germinans (L.) L., we assessed the genetic structure considering all SNPs, using the discriminant analysis of principal components (DAPC) 47 and ADMIXTURE v.1.3.0 48 . For DAPC analyses, we considered the number of groups (K) varying from 1-50 and the Bayesian information criteria for inferring K. Additionally, we used the optim.a.score function to avoid overfitting during the discrimination steps. For the ADMIXTURE analyses, we performed three separate runs for K varying from 1-15 using the block-relaxation method for point estimation; computing was terminated when estimates increased by <0.0001, and the most likely K was determined by cross-validation. (2019) 9:19936 | https://doi.org/10.1038/s41598-019-56469-w www.nature.com/scientificreports www.nature.com/scientificreports/ We used two programs to minimize false-positive signs of natural selection: LOSITAN 49 , assuming an infinite allele model of mutation, using a confidence interval of 0.99, a false-discovery rate (FDR) of 0.1, the neutral mean FST and forcing of the mean FST options; and pcadapt 3.0.4 50 , which simultaneously identifies the population structure and the loci excessively related to this structure, using an FDR < 0.1.
Putative evidence of selection was considered only for SNP loci that were conservatively identified by pcadapt and five independent runs of LOSITAN to avoid false-positives 51 . As selection is presumably stronger in coding regions of the genome and there is no reference genome for the species, we used the de novo assembled transcriptome, characterized herein, as a reference to identify candidate loci within putative coding regions. We performed a reciprocal alignment between nextRAD sequences (75 bp) and longer expressed sequences (≈300-11600 bp) using blastn v.2.2.31 36 , with a threshold of at least 50 aligned nucleotides, a maximum of one mismatch and no gaps.

Results
comparative physiology using a common garden experiment. Seedlings from equatorial and subtropical sites diverged in key ecophysiological traits in the common garden experiment (Figs. 2-4). The inclination angle of the leaf lamina and the average size of individual leaves were smaller in equatorial than in subtropical plants, but total leaf area and specific leaf area did not differ between the groups (Fig. 2, Supplementary Information Table S3, Supplementary Information Fig. S3). Additionally, equatorial plants showed lower stomatal conductance and transpiration rates than did subtropical plants (Fig. 3). Subtropical plants accumulated more biomass in leaves and roots than did equatorial plants. However, the stem dry mass ratio (MR) (stem dry biomass/ plant dry biomass) was greater in equatorial plants, whereas the leaf MR (leaf dry biomass/plant dry biomass) was greater in subtropical plants (Fig. 2). The dry biomass accumulated in the stems and the root MR (root dry biomass/plant dry biomass) were indistinguishable between the groups (Supplementary Information Table S3). Unexpectedly, 63% of the equatorial plants flowered after six months of growth (Supplementary Information Fig. S3g). Since this was not observed in any subtropical plants, flowering plants were not used in the biomass allocation analyses.
Mean stem vessel diameter was smaller in the equatorial seedlings than in the subtropical seedlings, indicating enhanced hydraulic safety (Fig. 2, Supplementary Information Fig. S3h,i). However, the vessel density, vessel grouping index, vessel lumen area in sapwood and total hydraulic conductivity of the stems were not significantly different between the groups (P-value >0.05) (Supplementary Information Table S3). Plants from contrasting latitudes exhibited different stem epidermis pigmentation, with equatorial seedlings reflecting more long-wavelength of red light (635-700 nm) and less medium wavelength of green light (520-560 nm) than subtropical seedlings (Fig. 4). characterization of the A. schaueriana transcriptome. In the absence of a reference genome, we obtained the first reference transcriptome for A. schaueriana from leaves, stems and flowers of adult individuals under field conditions ( Supplementary Information Fig. S4, Table S1). Over 209 million paired-end 72-bp reads showing high quality were de novo assembled into a reference, non-redundant transcriptome containing 49,490 sequences, of which 30,227 (61%) were putative protein-coding sequences. Over 91.9% of these reads were mapped to a single transcript, indicating minimum redundancy and a wide representation of sequenced data (Supplementary Information Table S4). Moreover, 91.8% of universal plant orthologous genes were present in the reference transcriptome (Supplementary Information Table S4). Sequences with complete ORFs represented approximately 42% (12,796) of all putative protein-coding transcripts (Supplementary Information Table S5, Fig. S5). Most of these putative protein-coding sequences (94.33%) showed significant similarity (e-value <1e-10) to proteins in the Plant RefSeq and TAIR databases ( Supplementary Information Fig. S5c). More than 80% of protein-coding sequences matched proteins from Sesamum indicum or Erythranthe guttata, which, as A. schaueriana, belong to the order Lamiales ( Supplementary Information Fig. S6d) comparative transcriptomics between trees from contrasting latitudes. To identify the effects of environmental variation on gene expression and, thus, on phenotypes at contrasting source sites in the field, we compared expression levels in trees under field conditions as comparable as possible, considering tidal conditions and seasonal climates (Supplementary Information Table S1). As expected, we observed a consistent source-site pattern in overall transcript expression levels from leaves and stems (Supplementary Information Fig. S6a,b). However, floral transcript expression levels did not show a clear pattern among samples of the same origin ( Supplementary Information Fig. S6c), leading to the identification of only one DET. Thus, we did not include flowers in the subsequent analyses ( Supplementary Information Fig. S6f). Conversely, 1,837 and 904 transcripts showed significantly different (FDR < 0.05) relative abundances between equatorial and subtropical samples in leaves and stems, respectively ( Supplementary Information Fig. S6d,e). Among the total 2,741 DETs, 1,150 (41.91%) were putative protein-coding transcripts.
The assignment of transcripts to GO terms was possible for 25,184 (83.31%) of the 30,227 putative protein-coding sequences, allowing GO enrichment analyses. Analyses were conducted separately for leaves and stems and for each of the following two sets of DETs: one showing higher expression levels in equatorial than in subtropical samples (which we refer to as "DET-Eq") and the other showing higher abundance in subtropical than in equatorial samples (which are refer to as "DET-St"). We focused on the biological processes associated with key aspects of the responses to contrasting climate conditions between the equatorial and subtropical sites (Table 1, Supplementary Information Fig. S1). The enriched processes among the DET sets included photosynthesis;  Information Tables S7-S11, Supplementary Information Fig. S6i-l).
Photosynthesis. Among the DET-St set, we observed various putative genes participating in the biosynthesis of the photosynthetic apparatus, chlorophyll and photoreceptors; the function of electron transporters and chloroplast movement coordination. Contrastingly, the DET-Eq set showed enrichment of transcripts similar to  www.nature.com/scientificreports www.nature.com/scientificreports/ proteins required for disassembling the light-harvesting complex of photosystem II in thylakoid membranes 55 ( Supplementary Information Table S11).
Cell wall biosynthesis. Transcripts similar to those encoding 33 distinct proteins and transcription factors that play central roles in the biosynthesis and deposition of cell wall components, such as cellulose, hemicellulose, lignin and pectin, were identified in DET-Eq (Supplementary Information Table S10).
Response to UV. Both the DET-St and DET-Eq sets showed enrichment of functions related to responses to UV radiation; however, the transcript annotations differed between these sets. The DET-St set included putative UV-B protectants and genes involved in UV-B-induced antioxidant biosynthesis, such as plastid copper/zinc superoxide dismutases, photosystem II repair proteins, and L-ascorbic acid. In contrast, the DET-Eq set showed enrichment of transcripts associated with photo-oxidative damage reduction and the positive regulation of anthocyanin biosynthesis in response to UV. Antioxidants induced by UV irradiation 56 , such as putative iron superoxide dismutases and pyridoxine biosynthesis genes, were among the DET-Eq set ( Supplementary Information Table S11).
Response to temperature. In the DET-St set, we observed putative genes presenting critical roles in chloroplast protein translation during cold acclimation and that provide tolerance to chilling temperatures 57,58 . There included transcripts similar to the GLYCINE-RICH RNA-BINDING (RZ1A), which has chaperone activity during cold acclimation 59 , and to the cold-inducible ATP-dependent DNA HELICASE ATSGS1, required for DNA damage-repair 60 . Interestingly, DET-St included a putative AGAMOUS-LIKE 24 (AGL24) transcription factor that is involved in the vernalization-induced floral transition 61 . Contrastingly, various transcripts similar to heat shock-inducible chaperones and to ADENINE NUCLEOTIDE ALPHA HYDROLASE-LIKE (AT3G53990), involved in chaperone-mediated protein folding 62 , were among the DET-Eq set, potentially enhancing tolerance to heat in equatorial plants. Additionally, a transcript similar to the ETHYLENE-RESPONSIVE ELEMENT BINDING (RAP2.3), which confers resistance to heat and hydrogen peroxide 63 , was observed in this group (Supplementary Information Table S11).
Response to water stress. Transcripts associated with the response and tolerance to water deficits and with cellular ion homeostasis and osmotic adjustment were enriched among DET-Eq; for instance, a transcript similar to the ETHYLENE-RESPONSIVE TRANSCRIPTION FACTOR (RAP2.4), which regulates the expression of several drought-responsive genes, including aquaporins 64,65 . Accordingly, the putative aquaporin PLASMA MEMBRANE INTRINSIC (PIP1;4) 66 was found in this set. In DET-Eq, we observed putative genes participating in the synthesis of raffinose, an osmoprotective soluble trisaccharide 67 , as well as transcripts similar to osmosensitive ion channels belonging to the EARLY-RESPONSIVE TO DEHYDRATION STRESS FAMILY. Correspondingly, we observed an ion channel, SLAC1 HOMOLOGUE 3 (SLAH3), required for stomatal closure induced by drought stress 68 , and the putative NINE-CIS-EPOXYCAROTENOID DIOXYGENASE 3 (NCED3), which increases plant resistance to water deficits through the accumulation of abscisic acid (ABA), leading to stomatal closure. Possibly as a consequence of decreased stomatal conductance, a transcript similar to the photorespiration gene, ALANINE-GLYOXYLATE AMINOTRANSFERASE 2 HOMOLOG 3 (AT3G08860) 69 , was also observed in DET-Eq (Supplementary  Information Table S11).
We confirmed the results obtained from RNA-Seq data computational analyses by qRT-PCR using ten DETs from leaf samples ( Supplementary Information Fig. S7, Table S12, Supplementary Note).

Detection of SNPs with signs of selection.
RNA-Seq enabled the assembly of a reference transcriptome and the identification of biological processes influenced by the environmental divergence over the South American coastline. To complement these analyses, as one cannot disentangle the effects of plasticity and adaptive selection in the differential expression 70,71 , we searched for gene sequence variations among trees sampled along the Atlantic coast of South America (Fig. 1, Supplementary Information Table S2). After quality filtering of the sequenced data, we selected 77 individuals without evidence of interspecific hybridization with A. germinans for downstream analyses. We identified a set of 6,170 high-quality unlinked biallelic SNP loci with a minor allele frequency ≥0.05 and ≥8x coverage. The overall genetic structure of genome-wide SNPs corroborated the results of a previous study based on putatively neutral loci 22 , dividing the species into two populations: north and south of the northeast extremity of South America (Fig. 5).
We observed 122 loci showing considerable departures from neutral expectations of interpopulational divergence, as conservatively detected 72 by pcadapt and LOSITAN. Fifteen of these loci aligned to A. schaueriana transcripts that were similar to gene models in A. thaliana and S. indicum (Supplementary Information  Table S13), enabling screening for their potential functional roles. However, five reference proteins did not have informative functions described for the model plant, hindering functional inferences. Conversely, among the remaining annotated candidates, we found five putative genes involved in processes related to contrasting equatorial and subtropical environmental variables (Fig. 6). One candidate locus was detected in the putative transcription factor RAP2.4, which is induced in response to water and salt stress 64 and regulates the expression of aquaporins 65 . Two other candidates showed similarity to the transcription factors ZINC-FINGER PROTEIN 1 (ZFN1), involved in the regulation of the expression of several water stress genes 73 , and the HYPOXIA RESPONSE ATTENUATOR (HRA1), which is strongly induced in response to low oxygen levels 74 . A putative UDP GLUCOSYL TRANSFERASE, an enzyme that catalyzes the final step of anthocyanin biosynthesis wherein pigmentation is produced 75 , also showed evidence of positive selection. Additionally, one candidate locus was in a transcript similar to that of the TETRATRICOPEPTIDE REPEAT gene (AT2G20710, TPR), which might play a role in the biogenesis of the photosynthetic apparatus 76

Discussion
In this study, we integrated genomic and ecophysiological approaches to investigate the foundations of adaptive variations in a dominant tree from the Atlantic coast of South America. We tested the hypothesis that latitudinal variations in climate regimens shape genetic and phenotypic variation involved in resource-use strategies in widespread coastal trees, using the mangrove Avicennia schaueriana as a model. Overall, our results supported this hypothesis. Using a common garden experiment, we detected differences between plants from contrasting latitudes in key ecophysiological traits involved in carbon and water balances, indicating the heritability of trait divergence. Accordingly, transcriptomic changes between plants sampled in contrasting latitudes showed enrichment in processes associated with central aspects of the environmental divergence across latitudes, such as responses to temperature, solar radiation, water deficit, photosynthesis and cell wall biosynthesis. The relevance of these biological processes in the field was corroborated by the identification of SNP loci with signs of selection present within putative coding transcripts similar to genes involved in responses to osmotic stress, accumulation of anthocyanins and photosynthesis, unveiling a genetic basis for environmental adaptation. Figure 6 summarizes the integration of three independent approaches that converge to suggest population adaptation to contrasting environments over the species range.

Evidence of a conservative resource-use strategy in the equatorial population of A. schaueriana.
Traits exhibited by equatorial plants compared to subtropical plants during the common garden experiment, such as a smaller leaf size and angle, higher levels of red light-reflecting pigments, narrower vessels and lower rates of stomatal conductance, limit carbon acquisition 77 and may have imposed constraints on carbon gain in equatorial plants. Accordingly, these plants also accumulated less biomass than subtropical plants. Despite causing limitations in growth, these characteristics allow plants to maintain suitable leaf temperatures for photosynthesis while reducing the risk of cavitation, UV exposure and water loss through the minimization of evaporative cooling 77 . Given the relevance of these traits to water and carbon balances, especially in high-salinity environments, the prevalence of more conservative resource-use traits among equatorial samples in the common garden experiment suggested selection favoring drought tolerance, likely being advantageous during hot and dry seasons in equatorial wetlands (August-December). These seasons frequently present a combination of stressful conditions, including high UV exposure, highly fluctuating soil salinity resulting from a wide daily tidal range, and a high vapor pressure deficit (VPD) resulting from high temperature (>30 °C) and an air humidity below 60%. Accordingly, equatorial plants also showed lower transpiration rates than subtropical plants in the common garden. Additionally, 63% of the six-month-old equatorial plants flowered from July-August in the common garden, whereas subtropical plants did not flower. This period is consistent with phenological observations reported for A. schaueriana in both equatorial 78 and southern subtropical forests 79 , and could indicate a genetic basis for the observed variation. Early flowering is a phenotype with complex genetic components, rarely studied in non-model organisms but renowned as an adaptive mechanism for maximizing the chances of reproduction under stressful environments 80 .
In their native environment, equatorial plants showed increased expression levels of putative heat-shock proteins, drought-induced ion transporters, aquaporins and genes that play central roles in stomatal closure and photorespiration compared to subtropical plants, likely contributing to heat and drought tolerance. These findings provided multiple lines of evidence of heat and water stress responses. Higher expression of aquaporins and genes involved in the accumulation of organic solutes may contribute to lowering the cellular water potential while improving drought tolerance in equatorial plants compared to subtropical plants 67,81 . Additionally, higher expression of several transcripts associated with secondary cell wall biosynthesis and thickening may enhance the rigidity of cells and reduce the risk of collapse during dehydration-rehydration cycles in equatorial trees 82 . Moreover, equatorial plants showed lower expression of several putative photosynthesis genes than subtropical plants. In response to drought, high-light and heat stress plants minimize the photooxidative damage by reducing photosynthetic activity via repression of photosynthesis genes [83][84][85][86] . Remarkably, evidence of selection was detected in two putative transcription factors, RAP2.4 and ZFN1, which play key roles in the regulation of osmotic stress-response 64,73 . These findings support the hypothesis that dry seasons marked by low rainfall and high VPD, which are caused by the combination of high temperatures and low air humidity 87 in the equatorial region, induce a pivotal selective pressure for coastal trees populations.

Evidence of an acquisitive resource-use strategy in the subtropical population of A. schaueriana.
Subtropical plants showed higher stomatal conductance and transpiration rates, higher levels of green light-reflecting pigments, larger leaf area, wider leaf lamina angle and larger xylem vessel diameter than equatorial plants in the common garden experiment. These characteristics enhance light energy absorbance and carbon www.nature.com/scientificreports www.nature.com/scientificreports/ acquisition, at the expense of a higher cavitation risk 88,89 . Conversely, this apparent riskier strategy may compensate for seasonal declines in growth resulting from decreasing temperature, photoperiod and light quality at higher latitudes 90 . Although low temperatures reduce enzymatic activity 91 and, thus, plant growth, the severity of low-temperature stress in the southernmost subtropical mangrove forests of the Americas is likely insufficient to favor the selection of freezing-tolerant adaptations, in contrast to results reported for mangroves at their northernmost edge on the Atlantic coast of the Americas 92 . At the southernmost limit of American mangrove forests, the minimum air temperature does not drop below 0 °C and remains within the expected mangrove physiological thresholds 93 . Additionally, the small population size of A. schaueriana at this location 27 and the arrival of maladapted propagules from northerly populations likely reduce the potential strength of selection favoring cold-adaptive traits.
Corroborating the observed differences in morphophysiological traits, we found divergence at the molecular level that may also be related to the increasing amplitude of variation in photoperiod and light quality towards high latitudes. Plants optimize the use of light energy and adjust photosynthetic activity through the regulation of light-harvesting and photosystem-component genes 83 . Thus, the higher expression levels of transcripts associated with photosynthesis in subtropical than in equatorial plants may have facilitated the absorption of light energy in subtropical plants during winter. Although solar irradiance levels were indistinguishable between sampling sites at the time of sampling, transcriptomic changes in putative UV-inducible antioxidant and photodamage repair genes suggest the use of distinct strategies to respond to differential seasonality in photoperiod and solar radiation levels between subtropical and equatorial latitudes. Notably, we observed signatures of selection in two transcripts, one showing similarity to a UDP-GLUCOSYL TRANSFERASE, a key enzyme in the anthocyanin biosynthesis pathway, and the other in a transcript similar to a TPR gene required for chlorophyll biosynthesis and for the biogenesis of the photosynthetic apparatus 76 . These results imply that solar radiation, in addition to VPD, may act as an environmental factor driving selection in A. schaueriana.
Although soil oxygen availability affects plant growth in intertidal areas, we did not focus our experiments on the relevance of hypoxia in shaping adaptive divergence in coastal wetlands 94 . However, evidence of selection was detected in a transcription factor highly induced by oxygen deprivation (HRA1) 74 . Additionally, the HRA1 putative homolog also showed 1.75-fold higher expression in subtropical leaves relative to that in equatorial leaves, even though sampling was conducted at low tide at both sites. As tidal amplitude decreases with increasing latitude along the Atlantic coast of South America 95 (Fig. 6), trees are exposed to increasing anoxia conditions southwards. These findings suggest that subtropical populations may have better stress preparedness for hypoxia than equatorial populations. climate change and conservation implications. Our results provide compelling evidence that adaptations to contrasting seasonality in freshwater availability and solar radiation are involved in the organization of diversity of the dominant coastal tree A. schaueriana. The functional divergence described herein might differentially affect the sensitivity of northerly and southerly populations to a rapidly changing climate. It has been suggested that the southernmost mangrove forests in the Americas could expand poleward in the near future 96 . We expect that the observed acquisitive resource-use characteristics may indeed favor subtropical tree growth under increased CO 2 concentrations, temperatures and rainfall, as predicted for this region by 2100 97 . However, great landward relocation of subtropical populations will be necessary for their persistence in the face of a rising sea level 98 due to the narrowing of intertidal zones towards mid-latitude areas on the Atlantic coast of South America. Even though subtropical populations apparently show better preparedness for hypoxia than equatorial populations, this scenario is aggravated by the dense coastal occupation by urban and agricultural areas, which may preclude the landward migration of subtropical mangroves in South America. Contrastingly, equatorial populations frequently have wider plains that are potentially available for expansion and lower human occupation, but might be threatened by reduced rainfall and increased VPD during more intense El Niño-Southern Oscillation events 97 . Increased temperature will stimulate both respiration 99 and photorespiration 100 and may eventually offset benefits in carbon acquisition caused by increased CO 2 concentrations 101 . The critical temperature threshold for photosynthesis may be overcome more frequently, possibly reducing carbon assimilation and productivity 102 . With a conservative resource-use strategy, further limitations in net carbon assimilation could, in extreme events, lead to biomass loss or tree mortality triggered by carbon starvation 103 . This scenario could be especially severe to the semiarid northeast South American coast, particularly in the face of intense human use of adjacent plains and decreasing trends in terrestrial freshwater availability and coastal freshwater input 104,105 .
Our results corroborate previous studies of marine species from the Northern Atlantic that associated limited north-south population connectivity with adaptation to latitudinal environmental dissimilarity [19][20][21] . As widespread trees distributed along the Atlantic coast of South America show an overlapping north-south organization of genetic variation [22][23][24] , we hypothesize that they may also show divergent adaptations to heterogeneous resource availability over their distribution ranges. Conservation efforts of coastal ecosystems should focus on resilient, habitat-forming species that give shelter and act as climate rescuers for several stress-sensitive species 106 . We recommend that populations of coastal trees occurring north and south from the northeast extremity of South America should warrant attention as distinct conservation management units 107 for the long-term persistence of coastal biodiversity and the ecosystem services they provide.

conclusions
Studies on local adaptations in tropical trees are frequently limited by difficulties in implementing traditional approaches, such as reciprocal transplants, due to their long generation time and the lack of basic biological information, such as knowledge of their evolutionary history, reproductive mode and phenological patterns. Despite limitations, we provided integrated ecophysiological and genomic evidence that supports the hypothesis that seasonal environmental heterogeneity in solar radiation and freshwater availability may shape the variation of traits associated with resource-use in A. schaueriana, a widespread tree along the Atlantic coastline of South America. The non-neutral trait diversity revealed by the interdisciplinary approach used here provides a perspective into the molecular and phenotypic scales at which environmental selection may shape functional variation of this dominant species, on a continental scale. Integrating high-throughput sequencing and ecophysiological approaches, has shown to be a promising strategy in the study of adaptation in non-model, long-lived species. To provide more realistic predictions of how dominant coastal trees may respond to current global climate changes we encourage the development of further studies accounting for phenotypically 14 and genetically 108 informed distribution modeling. Since freshwater availability has been decreasing in coastal areas worldwide 104 , strongly compromising the productivity of coastal plants 105 , such studies should focus on biological variation involved in the balance between carbon gain and water loss. This knowledge can improve predictions on the future of ecosystems they form and generate key information for forest conservation and management efforts 14,109 . Without the realization and dissemination of such studies, the success of conservation plans for tropical forests and their potential to mitigate climate change will likely be seriously compromised.

Data availability
Expression data and sequences that support the findings have been deposited in GenBank with the primary accession code GSE116060. A Variant Call Format file and its complementary file, both required for all of the genome-wide SNP diversity analyses are provided as Supplementary Files (Supplementary Datasets 1, 2).