Lígia M Ayres-Ostrock

Myriam Valero

Stéphane Mauger

Mariana C Oliveira

Estela M Plastino
email: emplasti@usp.br

Marie-Laure Guillemin

Dual influence of terrestrial and marine historical processes on the phylogeography of the Brazilian intertidal red alga Gracilaria caudata 1

Keywords: COI mtDNA, Gracilaria caudata, microsatellite, phylogeography, population genetics, refugia, Rhodophyta, South-western Atlantic Coast Linkage disequilibrium

In this study, we explored how both past terrestrial and marine climate changes have interacted to shape the phylogeographic patterns of the intertidal red seaweed Gracilaria caudata, an economically important species exploited for agar production in the Brazilian north-east. Seven sites were sampled along the north-east tropical and south-east sub-tropical Brazilian coast. The genetic diversity and structure of G. caudata was inferred using a combination of mitochondrial (COI and cox2-3), chloroplast (rbcL) and 15 nuclear microsatellite markers. A remarkable congruence between nuclear, mitochondrial and chloroplast data revealed clear separation between the north-east (from 03°S to 08°S) and the south-east (from 20°S to 23°S) coast of Brazil. These two clades differ in their demographic histories, with signatures of recent demographic expansions in the north-east and divergent populations in the south-east, suggesting the maintenance of several refugia during the last glacial maximum due to sea-level rise and fall. The Bahia region (around 12 o S) occupies an intermediate position between both clades. Microsatellites and mtDNA markers showed additional levels of genetic structure within each sampled site located south of Bahia. The separation between the two main groups in G. caudata is likely recent, probably occurring during the Quaternary glacial cycles. The genetic breaks are concordant with (1) those separating terrestrial refugia, (2) major river outflows and (3) frontiers between tropical and subtropical regions. Taken together with previously published eco-physiological studies that showed differences in the physiological performance of the strains from distinct locations, these results suggest that the divergent clades in G. caudata correspond to distinct ecotypes in the process of incipient speciation and thus should be considered for the management policy of this commercially important species.

INTRODUCTION

Historical processes have likely influenced hydrographic-climatic patterns, in addition to modifying coastline profiles and marine species' range dynamics [START_REF] Haq | Chronology of fluctuating sea levels since the Triassic[END_REF][START_REF] Hewitt | Some genetic consequences of ice ages, and their role in divergence and speciation[END_REF], 2000, 2004). Likewise, these processes may have shaped genetic diversity: the genetic structure of contemporary populations is the result of both short-term ecological processes and long-term evolution governed in part by environmental changes [START_REF] Benzie | Genetic structure of marine organisms and SE Asian biogeography[END_REF]. Phylogeographic approaches, generally using uniparentally inherited molecular markers, are powerful tools for inferring past range contraction and expansion and to establish evolutionary origin of genetic characteristics of presentday populations [START_REF] Avise | Phylogeography. The history and formation of species[END_REF]. The late Pleistocene glaciations are recognised as important factors that favoured the divergence of populations and ultimately speciation via repeated isolation in allopatric refuges [START_REF] Hewitt | Some genetic consequences of ice ages, and their role in divergence and speciation[END_REF], 2004[START_REF] Maggs | Evaluating signatures of glacial refugia for North Atlantic benthic marine taxa[END_REF][START_REF] Neiva | Climate oscillations, range shifts and phylogeographic patterns of North Atlantic Fucaceae[END_REF][START_REF] Bringloe | Mitochondrial DNA sequence data reveal the origins of postglacial marine macroalgal flora in the Northwest Atlantic[END_REF]. After the Last Glacial Maximum (LGM), it is assumed that climatic amelioration allowed species to expand their range poleward from the populations located at low latitude margins [START_REF] Davis | Range shifts and adaptive responses to Quaternary climate change[END_REF]. Thus, current patterns of genetic diversity have been largely determined by population responses at the margins of species' distribution ranges [START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF]. Longterm persistent populations, also known as 'rear-edge' populations, are usually restricted to areas where survival was possible under glacial maximum conditions. These low-latitude refugia show generally high genetic variation despite bottlenecks in these isolated populations during glacial periods [START_REF] Hewitt | Some genetic consequences of ice ages, and their role in divergence and speciation[END_REF], 2004, Assis et al. 2016). On the other hand, 'leading-edge' populations result from rare long-distance dispersal events followed by exponential population growth and are generally characterised by highly reduced genetic diversity [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF].

In tropical and sub-tropical zones, the effect of glacial-interglacial cycles was less dramatic than in temperate ones and changes in species distribution, mostly related to a cooler and more arid climate during glacial periods, generally do not reflect the classical contraction and expansion reported at higher latitudes [START_REF] Lessa | Genetic footprints of demographic expansion in North America, but not Amazonia, during the Late Quaternary[END_REF]. In Brazil, phylogeographic studies have mainly focused on terrestrial and freshwater habitats (e.g. [START_REF] Costa | Phylogeography and systematic notes on two species of gracile mouse opossums, genus Gracilinanus (Marsupialia: Didelphidae) from Brazil[END_REF][START_REF] Rocha | Patterns of distribution and processes of speciation in Brazilian reef fishes[END_REF][START_REF] Carnaval | Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest[END_REF], whereas marine organisms have received little attention [START_REF] Pinheiro | Island biogeography of marine organisms[END_REF][START_REF] Nauer | Phylogeography of the Hypnea musciformis species complex (Gigartinales, Rhodophyta) with the recognition of cryptic species in the western Atlantic Ocean[END_REF]. Nonetheless, studies on marine population genetics are crucial to assure a better sustainable management and future economic stability of these natural resources [START_REF] Palumbi | Population genetics, demographic connectivity, and the design of marine reserves[END_REF][START_REF] Waples | Integrating genetic data into management of marine resources: how can we do it better?[END_REF][START_REF] Couceiro | Management and conservation of the kelp species Laminaria digitata: using genetic tools to explore the potential exporting role of the MPA "Parc naturel marin d'Iroise[END_REF].

The Brazilian coast extends for about 8000 km and stretches across 38 degrees of latitude (from 5°16'20"N to 33°44'32"S, Figure 1). This long coastline provides a variety of ecological conditions resulting in very different abiotic environments (Gomes da Silva et al. 2016).

Phylogeographic patterns in marine fishes and invertebrates are generally more related to the impact of the strong currents affecting the region and acting as barriers to gene flow (e.g., the bifurcation into two branches of the South Equatorial Current at latitudes 10 o -15 o S; see Figure 1; [START_REF] Santos | Population genetic structuring of the king weakfish, Macrodon ancylodon (Sciaenidae), in Atlantic coastal waters of South America: deep genetic divergence without morphological change[END_REF][START_REF] Da Silva Cortinhas | Genetic structuring among silverside fish (Atherinella brasiliensis) populations from different Brazilian regions[END_REF][START_REF] Dos Santos Freitas | The mitochondrial control region reveals genetic structure in southern kingcroaker populations on the coast of the Southwestern Atlantic[END_REF][START_REF] Peluso | Contemporary and historical oceanographic processes explain genetic connectivity in a Southwestern Atlantic coral[END_REF] than to demographic processes linked to Pleistocene climatic changes. Interestingly, intertidal species such as some seaweeds located at the interface of the terrestrial and marine realms may have been affected by past marine and terrestrial climatic fluctuations [START_REF] Neiva | Species distribution models and mitochondrial DNA phylogeography suggest an extensive biogeographical shift in the high-intertidal seaweed Pelvetia canaliculata[END_REF][START_REF] Cardoso | Phylogeography of the sand dune ant Mycetophylax simplex along the Brazilian Atlantic Forest coast: remarkably low mtDNA diversity and shallow population structure[END_REF]. The existence of three main refugia for the rainforest located along the Brazilian Atlantic coast have been postulated (i.e. in the region of Pernambuco-PE, Bahia-BA, and São Paulo-SP, Figure 1) for which long term persistence of genetically isolated lineages have been associated (Carnaval andMoritz 2008, Carnaval et al. 2009). Additionally, along the Brazilian coast, the Holocene sea-level history suggests regional variations alternating periods of mean sea-level rise and fall. It may be assumed that these oscillations would have had a strong influence on the distribution and diversity of coastal species [START_REF] Cardoso | Phylogeography of the sand dune ant Mycetophylax simplex along the Brazilian Atlantic Forest coast: remarkably low mtDNA diversity and shallow population structure[END_REF][START_REF] Leite | Neotropical forest expansion during the last glacial period challenges refuge hypothesis[END_REF].

Phylogeographic studies on seaweeds have increased significantly over the past two decades (for review see [START_REF] Hu | Seaweed phylogeography from 1994 to 2014: an overview[END_REF]. Most of these studies have been carried out in the Northern Hemisphere and only a few have focused on South America [START_REF] Tellier | Phylogeographic analyses of the 30° S south-east Pacific biogeographic transition zone establish the occurrence of a sharp genetic discontinuity in the kelp Lessonia nigrescens: Vicariance or parapatry?[END_REF], 2011[START_REF] Guillemin | Phylogeography of seaweeds in the South East Pacific: complex evolutionary processes along a latitudinal gradient[END_REF], 2018[START_REF] Montecinos | Species replacement along a linear coastal habitat: phylogeography and speciation in the red alga Mazzaella laminarioides along the south east pacific[END_REF]. Recently, [START_REF] Nauer | Description of Hypnea pseudomusciformis sp. nov., a new species based on molecular and morphological analyses, in the context of the H. musciformis complex (Gigartinales, Rhodophyta)[END_REF][START_REF] Nauer | Phylogeography of the Hypnea musciformis species complex (Gigartinales, Rhodophyta) with the recognition of cryptic species in the western Atlantic Ocean[END_REF] highlighted the existence of intraspecific diversity within the intertidal red alga Hypnea pseudomusciformis along the Brazilian coast, with divergence between samples from the north-central part of the distribution and the ones located on the more south-eastern part of the coast. Here, we report on the phylogeographic patterns of the intertidal red alga Gracilaria caudata. The order Gracilariales a taxonomically challenging group of marine benthic red algae, includes species of major ecological and economical importance (worldwide invasive species: [START_REF] Bellorin | Gracilaria vermiculophylla: a western Pacific species of Gracilariaceae (Rhodophyta) first recorded from the eastern Pacific[END_REF][START_REF] Krueger-Hadfield | Invasion of novel habitats uncouples haplo-diplontic life cycles[END_REF]; harvested and cultivated agarophyte species: [START_REF] Valero | Perspectives on domestication research for sustainable seaweed aquaculture[END_REF]. . Morphological species delineation is particularly difficult due to the relatively small number of diagnostic characters and their great phenotypic plasticity (Gurguel et al. 2004). Recent molecular-based analyses revealed the occurrence of numerous cryptic species in this group [START_REF] Cohen | Molecular identification of two sibling species under the name Gracilaria chilensis (Rhodophyta, Gracilariales)[END_REF][START_REF] Guillemin | Genetic variation in wild and cultivated populations of the haploid-diploid red alga Gracilaria chilensis: how farming practices favor asexual reproduction and heterozygosity[END_REF][START_REF] Destombe | Delineation of two sibling red algal species, Gracilaria gracilis and Gracilaria dura (Gracilariales, Rhodophyta), using multiple DNA markers: resurrection of the species G. dura previously described in the Northern Atlantic 200 years ago[END_REF][START_REF] Lyra | Delimitating cryptic species in the Gracilaria domingensis complex (Gracilariaceae, Rhodophyta) using molecular and morphological data[END_REF]. The difficulty in distinguishing species also explains the current debate about the classification of the Gracilariales (see [START_REF] Gurgel | Systematics of the Gracilariales (Rhodophyta) including new subfamilies, tribes, subgenera, and two new genera, Agarophyton gen. nov. and Crassa gen[END_REF][START_REF] Iha | Organellar genomics: a useful tool to study evolutionary relationships and molecular evolution in Gracilariaceae (Rhodophyta)[END_REF]. Gracilaria caudata is a haploid-diploid intertidal rocky shore species characterized by forming dense beds. The species is widespread and occurs from southern Florida (USA-27°N) to the southern part of the Brazilian coast (Santa Catarina State, SC-27 o S; [START_REF] Plastino | Gracilaria caudata J. Agardh (Gracilariales, Rhodophyta)restoring an old name for a common western Atlantic alga[END_REF]. This species is an important source of agar production in Brazil [START_REF] Hayashi | Cultivation of red seaweeds: a Latin American perspective[END_REF]. Over the last 50 years, intensive and uncontrolled exploitation of natural beds has led to a significant decline of native populations [START_REF] Hayashi | Cultivation of red seaweeds: a Latin American perspective[END_REF]. Today, artisanal mariculture of G. caudata in the north-east State of Ceará, Brazil offers an important source of income for the local fishing community [START_REF] Costa | The Gracilariaceae Germplasm Bank of the University of São Paulo, Brazil-a DNA barcoding approach[END_REF]. Over the past two decades, the taxonomy, physiology, ecology, life history and mariculture of this species of economic interest have been intensively studied [START_REF] Plastino | Gracilaria caudata J. Agardh (Gracilariales, Rhodophyta)restoring an old name for a common western Atlantic alga[END_REF][START_REF] Costa | The Gracilariaceae Germplasm Bank of the University of São Paulo, Brazil-a DNA barcoding approach[END_REF][START_REF] Araújo | Intraspecific variation in Gracilaria caudata (Gracilariales, Rhodophyta): growth, pigment content, and photosynthesis[END_REF][START_REF] Faria | Physiological assessment of the mariculture potential of a Gracilaria caudata (Gracilariales, Rhodophyta) variant[END_REF][START_REF] Faria | Ecotypes of Gracilaria caudata (Gracilariales, Rhodophyta): physiological and morphological approaches considering life history phases[END_REF][START_REF] Trigueiro | Cultivation of three color strains of Gracilaria domingensis in an integrated organic system[END_REF]. Recently, a comparative study of individuals from three distant geographical areas along the Brazilian coast done under laboratory conditions, showed clear differences in the physiological performance of G. caudata strains among regions with strains from the north-east state of Ceará (CE) presenting higher growth rates and better photosynthetic performances than the ones sampled in the south (BA and SP) [START_REF] Faria | Ecotypes of Gracilaria caudata (Gracilariales, Rhodophyta): physiological and morphological approaches considering life history phases[END_REF]. Moreover, differences in growth rate and sensitivity to UVB radiation were observed between strains from CE and SP under controlled laboratory conditions [START_REF] Araújo | Intraspecific variation in Gracilaria caudata (Gracilariales, Rhodophyta): growth, pigment content, and photosynthesis[END_REF]). These differences suggest the existence of intraspecific diversity, supporting the hypothesis of ecotypic differentiation within this species [START_REF] Araújo | Intraspecific variation in Gracilaria caudata (Gracilariales, Rhodophyta): growth, pigment content, and photosynthesis[END_REF][START_REF] Faria | Ecotypes of Gracilaria caudata (Gracilariales, Rhodophyta): physiological and morphological approaches considering life history phases[END_REF], and raising questions on gene flow patterns and possible reproductive isolation along the extensive range of G. caudata. These results are also critical for developing effective management strategies in priority areas for conservation of coastal and marine biodiversity in Brazil [START_REF] Prates | Coastal and marine conservation priorities in Brazil. Priorities for coastal and marine conservation in South America[END_REF].

Using mitochondrial and chloroplast DNA sequences and nuclear microsatellite markers, we examine the population genetic structure of G. caudata throughout its whole distribution range in Brazil. We use this combined information to test for 1) the existence of differentiated genetic groups that echo the ecotypic differences revealed by eco-physiological experiments and 2) the processes responsible for population structuring. We hypothesised that both past isolation in refugia and oceanic barriers to gene flow have affected genetic structuring of G. caudata and that the main barriers to gene flow coincide with previously reported biogeographical breaks in the South-West Atlantic and/or with strong marine circulation patterns.

MATERIAL AND METHODS

Study Species

Gracilaria caudata is a commonly encountered seaweed in Brazil. The species grows on rocky substrate, often partially buried in sand and forming dense beds. It occurs mostly in protected bays and turbid waters, extending from the intertidal to the subtidal fringe [START_REF] Plastino | Gracilaria caudata J. Agardh (Gracilariales, Rhodophyta)restoring an old name for a common western Atlantic alga[END_REF]. In Gracilaria, the tetrasporophytic, female gametophytic and male gametophytic individuals are isomorphic. The thallus is an erect system of cylindrical branches that grow from the holdfast, fixing the individual to the substrate. Male and tetrasporophytic individuals can be distinguished by their reproductive structures, which are readily visible under a dissecting microscope, while female individuals are recognized by the presence of cystocarps (when fertilized) and can be detected by direct observation. In Gracilaria, as in most red algae, none of the propagules released in the seawater during the sexual life cycle (i.e, spermatia, haploid and diploids spores) are motile (Kain and Destombe 1995) and even if thalli are detached and found washed up on beach, they are not buoyant. As expected, gene flow has been estimated to be generally restricted to a few meter or kilometer in these algae [START_REF] Engel | Performance of non-motile male gametes in the sea: analysis of paternity and fertilization success in a natural population of a red seaweed, Gracilaria gracilis[END_REF], 2004[START_REF] Guillemin | Genetic variation in wild and cultivated populations of the haploid-diploid red alga Gracilaria chilensis: how farming practices favor asexual reproduction and heterozygosity[END_REF]; but in [START_REF] Guillemin | Tracing the trans-Pacific evolutionary history of a domesticated seaweed (Gracilaria chilensis) with archaeological and genetic data[END_REF] an event of rare long-distance dispersal caused by trans-oceanic colonization of non-buoyant algae species in association with rafting seaweeds is suggested for G. chilensis.

Sample collection

A total of 735 individuals of G. caudata were randomly collected from seven sites (Figure 1) covering most of the species' distribution range along the Brazilian coast including priority areas for conservation of coastal and marine biodiversity in Brazil [START_REF] Prates | Coastal and marine conservation priorities in Brazil. Priorities for coastal and marine conservation in South America[END_REF]. Two sampling strategies were adopted depending on algal density and spatial distribution. When possible (i.e. in CE, RN, PB, PE and BA), three transects of 20 meters in length were done perpendicularly to the shoreline during the low tide. For each transect, the apical branches of 35 individuals were sampled, resulting in 105 per site site. At sites with low densities of G. caudata (i.e. ES and SP), the sampling was done randomly and at least 50 individuals were collected. Each sample corresponds to an individual holdfast. All specimens were collected in between high and low intertidal areas uncovered at low tide.

Life stage and sex of individuals were determined by reproductive structures (by eyes for fertilized females or under a stereoscopic microscope for males and tetrasporophytes).Only diploid individuals (a total of 411 samples, i.e. more than 50 individuals per site) were preserved for genotyping. For the phylogeographic study, a subsampled of 18 tetrasporophytic individuals per site were sequenced.

DNA extraction and PCR amplification of mitochondrial, chloroplast and nuclear microsatellite markers DNA extractions followed the procedures described in Ayres-Ostrock et al. (2016). PCR amplifications of the mitochondrial cytochrome c oxidase I gene (COI) and the intergenic spacer located between the cytochrome oxidase subunits 2 and 3 genes (cox2-3) were performed following the protocols described in [START_REF] Saunders | Applying DNA barcoding to red macroalgae: a preliminary appraisal holds promise for future applications[END_REF] and [START_REF] Zuccarello | A mitochondrial marker for red algal intraspecific relationships[END_REF]. To supplement the mitochondrial marker data set, some individuals were sequenced for the chloroplast gene for the large subunit of ribulose-l, 5-bisphosphate carboxylase/oxygenase (rbcL). PCR amplifications were carried out following the protocol described in [START_REF] Freshwater | A gene phylogeny of the red algae (Rhodophyta) based on plastid rbcL[END_REF].

All PCR products were purified using Illustra™ GFX™ MicroSpin™ columns (GE Healthcare, Chicago, USA) and sequenced using the BigDye™ Terminator v3.1 Cycle Sequencing kit (Thermo Fisher Scientific, Waltham, USA), with the forward and the reverse amplification primers. PCR amplifications and genotyping of the 15 microsatellite loci were performed according to Ayres-Ostrock et al. (2016).

Mitochondrial and chloroplast DNA markers: sequence alignment, phylogenetic reconstructions, diversity and network analyses

Sequences were edited and aligned using ClustalW/Bioedit v 7 [START_REF] Hall | BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT[END_REF] and Mega v 5 [START_REF] Tamura | MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods[END_REF]) and deposited in GenBank under accession numbers MF995393-MF995549 for COI, MG452409-MG452552 for cox2-3 and MF995550-MF995561 for rbcL.

Maximum likelihood (ML) phylogenetic trees were constructed independently for rbcL and for the concatenated sequences of the two mitochondrial markers (COI and cox2-3) with Iq-Tree software using 2000 iterations for the optimal-tree search and as bootstrap pseudo-replicates [START_REF] Minh | Ultrafast approximation for phylogenetic bootstrap[END_REF][START_REF] Nguyen | IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies[END_REF][START_REF] Kalyaanamoorthy | ModelFinder: fast model selection for accurate phylogenetic estimates[END_REF]. Phylogenetic trees were visualised in Figtree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). For rbcL, four sequences of G. secunda and eight of G. caudata from the USA, Mexico, Venezuela and Brazil were downloaded from GenBank and added to our data set. For both trees, G. birdiae was used as outgroup. Moreover, a haplotype network was reconstructed for the concatenated COI+cox2-3 dataset using the median-joining algorithm implemented in Network v 5.0.0.1 [START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF].

Pairwise values of F ST among locations were calculated and their significance was tested using 1000 permutations. Moreover, a nested analysis of molecular variance (AMOVA) was implemented to test for the partition of genetic variance within sites, between sites within haplogroups and among haplogroups [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data[END_REF]). The isolation-by-distance model (Mantel test) was used to test for a relationship between genetic distance (F ST ) and geographic distance (in km), and its significance was tested using 1000 permutations.

For each site, four diversity indices were calculated for the COI and cox2-3 markers using Arlequin: the number of haplotypes (h), the number of polymorphic sites (S), gene diversity (H) (Nei, 1987) and nucleotide diversity (π) [START_REF] Nei | Mathematical model for studying genetic variation in terms of restriction endonucleases[END_REF]. A Mann-Whitney U test was used to evaluate differences in genetic diversity among regions (north-east vs. south-east) using Statistica 10 software.

For the 15 microsatellite loci, allele size was scored manually according to Ayres-Ostrock et al. (2016). Initially, the frequency of null alleles was estimated using MICRO-CHECKER software (van Oosterhout et al. 2004). The frequency of different multi-locus genotypes was calculated using the Monte Carlo procedure implemented in GenClone 2.0 (Arnaud-Haond and Belkhir 2007). Genotype diversity was calculated as R= (G-1)/(N-1), where G is the number of distinct genotypes identified and N is the number of individuals [START_REF] Dorken | Severely reduced sexual reproduction in northern populations of a clonal plant, Decodon verticillatus (Lythraceae)[END_REF]. Linkage disequilibrium (LD) was assessed using the association index (Brown et al. 1980, modified by Agapow and[START_REF] Agapow | Indices of multilocus linkage disequilibrium[END_REF]) and r d computed using Multilocus v1.3b [START_REF] Agapow | Indices of multilocus linkage disequilibrium[END_REF]. Their significance was tested using 1000 permutations [START_REF] Burt | Molecular markers reveal cryptic sex in the human pathogen Coccidioides immitis[END_REF].

Single and multi-locus estimates of gene diversity were calculated as mean number of alleles per locus (N a ), expected heterozygosity (H e , sensu [START_REF] Nei | Estimation of average heterozygosity and genetic distance from a small number of individuals[END_REF]) and observed heterozygosity (H o ) using Genetix 4.05 software (Belkhir et al. 1996(Belkhir et al. -2004)). Single-and multi-locus estimates of deviation from random mating (F IS ) were calculated according to [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] and significant was tested by running 1000 permutations using Genetix.

Genetix was also used to calculate the genetic differentiation (F ST ) between sites, and to compute the Mantel test, which evaluates the existence of a correlation between genetic and geographic distances. The possibility that G. caudata populations underwent a recent bottleneck event was evaluated using Bottleneck software v1.2.02 [START_REF] Piry | Computer note. BOTTLENECK: a computer program for detecting recent reductions in the effective size using allele frequency data[END_REF]. Moreover, an AMOVA was implemented to test for the partition of genetic variance within sites, among sites within haplogroups and among haplogroups (haplogroups as defined with mitochondrial sequence dataset; [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data[END_REF].

The Bayesian clustering method as implemented in Structure v2.3.4 software [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF] was used to determine the existence of genetic groups within G. caudata populations categorising them into K sub-populations. A range of clusters (K), from one to eight was tested. Each run, replicated 10 times, consisted of 400,000 iterations after a 'burn-in' of 200,000 iterations. K was determined by the method developed by [START_REF] Evanno | Detecting the number of clusters of individuals using the software structure: a simulation study[END_REF]. Structure results were combined using the greedy algorithm with 100,000 random input orders in Clummp [START_REF] Jakobsson | CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure[END_REF] and visualised with the Distruct programme [START_REF] Rosenberg | DISTRUCT: a program for the graphical display of population structure[END_REF]. The discriminant analysis of principal components (DAPC, [START_REF] Jombart | Discriminant analysis of principal components: a new method for the analysis of genetically structured populations[END_REF]) calculated with the R software (R Development Core Team 2011) was used to investigate the relatedness across sites.

RESULTS

Mitochondrial and chloroplast diversity

Between 18 and 26 tetrasporophytes were sequenced per site for both the COI and cox2-3 markers producing 301 mitochondrial sequences in total (Table 1). After editing, an alignment of 630 base pairs (bp) was built from 157 individuals for COI. Twenty-one mitotypes and 21 polymorphic sites were observed (GenBank accession numbers MF995393-MF995549). A total of 144 sequences of cox2-3 (358 bp) revealed five polymorphic sites and six mitotypes (GenBank accession numbers MG452409-MG452552). Once concatenated, the COI+cox2-3 sequences gave 26 mitotypes over the 983 bp studied on 144 individuals. After editing, an alignment of 1162 bp was built from 24 individuals for rbcL (i.e. 20 G. caudata and 4 G. secunda) and four chlorotypes and four polymorphic sites were observed in G. caudata.

Maximum likelihood (ML) tree reconstruction for the concatenated COI+cox2-3 sequence revealed the presence of two supported haplogroups distributed in distinct regions (Figure 2b). The first group corresponded to a north-eastern lineage and included all individuals sampled in CE, RN, PB and PE (i.e., ranging from 03 o 24'36"S to 08 o 07'59"S). The second group corresponded to a south-eastern lineage and included all individuals sampled in ES and SP (i.e., ranging from 20 o 48'31"S to 23 o 55'01"S). The sequences of individuals sampled in BA were retrieved as a polytomy from which south-eastern (SE) and north-eastern (NE) clades emerged (Figure 2b). In the ML tree reconstructed for rbcL, G. secunda from Mexico and USA appeared as a sister species to G. caudata (Figure 2a). Within G. caudata, the ML tree was very poorly resolved. However, G. caudata from ES, SP and Santa Catarina State (SC, Brazil) formed a monophyletic, albeit non-supported, branch separated from all the other samples collected from the north-eastern coast of Brazil, the USA and Venezuela. For the rbcL, a single base pair difference was observed between samples collected in the north-east and south-east regions (0.09% of divergence between MF995550 and MF995560).

The haplotype network for the concatenated COI+cox2-3 sequences supports the phylogenetic results (Figure 3) and shows two main haplogroups. Seventeen haplotypes (i.e., C1-C17 in Figure 3)

were found in the north-eastern part of the country in CE, RN, PB and PE. The two most widespread haplotypes, C1 and C5, were distributed in all four sites and the C13 haplotype was observed in the two adjacent sites of PB and PE. On the other hand, nine haplotypes (i.e., C18-C26 in Figure 3) were found in BA and in the south-eastern part of the country (ES and SP). No haplotypes were shared between the north-eastern and the BA plus south-eastern coasts or between sites sampled within the south-eastern haplogroup (Figure 3). Haplotype C25 was restricted to ES, C26 was restricted to SP, and seven haplotypes were private to BA. Concatenated mitotypes differed by 1-10 bp and only 2 bp separated the north-eastern and BA haplogroups (i.e., difference between C9 and C18 or C23, Figure 3).

A reduction in both gene (H) and nucleotide diversity (π) was observed from the north-eastern (NE) haplogroup to the BA and south-eastern monophyletic group (SE) for both mitochondrial markers (p <0.05 for COI and cox2-3, Table 1). The AMOVA indicated that total genetic variance was mainly explained by the variance among haplogroups (COI 58.72%, cox2-3 71.22%; Table 2b).

The variance within sites (COI 20.90%, cox2-3 16.18%) and the variance among sites within haplogroups (COI 20.38%, cox2-3 12.60%), although significant, were much lower (Table 2b). F ST values were significant between most pairs of sites, except for CE-RN (p >0.654) and PB-PE (p >0.384) for both markers, and for BA-ES (p >0.261) for the cox2-3 (Table 2a). The Mantel tests showed a clear correlation between genetic and geographic distances for COI and for cox2-3 (p <0.005) (Figure S1).

Microsatellite data: summary statistics, population differentiation and clustering analysis

A total of 411 diploid individuals (tetrasporophytes) of G. caudata were genotyped in this study and samples genotyped per site varied from 49 to 79 (Table 3). Fifteen polymorphic microsatellite markers were selected for the analyses, among the 17 loci previously developed to assess the genetic diversity of G. caudata (Ayres-Ostrock et al., 2016). The loci GraC_09 and GraC_11 were excluded due to very low amplification percentages (<57% and <47%, respectively).

The frequency of null alleles was significant for five loci (GraC_03, GraC_04, GraC_05, GraC_06 and GraC_12) and ranged from 0.070 to 0.349. Except for the RN site, which showed no evidence of null alleles for any locus, loci for which null alleles were detected and the frequency of null alleles varied with the site. The site with null alleles detected in the highest number of loci was SP (three loci) (Table 4). The number of alleles found in each site varied from one in CE, RN, PB, PE and SP for the locus GraC_10 to 16 for GraC_04 in the site of RN. Five loci were moderately polymorphic (4-6 alleles) and seven were highly variable (10-26 alleles) (Table 4). Evidence of linkage disequilibrium (LD) was found only in three of the seven sampled sites of G. caudata: CE, BA and SP. The locus GraC_03 showed strong and significant linkage disequilibrium with loci GraC_05 (SP, =0.492 and p <0.036) and GraC_07 (SP, =0.639 and p <0.036), as did the locus r d r d

GraC_08 with locus GraBC_04 (CE, =0.262 and p <0.011). All 15 microsatellite loci were used in r d the following analyses.

Regardless of the site under study, all individuals analysed had unique multi-locus genotypes and values of clonal richness (R) were equal to one. Observed heterozygosity (H o ) ranged from 0.288 (SP) to 0.479 (RN) and expected heterozygosity (H e ) ranged from 0.284 (SP) to 0.458 (RN) (Table 3). The average number of alleles per locus varied among sites, with the highest value encountered in RN (5.267) and the lowest value in SP (3.400) (Table 3). The geographic pattern of genetic variation observed with the microsatellite markers resembled that observed with the mitochondrial sequences.

Diversity was higher in the north-eastern area than in the south-eastern area (i.e., corresponding to north-eastern and BA+south-eastern mitochondrial haplogroups, respectively) and allelic richness (A R ) varied from 4.163 (NE) to 3.440 (BA+SE) and H e varied from 0.446 (NE) to 0.363 (BA+SE).

These differences were not significant (p >0.05) (Table 3).

Estimates of the inbreeding coefficient (F IS ) varied among loci from 0.000 in loci GC_08 (ES and SP), GBC_01 (SP), GBC_03 (CE, PE and BA) and GBC_04 (BA) to one in the locus GC_12

(SP) (Table 4). Multi-locus F IS values were generally close to zero and non-significant, except for ES (F IS =-0.125) and BA (F IS =0.115) (Table 3). The existence of a recent bottleneck in G. caudata was evaluated using Bottleneck software. Results indicate that all our sampled sites were in mutationdrift equilibrium. Unlike estimates previously obtained for mitochondrial markers, AMOVA analysis calculated using the microsatellites indicated that total genetic variance was mainly explained by the variance within sites 55.58% (p <0.0001, Table 5b). The variance among haplogroups (29.34%, p <0.05) and the variance among sites within haplogroups (15.07%, p <0.0001), although significant, were lower (Table 5b). Moreover, estimates of genetic differentiation (F ST ) among sites of G.

caudata ranged from 0.038 (PB-PE) up to 0.476 (SP-CE) (Table 5a). The Mantel test indicated that there was a clear correlation between genetic and geographic distances (Figure S2).

The Bayesian clustering method implemented in Structure software clearly showed that K =2

was the optimal number of clusters in our study (Figure 3c). For K =2, Structure results showed that all individuals from the north-eastern part of the coast (CE, RN, PB and PE) clustered together and the three remaining sites of BA, ES and SP formed the other group. Except in BA, a very low level of admixture was visible between the two genetic groups. When the Structure analysis was carried out within each genetic group, K =3 clusters were encountered in the north-eastern and the BA+south-eastern genetic group. Clustering within each geographic group generally separated the distinct sites. The level of admixture was much higher in the north-eastern (NE) than in the BA+SE (south-eastern) genetic group, where all individuals from the same site belonged to the same genetic cluster. In the north-eastern genetic group, a high genetic similarity was observed between the two neighbouring sites of PB and PE, with the individuals from PE showing traces of admixture between two genetic groups (Figure 3c). A similar pattern of genetic structure was also observed in the discriminant analysis of principal components (DAPC), where six clusters of genetically related individuals were identified, mostly related to their sampling sites, except for PB and PE (Figure S3).

Allele size variation depended on the sampling region for three loci (GraC_03, GraBC_02 and GraBC_03) (Figure 4). For GraC_03, private alleles were encountered in BA and in each of the two most south-eastern sites (ES and SP) and only a low number of these (N=4) were observed in the north-eastern part of G. caudata distribution at a very low frequency (Figure 4). For GraBC_02, one frequent allele was observed in the three most south-eastern sites (i.e., BA, ES and SP; most common allele size =258 bp), but another allele, present in high frequency, characterised the four most northeastern sites of CE, RN, PB and PE (most common allele size =264 bp) (Figure 4). For GraBC_02, alleles 258 bp and 264 bp were not shared between the two genetic groups (i.e., between the northeastern and BA+south-eastern genetic group). For GraBC_03, the five north-eastern sites of CE, RN, PB, PE and BA shared the same common allele, whereas another allele was encountered in samples collected in the States of Espírito Santo (ES) and SP. In general, the sites of ES and SP from the south-east presented a higher frequency of private alleles (0.177) than the north-east states (0.073)

and BA (0.068) (Table 3). Values of divergence observed between these haplogroups (rbcL 0.09% and COI 0.72%) are two times lower than those obtained between the two northern and southern phylogeographic groups of the intertidal Hypnea pseudomusciformis studied in the same geographic region (rbcL 0.2% and COI 1.3%) [START_REF] Nauer | Description of Hypnea pseudomusciformis sp. nov., a new species based on molecular and morphological analyses, in the context of the H. musciformis complex (Gigartinales, Rhodophyta)[END_REF][START_REF] Nauer | Phylogeography of the Hypnea musciformis species complex (Gigartinales, Rhodophyta) with the recognition of cryptic species in the western Atlantic Ocean[END_REF]. Nevertheless, these values are quite low when compared with pairs of recently diverged red macroalgae species (Iridaea cordata cryptic species found on both side of the Drake Passage, at the tip of South America in Tierra del Fuego and along the Antarctic Peninsula: rbcL=3.17% and COI=8.31%, Ocaranza-Barrera et al. 2018; Bostrychia intricata cryptic species N3 and N4 found in New Zealand Sub-Antarctic islands and New Zealand main islands:

DISCUSSION

rbcL=0.30% and COI=3.95%, [START_REF] Muangmai | Contrasting patterns of population structure and demographic history in cryptic species of Bostrychia intricata (Rhodomelaceae, Rhodophyta) from New Zealand[END_REF], driven and linked to historical events and changing environments.

Estimates of divergence time in G. caudata, based on the percentage of divergence per million years obtained in I. cordata and B. intricata, suggest that the effective separation between the north-eastern (NE) and BA+south-eastern (SE) haplogroups occurred recently, during the Quaternary, between 100,000 and 1,500,000 years ago. Apart from this major phylogenetic break (between the PE and BA regions), mtDNA and microsatellite markers showed an additional level of genetic structure, with each site located south of BA (i.e., ES, 20°S and SP, 23ºS) composed of a single unique mtDNA haplotype and constituting strongly divergent nuclear group.

A phylogeographic pattern dually influenced by marine and terrestrial climatic fluctuations in the past

The effect of the Pleistocene glacial cycles on the phylogeography of tropical species is much less studied than that of temperate ones [START_REF] Hewitt | Genetic consequences of climatic oscillations in the Quaternary[END_REF]). The effects of climatic change during the Pleistocene in tropical and sub-tropical regions in southern America restricted species' distribution, but allowed the persistence of populations (refugia), which remained isolated through multiple glacial-interglacial oscillations [START_REF] Lira | Long-term population isolation in the endangered tropical tree species Caesalpinia echinata Lam. revealed by chloroplast microsatellites[END_REF]. Therefore, strong signals of population expansion during interglacial periods are generally not detected [START_REF] Lessa | Genetic footprints of demographic expansion in North America, but not Amazonia, during the Late Quaternary[END_REF]). In the terrestrial realm, studies suggest that there were three Late Quaternary (with the most recent glacial period occurring between about 120,000 and 20,000 years before present) main refugia in the rainforest located along the Brazilian Atlantic coast (i.e. in the regions of PE, BA and SP), which fostered the long-term persistence of genetically isolated lineages (Carnaval andMoritz 2008, Carnaval et al. 2009).

Interestingly, these three main terrestrial refugia (see Figure 1) fit well with the two main phylogeographic breaks revealed in this study and correspond to areas where a high level of private alleles and haplotypes were detected in G. caudata.

The phylogeographic structure of marine species has been shaped by sea-level fluctuations during the Pleistocene [START_REF] Suguio | Flutuações do nível relativo do mar durante o Quaternário Superior ao longo do litoral brasileiro e suas implicações na sedimentação costeira[END_REF][START_REF] Martin | Fluctuating Holocene sea levels in eastern and southeastern Brazil: evidence from multiple fossil and geometric indicators[END_REF] as well as ocean current circulation patterns. Similar to our study, several phylogeographic studies have described one or two major genetic breaks along the Brazilian coast in invertebrate and vertebrate marine species (in the king weakfish, Macrodon ancylodon, [START_REF] Santos | Population genetic structuring of the king weakfish, Macrodon ancylodon (Sciaenidae), in Atlantic coastal waters of South America: deep genetic divergence without morphological change[END_REF] [START_REF] Francisco | Population genetic structure, introgression, and hybridization in the genus Rhizophora along the Brazilian coast[END_REF]).

However, the exact limits of the breaks do not always coincide among the studied species. In most studies, the first genetic barrier is located at 4°S-6°S, separating the populations of northern Brazil from those of central Brazil [START_REF] Santos | Population genetic structuring of the king weakfish, Macrodon ancylodon (Sciaenidae), in Atlantic coastal waters of South America: deep genetic divergence without morphological change[END_REF][START_REF] Francisco | Population genetic structure, introgression, and hybridization in the genus Rhizophora along the Brazilian coast[END_REF][START_REF] Peluso | Contemporary and historical oceanographic processes explain genetic connectivity in a Southwestern Atlantic coral[END_REF]) and the second one occurs at 21°S-23°S, separating the tropical populations located along the coast where the Brazil Current flows from the sub-tropical populations influenced by cooler South Atlantic Central Waters (SACW) [START_REF] Santos | Population genetic structuring of the king weakfish, Macrodon ancylodon (Sciaenidae), in Atlantic coastal waters of South America: deep genetic divergence without morphological change[END_REF][START_REF] Da Silva Cortinhas | Genetic structuring among silverside fish (Atherinella brasiliensis) populations from different Brazilian regions[END_REF][START_REF] Dos Santos Freitas | The mitochondrial control region reveals genetic structure in southern kingcroaker populations on the coast of the Southwestern Atlantic[END_REF][START_REF] Peluso | Contemporary and historical oceanographic processes explain genetic connectivity in a Southwestern Atlantic coral[END_REF]. All these studies suggest ocean currents (Figure 1) play a role as effective barriers to gene flow.

In intertidal organisms, such as G. caudata, the phylogeographic structure is probably best explained by the interaction between terrestrial history and marine current patterns. Other features such as the outflows from the Doce and São Francisco rivers may also act as biogeographical barriers separating populations located along the Brazilian coast [START_REF] Schmid | The Vitória eddy and its relation to the Brazil Current[END_REF][START_REF] Carnaval | Stability predicts genetic diversity in the Brazilian Atlantic forest hotspot[END_REF], Figure 1). Indeed, even if G. caudata has been described as euryhaline species [START_REF] Yokoya | Effects of salinity on the growth rate, morphology and water content of some Brazilian red alga of economic importance[END_REF], experiments have shown that it was not able to grow at salinities lower than 15 PSU (de [START_REF] De Miranda | Effects of temperature, salinity and irradiance on carposporeling development of Hidropuntia caudata (Gracilariales, Rhodophyta)[END_REF]) and it can be assumed that strong freshwater outflows may have a significant impact on the fitness of these intertidal populations. The possible role of the São Francisco River as a promoter of differentiation has been commonly reported in terrestrial organisms [START_REF] Pellegrino | Phylogeography and species limits in the Gymnodactylus darwinii complex (Gekkonidae, Squamata): genetic structure coincides with river systems in the Brazilian Atlantic Forest[END_REF][START_REF] Carnaval | Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest[END_REF][START_REF] Carnaval | Stability predicts genetic diversity in the Brazilian Atlantic forest hotspot[END_REF]). However, there is congruence between river systems and the borders of putative long-term refugia in the area [START_REF] Carnaval | Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest[END_REF], and some authors have proposed that primary divergence was due to isolation in climatic refugia, rivers being only secondary barriers maintaining the resulting genetic structure [START_REF] Carnaval | Prediction of phylogeographic endemism in an environmentally complex biome[END_REF]; but see [START_REF] Thomé | Barriers, rather than refugia, underlie the origin of diversity in toads endemic to the Brazilian Atlantic Forest[END_REF]. Additionally, along the Brazilian coast, the Holocene sea-level history suggests regional variations, with alternating periods of mean sea-level rise and fall resulting in huge reductions in population size and decreased connectivity. The frequency, time and impact of these oscillations during the past 7000 years have been well studied [START_REF] Angulo | The Brazilian sea-level curves: a critical review with emphasis on the curves from the Paranaguá and Cananéia regions[END_REF][START_REF] Martin | GC "The Brazilian sea-level curves: a critical review with emphasis on the curves from the Paranaguá and Cananéia regions[END_REF], Lessa and Angulo 1998[START_REF] Martin | Fluctuating Holocene sea levels in eastern and southeastern Brazil: evidence from multiple fossil and geometric indicators[END_REF]) and these oscillations have had a strong influence on the distribution and diversity of several coastal species (e.g., in the sand dune ant Mycetophylax, [START_REF] Cardoso | Phylogeography of the sand dune ant Mycetophylax simplex along the Brazilian Atlantic Forest coast: remarkably low mtDNA diversity and shallow population structure[END_REF]; in five small mammal species widespread in the Brazilian Atlantic forest, [START_REF] Leite | Neotropical forest expansion during the last glacial period challenges refuge hypothesis[END_REF]; and in the red mangrove Rhizophora mangle, [START_REF] Pil | Postglacial north-south expansion of populations of Rhizophora mangle (Rhizophoraceae) along the Brazilian coast revealed by microsatellite analysis[END_REF]. Depending on the study, sea-level fluctuations had amplitudes of 2 to 4 m when compared to today's sea-level, and durations of 400 to 500 years [START_REF] Suguio | Flutuações do nível relativo do mar durante o Quaternário Superior ao longo do litoral brasileiro e suas implicações na sedimentação costeira[END_REF][START_REF] Martin | Fluctuating Holocene sea levels in eastern and southeastern Brazil: evidence from multiple fossil and geometric indicators[END_REF]. These oscillations were characterized by rates in sea level change of 16 up to 32 mm per year, three times faster than the fastest sea-level rise period (10 mm/year, between 15000 and 8000 yr BP - [START_REF] Angulo | A critical review of mid-to late-Holocene sealevel fluctuations on the eastern Brazilian coastline[END_REF], that could have strongly affected marine species distribution along the RN-SC coast and in adjacent regions (Fleming et al. 1998). Strikingly, receding sea levels exposed as much as 92% of today's Brazilian continental shelf [START_REF] Araújo | Quaternary climate changes explain diversity among reptiles and amphibians[END_REF][START_REF] Carnaval | Stability predicts genetic diversity in the Brazilian Atlantic forest hotspot[END_REF], with significantly more land being exposed south of BA (i.e., Abrolhos Bank,18 o S) [START_REF] Leite | Neotropical forest expansion during the last glacial period challenges refuge hypothesis[END_REF]. There are many reports worldwide of genetic evidence of bottlenecks in coastal marine taxa in tropical regions due to the global reduction in sea level during the Pleistocene [START_REF] Ludt | Shifting seas: the impacts of Pleistocene sea-level fluctuations on the evolution of tropical marine taxa[END_REF]. Such climatic events may explain the reduced mitochondrial and nuclear genetic diversity and high genetic structure observed between sites of G.

caudata located south of BA (ES and SP) compared with the north-eastern sites. Moreover, the fluctuation of the sea levels along the Vitoria-Trindade Ridge was reported as a possible cause of the high level of endemism in this coral reef environment because of successive bottlenecks [START_REF] Pinheiro | Island biogeography of marine organisms[END_REF]. Our hypothesis is that such successive bottlenecks will impact more strongly the genetic structure of intertidal benthic organisms rather than subtidal pelagic species.

Secondary contact zone, but low level of admixture in the Bahia (BA) region

The Bahia (BA) region presents a high number of divergent haplotypes, suggesting that one large or several refugia may have withstood the successive glacial cycles of the Pleistocene in this region. Mitochondrial markers point out the existence of large and stable populations in BA, historically isolated from both the south-eastern and the north-eastern parts of the G. caudata distribution in Brazil. The existence of major long-standing BA refugia in the coastal Brazil Atlantic forest has also been demonstrated for terrestrial plants and animals and is associated with higher species and genetic endemism, compared with areas located south of the Doce River (i.e., our sites of ES and SP) (Carnaval andMoritz 2008, Carnaval et al. 2009). In G. caudata, the mtDNA tree shows that individuals sampled in this region constitute a basal polytomy to both the north-eastern and south-eastern monophyletic branches. This result is supported by the private BA haplotypes, which occupy an intermediate position in the mtDNA haplotype network, connecting the northern-most and southern-most haplotypes. However, depending on the locus, diagnostic microsatellite alleles of the north-eastern and south-eastern regions were both observed in the BA population (Figure 4) and BA was the only case in which a low level of admixture between the north-eastern and south-eastern lineages was detected. All these results suggest that BA may correspond not only to a refugium, but also to a recent secondary contact zone between the highly differentiated north-eastern and southeastern regions. The occurrence of a contact zone in the BA region can also be surmised from a previous phylogenetic study conducted in another red alga species (H. pseudomusciformis, [START_REF] Nauer | Phylogeography of the Hypnea musciformis species complex (Gigartinales, Rhodophyta) with the recognition of cryptic species in the western Atlantic Ocean[END_REF]. As in G. caudata, a clear divergence between the samples from north-east and the samples from the south-east of Brazil was detected using COI mitochondrial sequences and, in both studies, the BA region is an intermediate region that is either grouped with the north-eastern (in H. pseudomusciformis) or with the south-eastern cluster (in G. caudata).

Both historical isolation and more recent micro-evolutionary processes have probably played a major role in shaping G. caudata genetic structure in the BA region. Although we did not observe the gradient of admixture and/or clines of allele frequencies expected in this contact zone, more extensive sampling in this zone is required to clarify the size of the contact zone and level of introgression between the two main lineages. It is possible that the genetic structure is currently maintained by the combination of the low dispersal capacity of the species, the existence of various barriers to gene flow and adaptation to sharp changes in environmental conditions (e.g., salinity, sea surface temperature, oxygen and nutrient concentration; [START_REF] Briani | The influence of environmental features in the content of mycosporine-like amino acids in red marine algae along the Brazilian coast[END_REF]). In seaweeds with low dispersal capacity, even weak and/or transient barriers can easily promote and maintain strong genetic divergence [START_REF] Montecinos | Species replacement along a linear coastal habitat: phylogeography and speciation in the red alga Mazzaella laminarioides along the south east pacific[END_REF]. Studies on different intertidal red seaweed species revealed restricted gene flow even at short distance. Particularly between tide pools located at different levels of the shore and exposed to contrasting levels of physiological stress [START_REF] Engel | Performance of non-motile male gametes in the sea: analysis of paternity and fertilization success in a natural population of a red seaweed, Gracilaria gracilis[END_REF], 2004[START_REF] Krueger-Hadfield | Genetic population structure and mating system in Chondrus crispus (Rhodophyta)[END_REF], 2013, 2015[START_REF] Maggs | Speciation in red algae: members of the Ceramiales as model organisms[END_REF]. Such patterns of reduced gene flow, if congruent with contrasted environmental conditions, may also initiate ecotypic differentiation. [START_REF] Briani | The influence of environmental features in the content of mycosporine-like amino acids in red marine algae along the Brazilian coast[END_REF] reported that the content of mycosporine-like amino acids (MAAs) in 39 different red algae including G. caudata varies sharply between the north-east and south-east of the Brazilian coast. Algae MAAs content varied mainly according to the level of UV and irradiance, however, other environmental factors such as pH and phosphate and nitrate concentrations in the water column also had some influence [START_REF] Briani | The influence of environmental features in the content of mycosporine-like amino acids in red marine algae along the Brazilian coast[END_REF]. The genetic break revealed in G. caudata thus fits nicely with the sharp abiotic environmental break detected between tropical and subtropical regions [START_REF] Briani | The influence of environmental features in the content of mycosporine-like amino acids in red marine algae along the Brazilian coast[END_REF].

Ecotypic differentiation and conservation management issues

Differences in physiological performance (e.g., growth rate, photosynthesis parameters, pigment content and sensitivity to UVB radiation) were observed between strains of G. caudata from CE (corresponding to the NE clade) and SP (corresponding to the SE group) grown in controlled laboratory conditions [START_REF] Araújo | Intraspecific variation in Gracilaria caudata (Gracilariales, Rhodophyta): growth, pigment content, and photosynthesis[END_REF][START_REF] Faria | Ecotypes of Gracilaria caudata (Gracilariales, Rhodophyta): physiological and morphological approaches considering life history phases[END_REF]. Possible ecotypic differentiation may have arisen along the G. caudata distribution range, linked to local adaptation to habitats specific to the north-east and the regions of lower latitude (such as MAAs content, see [START_REF] Briani | The influence of environmental features in the content of mycosporine-like amino acids in red marine algae along the Brazilian coast[END_REF]). These results were also confirmed by the Mantel test that showed a significant correlation between genetic and geographic distances, suggesting that the population structure and the subsequent occurrence of ecotypic differentiation could be explained by an isolation-by-distance model.

Interestingly, differences in seawater temperature linked to ocean circulation patterns have also been detected along the Brazilian coast (average ranging between 20 and 28°C in the northeast region and between 6 and 24°C in the southeast region); and are regarded as a potential physiological barrier to gene flow for some species of fish [START_REF] Galetti | Assessing genetic diversity of Brazilian reef fishes by chromosomal and DNA markers[END_REF][START_REF] Santos | Population genetic structuring of the king weakfish, Macrodon ancylodon (Sciaenidae), in Atlantic coastal waters of South America: deep genetic divergence without morphological change[END_REF][START_REF] Machado | Population genetic structure and demographic history of the spadefish, Chaetodipterus faber (Ephippidae) from Southwestern Atlantic[END_REF]). Temperature ecotypes have been described along the Brazilian coast for another Gracilaria species: G. birdiae [START_REF] Ursi | Intraspecific variation of photosynthesis, respiration and photoprotective carotenoids in Gracilaria birdiae (Gracilariales: Rhodophyta)[END_REF]. Differences in sensitivity to UVB radiation were also detected between G. birdiae strains from CE and SP states (Ayres and Plastino 2014). Physiological divergence between habitats, linked to distinct resilience capacity to temperature fluctuation and emersion (i.e., desiccation) stress, has been shown in other intertidal seaweeds, for example in the brown alga Fucus vesiculosus [START_REF] Nicastro | Shift happens: trailing edge contraction associated with recent warming trends threatens a distinct genetic lineage in the marine macroalga Fucus vesiculosus[END_REF][START_REF] Saada | Taking the heat: distinct vulnerability to thermal stress of central and threatened peripheral lineages of a marine macroalga[END_REF]. Ecotypic divergence has been shown to lead to rapid speciation in several organisms including macroalgae (i.e., adaptation to low salinity in the Baltic sea in Fucus in the time span of a few hundreds of years, [START_REF] Pereyra | Rapid speciation in a newly opened postglacial marine environment, the Baltic Sea[END_REF].

Nevertheless, a recent study, based on experimental crosses, demonstrated that individuals of G.

caudata from north-eastern (i.e, CE) and south-eastern (i.e, ES and SP) clades produced viable and fertile descendants, suggesting that complete reproductive barriers may not yet have evolved between the two clades [START_REF] Chiaramonte | Gracilaria caudata (Gracilariales, Rhodophyta) is reproductively compatible along the whole Brazilian coast[END_REF].

In conclusion, our study demonstrates the existence of divergent clades -which based on previous physiological work may correspond to distinct ecotypes [START_REF] Araújo | Intraspecific variation in Gracilaria caudata (Gracilariales, Rhodophyta): growth, pigment content, and photosynthesis[END_REF][START_REF] Faria | Ecotypes of Gracilaria caudata (Gracilariales, Rhodophyta): physiological and morphological approaches considering life history phases[END_REF]) -in G. caudata and thus should be considered for the management policy of this commercially important species. Harvesting of natural Gracilaria beds in Brazil has gradually diminished since the 1970s and has been replaced by small-scale cultivation [START_REF] Hayashi | Cultivation of red seaweeds: a Latin American perspective[END_REF].

Nevertheless, since the beginning of this activity in the 1960s, over-harvesting of the natural beds in the north-eastern region has compromised one of the main sources of genetic diversity.

Acknowledging the existence of multiple phylogeographical lineages along the Brazilian coast is important not only for understanding the recent historical processes shaping genetic diversity in these tropical and sub-tropical regions, but also for developing effective conservation strategies, We are also most grateful to the Biogenouest Genomics core facility for its technical support.

We thank the anonymous referees for valuable suggestions.

BIOSKETCH boundaries, redrawn from the most historically stable Atlantic forest areas proposed in [START_REF] Carnaval | Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest[END_REF]. Please note that the number and localization of refugia predicted for the Brazilian Atlantic forest south of the Doce River was not as well defined as for the northern part of the coast [START_REF] Carnaval | Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest[END_REF] and that the existence of three main refugia have been proposed: 1) Pernambuco (i.e., PE), 2) Bahia (i.e., BA) and 3) a small refuge with a localization difficult to predict that could be at the boundaries of the states of Espírito Santo and Rio de Janeiro (i.e., ES+SP); the last refugium corresponding to a unique center of endemism identified for mammals, birds, butterflies and plants.

  Throughout the Brazilian coast, two major genetic clusters, distributed in strict parapatry, were revealed within the red alga G. caudata based on both phylogenetic inferences and population genetic studies. There was a remarkable congruence between nuclear, mitochondrial and chloroplast data showing a clear separation between the north-east (i.e., ranging from 03 o S to 08 o S) and the south-east of the country (i.e., ranging from 20°S to 23 o S), supporting the existence of a long-term divergence along the Brazilian coast. The BA region (around 12 o S) occupies an intermediate position between both clades. These lineages also differ by their demographic histories: signatures of recent demographic expansions in the north-east cluster, whereas the BA and south-eastern haplogroups are composed of highly divergent populations, suggesting the maintenance of several different refugia.

  particularly in an environment subject to important anthropogenic factors like overharvesting, habitat fragmentation and degradation.ACKNOWLEDGMENTSThis research was supported by the São Paulo Research Foundation (FAPESP: 2010/50175-3; 2011/10189-8), the Brazilian National Council for Scientific and Technological Development (CNPq: 300148/93-3; 301491/2013-5), the programme of international cooperation USP/COFECUB between the University of São Paulo and the Comité Français d'Evaluation de la Coopération Universitaire avec le Brésil and the International Research Network DEBMA 'Diversity, Evolution and Biotechnology of Marine Algae' (GDRI CNRS 0803). This study was funded in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) -Finance Code 001.

FIGURE 1

 1 FIGURE 1 Map of the Brazilian coast showing the direction of major ocean currents: North Brazil

FIGURE 2

 2 FIGURE 2 Maximum likelihood (ML) trees based on: (a) rbcL data and (b) concatenated mtDNA

  FIGURE S1Mantel test of the relationship between genetic and geographic distances, using the mitochondrial markers (a) COI and (b) cox2-3, for the seven sampled Gracilaria caudata localities along the Brazilian coast, expressed as FST versus the distance in kilometres (km). The red dots represent relationship between genetic and geographic distances of among the north-east populations. The blue dots represent the relationship between genetic and geographic distances of the Bahia location and the remaining populations. The black dots represent relationship between genetic and geographic distances of among the south-east and the remaining populations.

  ; in the silver fish Atherinella brasiliensis, da Silva Cortinhas et al. 2016; in the king fish Menticirrhus americanus, dos Santos Freitas et al. 2017;in the white shrimp Litopenaeus schmitt[START_REF] Maggioni | Population structure of Litopenaeus schmitti (Decapoda: Penaeidae) from the Brazilian coast identified using six polymorphic microsatellite loci[END_REF]; in the coral Mussismilia hispida[START_REF] Peluso | Contemporary and historical oceanographic processes explain genetic connectivity in a Southwestern Atlantic coral[END_REF]) and in the red mangrove tree Rhizophora mangle

TABLE 1

 1 Sampling sites for Gracilaria caudata and their genetic diversity for two mitochondrial markers (COI and cox2-3) in Brazil.

	Molecular diversity estimates calculated for COI (630 bp) and cox2-3 (358 bp): number of sequences (N), number of haplotypes (nH), gene
	diversity (H), nucleotide diversity (π), number of polymorphic sites (S), Tajima's D (D) and Fu's F S (F S ) statistics. Standard deviation (SD) is
	given in parentheses. Ceará (CE), Rio Grande do Norte (RN), Paraíba (PB), Pernambuco (PE), Bahia (BA), Espírito Santo (ES) and São Paulo
	(SP). North-eastern haplogroup: CE, RN, PB and PE; BA and south-eastern haplogroup: BA, ES and SP. All specimens were collected in-
	between high and low intertidal areas uncovered at low tide
	Site (State)	Coordinates (GPS)					COI								cox2-3
			N	nH	H (SD)	𝜋	(*10 -2 ) (SD)	S	D	F S	N	nH	H (SD)	𝜋	(*10 -2 ) (SD)	S	D	F S
	CE	03º24'36.0"S/39º01'50.0"W	26	10	0.750 (0.074)	0.270 (0.182)	12 -1.542* -4.389** 24	2	0.463 (0.069) 0.129 (0.131)	1	1,231	1,362
	RN	05º15'41.0"S/35º23'11.0"W	22	5	0.649 (0.064)	0.220 (0.157)	6 -0.482 -0.145	19	3	0.374 (0.129) 0.153 (0.147)	2	-0.094	-0.071
	PB	07º 17'62.9"S/34 º 48'08.5"W	22	5	0.528 (0.118)	0.186 (0.138)	5 -0.430 -0.582	22	2	0.311 (0.106) 0.087 (0.103)	1	0.236	0.647
	PE	08º07'58.9"S/34º53'57.3"W	23	4	0.557 (0.083)	0.173 (0.131)	4 0.018	0.347	22	3	0.536 (0.090) 0.232 (0.192)	3	0.025	0.890
	BA	12º44'28.0"S/38º09'01.0"W	23	4	0.525 (0.094)	0.097 (0.089)	3 -0.615 -0.964	20	4	0.363 (0.130) 0.230 (0.001)	4	-0.773	-0.443
	ES	20º48'31.0"S/40º36'39.0"W	23	2	0.087 (0.077)	0.027 (0.042)	2 -1.514* -0.153	18	1	0.000 (0.000) 0.000 (0.000)	0	0.000	0.000
	SP	23º55'01.0"S/46º19'16.8"W	18	1	0.000 (0.000)	0.000 (0.000)	0 0.000	0.000	19	1	0.000 (0.000) 0.000 (0.000)	0	0.000	0.000
		Total G. caudata 157	21	0.442 (0.073)	0.139 (0.106)	21 -0.652 -0.841	144	6	0.292 (0.075) 0.119 (0.109)	5	0.089	0.340
		North-eastern haplogroup 93	16	0.621 (0.085)	0.212 (0.152)	16 -0.609 -1,192	87	4	0.421 (0.099) 0.150 (0.143)		0.350	0.707
		BA + south-eastern haplogroup 41	5	0.204 (0.057)	0.041 (0.043)	5 -0.710 -0.558	37	3	0.121 (0.043) 0.076 (0.064)		-0.257	-0.443
	*p< 0.05; **p< 0.001
	*P < 0.05; **P < 0.001

TABLE 2 (

 2 a) Pairwise estimates of genetic differentiation (F ST ) between the seven sampled Gracilaria caudata sites for COI (below the diagonal) and cox2-3 (above the diagonal) and (b)

	analysis of molecular variance (AMOVA) within sites, among sites within haplogroup and among
	haplogroups for both molecular markers. Haplogroups were defined according to the haplotype
	network (see Figure 3; north-eastern haplogroup: CE, RN, PB and PE; Bahia and south-eastern
	haplogroup: BA, ES and SP). Ceará (CE), Rio Grande do Norte (RN), Paraíba (PB), Pernambuco
	(PE), Bahia (BA), Espírito Santo (ES) and São Paulo (SP). Degrees of freedom (d.f.) and sum of
	squares (SS).
	(a)
	CE	RN	PB	PE	BA	ES	SP
	CE (a)	0.004 0.357** 0.111* 0.701** 0.887** 0.922**
	RN -0.020		0.460** 0.199** 0.674** 0.878** 0.917**
	PB 0.209** 0.155*		0.039 0.776** 0.939** 0.956**
	PE 0.259** 0.211** 0.003		0.628** 0.815** 0.871**
	BA 0.636** 0.690** 0.723** 0.727**		0.073 0.709**
	ES 0.804** 0.850** 0.870** 0.874** 0.834**		1.000**
	SP 0.805** 0.852** 0.875** 0.879** 0.856** 0.948**
	*p< 0.05; **p< 0.001
	*P < 0.05; **P < 0.001

TABLE 3

 3 Genetic diversity of Gracilaria caudata calculated for 15 nuclear microsatellite loci. For each site, the number of individuals genotyped is indicated (n). Multi-locus mean estimates of the number of alleles per locus (N a ), expected heterozygosity (H e ), observed heterozygosity (H o ), allele richness (A R ), frequency of private alleles (F PRIV ) and deviation from random mating (F IS ). F IS values significantly different from zero are shown in bold. Standard deviation (SD) is given in parentheses. Ceará (CE), Rio Grande do Norte (RN), Paraíba (PB), Pernambuco (PE), Bahia (BA), Espírito Santo (ES) and São Paulo (SP). Haplogroups were defined using the mitochondrial data set. North-eastern haplogroup: CE, RN, PB and PE; BA and south-eastern haplogroup: BA, ES and SP.

	Site (State)	Coordinates (GPS)			Microsatellite Loci
			n	N a	H e (SD)	H o (SD)	F IS	A R	F PRIV
	CE	03º24'36.0"S/ 39º01'50.0"W	72	4.267	0.379 (0.248)	0.384 (0.282)	-0.006	3.803	0.047
	RN	05º15'41.0"S/ 35º23'11.0"W	52	5.267	0.458 (0.249)	0.479 (0.280)	-0.037	4.992	0.076
	PB	07º17'62.9"S/ 34º48'08.5"W	79	4.533	0.386 (0.242)	0.402 (0.298)	-0.034	4.077	0.103
	PE	08º07'58.9"S/ 34º53'57.3"W	53	4.000	0.366 (0.272)	0.379 (0.296)	-0.027	3.780	0.067
	BA	12º44'28.0"S/ 38º09'01.0"W	55	3.933	0.333 (0.275)	0.298 (0.250)	0.115	3.692	0.068
	ES	20º48'31.0"S/ 40º36'39.0"W	51	3.733	0.345 (0.301)	0.392 (0.376)	-0.125	3.582	0.196
	SP	23º55'01.0"S/ 46º19'16.8"W	49	3.400	0.284 (0.279)	0.288 (0.309)	-0.001	3.299	0.157

TABLE 4

 4 Genetic variability within the sampled sites of Gracilaria caudata. For each site and each locus, the number of individuals genotyped (n) is indicated. The number of alleles per locus (N a ), expected heterozygosity (H e ), observed heterozygosity (H o ), frequency of private alleles (Freq. N priv ) and estimates of deviation from random mating (F IS ) was calculated for the 15 selected microsatellite loci independently. F IS values significantly different from zero are shown in bold. Presence of null allele frequency was tested with MICRO-CHECKER software[START_REF] Van Oosterhout | MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data[END_REF]) and loci for which frequency of null alleles was significant are indicated. Mean and standard error (SE) computed over the 15 loci are given in parentheses.

TABLE 5 (

 5 a) Pairwise estimates of genetic differentiation (F ST ) between the seven sampled Gracilaria caudata sites for the 15 selected microsatellite loci and (b) analysis of molecular variance (AMOVA) within sites, among sites within haplogroup and among haplogroups for the nuclear microsatellite markers. Haplogroups were defined according to the haplotype network (see Figure3, north-eastern haplogroup: CE, RN, PB and PE; BA and south-eastern haplogroup: BA, ES and SP).

	Ceará (CE), Rio Grande do Norte (RN), Paraíba (PB), Pernambuco (PE), Bahia (BA), Espírito Santo
	(ES) and São Paulo (SP). Degrees of freedom (d.f.) and sum of squares (SS).
		CE	RN	PB	PE	BA	ES	SP
	CE
	RN	0.173*
	PB	0.140* 0.217*
	PE	0.085* 0.204* 0.038*
	BA	0.309* 0.393* 0.303* 0.297*
	ES	0.449* 0.456* 0.420* 0.447* 0.299*
	SP	0.477* 0.477* 0.417* 0.463* 0.381* 0.228*
								Microsatellites
	Source of variation			d.f.	SS	Variance	% variation	p-value
								components
	Among haplogroups			1	483.370	1.081	29.34	p< 0.05
	Among localities within haplogroups		5	335.866	0.555	15.07	p< 0.001
	Within localities				815	1670.016	2.049	55.58	p< 0.001
	Total					821	2489.252	3.686	-	-

Journal of Phycology

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest