Introduction

Global climate change has altered the distribution of many species and also affects the environmental conditions of many ecosystems (^{Hughes, 2000}). Some species have the ability to tolerate and adapt to changes in the environment (^{Aitken et al., 2008}). For example, in fish populations of the genus Fundulus, genes associated with acclimatization processes have more patterns of adaptive divergence than those genes that are not associated with favoring species adaptation (^{Whitehead et al., 2012}). Variations in gene expression within fish populations may be of great importance for evolutionary ecology, survival, and species distribution ranges (^{Rees et al., 2011}; ^{Townsend et al., 2003}). Gene expression is associated with adaptation that determines how organisms interact with the environment and responds under certain stress conditions (^{Bonduriansky & Day, 2009}).

An increase in water temperature affects physiological and biochemical processes in many organisms, particularly in fish, and this may limit their distribution (^{Sorensen et al., 2001}). The genus Fundulus has been useful for evaluating the response to conditions of thermal stress, because its phenotypic plasticity is high, and it is distributed in diverse ecosystems (^{Fangue et al., 2006}; ^{Picard & Schulte, 2004}; ^{Whitehead, 2010}). Species of Fundulus genus have wide physiologic tolerance to several environmental conditions, which allow them to have broad distribution in costal environments (^{Whitehead, 2010}). Populations of this genus have developed in coastal areas with high salinity and a wide temperature range (^{Schulte et al., 2000}), from -1.4 °C to 31 °C (^{Bulger & Tremaine, 1985}). Variation in the expression of temperature-associated genes (in particular heat shock proteins, Hsp) has been demonstrated through a latitudinal gradient in Fundulus heteroclitus (^{Fangue et al., 2006}). The Hsp avoid the denaturation of proteins and favor the conformation of the proteins to their natural form; they also participate in the elimination of damaged proteins (^{Bukau & Horwich, 1998}). The Hsp are important in the ecology of the fish genus Fundulus, allowing its adaptation to high temperatures (^{Parsell & Lindquist, 1993}), and they are good indicators of response to thermal stress. In F. heteroclitus, expression of Hsp genes increases when fish are subjected to thermal stress (^{Healy & Schulte, 2012}; ^{Healy et al., 2010}, ^{White et al., 1994}). This may help the fish to withstand variations in climate, mainly an increase in sea surface temperature such as that caused by global climate change (^{Kayhan & Duman, 2010}).

Two species of the genus Fundulus are endemic to the Yucatán Peninsula: F. grandissimus and F. persimilis inhabit coastal lagoons along the northern coast of the peninsula (^{Vega-Cendejas & Hernández, 2004}). These fish are a source of food for local and migratory birds and are part of the diet of several other species of fish (^{De Dios, 2014}). Organisms of F. grandissimus have a limited distribution, while organisms of F. persimilis have a wide distribution. Either inhabiting brackish environments or salt water (^{Miller, 2009}) with high temperatures, at depths of < 1 m (^{Vega-Cendejas & Hernández, 2004}). Individuals of F. grandissimus and F. persimilis in coastal lagoons in Yucatán Peninsula rarely co-occur, and their distribution and presence in the lagoon is affected by environmental conditions (^{Álvarez, 2016}).

In this study, we evaluated gene expression of 2 variant Hsp in wild F. grandissimus and F. persimilis from coastal lagoons on the northern Yucatán Peninsula. The aim was to determine the response capacity of these species to temperature variations in their natural habitat by assessing variations in Hsp gene sequences, and expression, also to determine how thermal tolerance influence the distribution and dispersion of each species.

Materials and methods

Adults of F. grandissimus (37 individuals) were collected in Chelem Lagoon and F. persimilis (38 individuals) were collected in 2 lagoons (Chelem and Carbonera) on the northern coast of the Yucatán Peninsula (Fig. 1). In Chelem Lagoon, the size of the lagoon and the distribution of the species allowed fish to be collected from 3 different sites: in the sampling points ChM1: 21°15’47.05” N, 89°44’30.69” W and ChL1: 21º15’31.72” N, 89°42’7.35” W we found F. perismilis organisms, in the sampling point ChL2: 21°14’56.682” N, 89°43’40.73” W, was found F. grandissimus organism. In Carbonera Lagoon, only 1 organism of F. persimilis was found in the site Ca1: 21°13’54.35” N, 89°53’41.70” W (Fig. 1). At each sampling site, temperature, salinity, depth and dissolved oxygen were measured at depth of < 1m using a multiparameter Yellow Springs Instruments (YSI) probe, model 556 MPS. All measurements were recorded during the capture of the organisms on August 26^th, 2016 between 11:00 am and 2:00 pm. We collected in this time period since the highest density of fish was observed at that time, later or earlier the organisms were not found. Recorded water surface temperature was 33.38 °C at ChL2, the only site where F. grandissimus was collected. It was 31.94 ºC at ChM1; 31.39 ºC at ChL1, and 31.34 °C at Ca1, sites at which F. persimilis was collected.

Figure 1 Study area and points where Fundulus grandissimus and F. persimilis were collected; Chelem and Carbonera lagoons, northern Yucatán Peninsula, Mexico. ChM (Chelem dock); ChL (Chelem Lagoon); CaL (Carbonera Lagoon). 

A small fragment of muscle tissue from each fish was preserved in nucleic acids preserving buffer (for 1 L of buffer, 40 ml of EDTA 0.5M, 25 ml of sodium citrate 1M, 700 g of ammonium sulfate) and kept cold in ice until processed. DNA was extracted from ~10 mg of tissue with the Qiagen tissue DNA extraction kit. The quality of genetic material obtained was confirmed by means of a 1% electrophoresis gel.

The region of Hsp 70 gene was amplified from genomic DNA by the polymerase chain reaction (PCR), using the following primers: Hsp 70-1F 5 'CTC AGA TCT TTT CCA CGT ACT CA' 3; Hsp 70-1R: 5 'CTC CAG TAG TGA AAT GAT GCA GT' 3 (^{Fangue et al., 2006}). The PCR products were verified on 1.5% agarose gels stained with SYBR^® Green. Sequences of PCR product were obtained with a 16-capillary Genetic Analyzer sequencer (model 3130; Applied Biosystems, Foster City, CA) and BigDye Terminator V3.1, Cycle Sequencing Kit (Thermo Fisher Scientific, Waltham, MA). The verification of the electropherograms was analyzed with the free software 4Peaks v1.8 (Nucleobytes, Amsterdam, Netherlands). Sequences were aligned with the software Clustal IX Ver 2.1 (^{Thompson et al., 1994}), using the parameters determined by the program, and visually checked in the software Geneious v9 (^{Kearse et al., 2012}). Estimates of genetic distances between sequences were conducted using the Hasegawa-Kishino-Yano plus Gamma model (^{Hasegawa et al., 1985}) according to according to Akaike Information Criteria (AIC); analyses were performed in MODELTEST software (^{Posada & Crandall, 1998}). The analysis involved 75 DNA sequences of the Hsp 70 gene. There was a total of 578 positions in the final dataset, genetic distances were conducted in software Geneious ver. 9 (^{Kearse et al., 2012}).

Genetic comparisons were performed using as a reference the genetically more closely related species F. heteroclitus reported by ^{Fangue et al. (2006)} (accession number: DQ202279). The sequences were compared with the NCBI (National Center for Biotechnology Information) database by means of a BLAST (Basic Local Alignment Search Tool) analysis using software Geneious ver. 9 (^{Kearse et al., 2012}).

Nucleotide sequences of Hsp70 gene were translate to amino acid sequences (using standard genetic code, translocation frame 3), to predicts protein secondary structure with the original Garnier Osguthorpe Robson algorithm (GOR I) provided by the EMBOSS suite; both analyses performed in software Geneious ver. 9 (^{Kearsen et al., 2012}). Mc Donald and Kreitman’s test (MKT) was used to detect neutral and non-neutral variation in Hsp70 gen nucleotide sequences using MKT website (^{Egea et al., 2008}).

Expression of 2 variants of the Hsp70 gene (Hsp70-1 and Hsp70-2) was assessed by qPCR. First total RNA was extracted from the samples using Trizol (Invitrogen Life Technologies, Burlington, ON, Canada). RNA concentrations were determined by spectrophotometry. RNA integrity was visualized in a 1.5% agarose gel and the RNA was stored at -81 °C. The cDNA strand was synthesized from 5 μg of total RNA, with oligo dT18 and the reverse transcriptase enzyme RevertAid^TM H Minus M-MuLV, according to the manufacturer's parameters (MBI Fermentas Inc., Burlington, ON, Canada). The expression levels of the 2 Hsp70 genes were determined with 1 μl to a concentration of 10 mM/μl of the following primers designed: Hsp70-1 FI 5'TCCCTACGAAGCAGACTCAGA, Hsp70-1 RI 5’TTGTCCTTGGTCATGGCTCTC, Hsp70-2 FI 5'GCTTGAAGGGC ATAATCAGCG, and Hsp70-2 RI 5’TGGTGCTGGTACTCGTCTTTC; The expression levels of the Hsp genes were compared with the EF1α (elongation factor 1α) reference gene with the primers designed: forward 5' GAA GGA AGC CGC TGA GAT GG3', reverse 3' CGG TCT GCC TCA TGT CAC GC5’. EF1α has been reported as a suitable housekeeping gene for RT-qPCR. All primers were designed in the software Primer 3 (^{Koressaar & Remm 2007}; ^{Untergasser et al., 2012}). The qPCR was performed in a thermal cycler Bio-RAD model CFX96, using 2 μl of cDNA, 1 μl of each primer and 2 x SYBR Green Master Mix (Applied Biosystems Inc., Foster City, CA, USA) for a total volumen of 25 μl under the following run conditions: 1 cycle at 50 °C for 2 min, 1 cycle at 94 °C for 10 min, 40 cycles at 95 °C for 15 s, and at a final temperature of 60 °C for 1 min (^{Fangue et al., 2006}). Nonspecific PCR products were verified by melt-curve analysis, and representative samples were electrophoresed to verify that only a single product was present. Control reactions were carried out with no reverse transcribed RNA, to determine genomic DNA contamination. Elongation factor 1α (EF1α) was used as a housekeeping gene to normalize the expression. The level expression of Hsp genes was calculated according to E = Peff (−∆Ct), where Peff is the primer efficiency calculated using LinRegPCR (^{Ramakers et al., 2003}).

Results

The Hsp 70 gene could be amplified in both Fundulus species with a length of 578 bp with a G+C (Guanine + Cytosine) content of 51%. The sequences of F. heteroclitus (Hsp 70; accession No. DQ202279) reported by ^{Fangue et al. (2006)} was used for reference.

Genetic distances between F. persimilis and F. grandissimus sequences were low (0.002), and the divergence with the reference F. heteroclitus species was 0.013 and 0.011 respectively. We found 8 samples with wider divergence with F. persimilis 0.0031, F. grandissimus 0.0053 and F. heteroclitus 0.014 (Table 1). The average haplotype diversity was Hd = 0.996 ± 0.27, with a nucleotide diversity of Pi = 0.48.

When MKT was performed in our data, 8 mutations in nucleotide sequences where found, 3 neutral polymorphisms and 5 non-neutral polymorphisms. Variation in sequences is to low for which the Ks/Kn rate could not be calculated. All mutation is in the most variable organisms (that are more closely related to F. persimilis (genetic distance of 0.0031). The 5 non-neutral mutations are distributed in 2 positions, the first in the nucleotide position 177 that correspond to the first position of the codon 60, inducing a change in amino acid composition from Leucine to Phenylalanine; the second position in the nucleotide 387 that correspond to the first position of the codon 129, from Tryptophan to Glycine. In the first position, a non-polar aminoacid is changed to a hydrophobic one; in the second position, tryptophan the bulkiest aminoacid is changed to glycine, the smallest one.

The Hsp70-1 gene expression differed significantly from the expression of EF1α in both F. persimilis (one-way Anova, F = 547.04, p < 0.0001; fold change 1.25; Fig. 2A) and F. grandissimus (F=103.87, p < 0.0001; fold change 1.54; Fig. 2A), and it was significantly higher in F. grandissimus than in F. persimilis (F = 534.7; p < 0.0001; Fig. 2B; average fold change 20.54). The Hsp70-2 gene expression also differed significantly from the expression of EF1α in both F. grandissimus (F = 284.65, p < 0.0001; fold change 1.39; Fig. 2A) and F. persimilis (F = 174.21, p < 0.0002; 1.26 fold change; Fig. 2A), and it also was significantly higher in F. grandissimus than in F. persimilis (F = 51.85; p < 0.0001; fold change 2.89; Fig. 2B).

Figure 2 A) Fold change values of Hsp gene expression in Fundulus species with respect to EF1αgene, and B) differences in the expression between species. 

Organism of F. grandissimus has an average fold change of 8.09 for Hsp70-1 and 2.84 for Hsp70-2, all samples were collected at 33.38º C, whereas F. persimilis had lower gene expression (a fold change of 14.52 for Hsp70-1 and 25.20 for Hsp70-2) and a range of distribution into areas with temperatures lower than 32 °C. Gene expression in F. persimilis differed significantly between collecting sites ChM1 and ChL1 (F = 168.96; p < 0.0001) but not between lagoons. Organisms of F. persimilis collected in ChM1 exhibit high expression levels of Hsp70-1 and Hsp70-2 that those collected in point ChL1 (Fig. 3).

Figure 3 Expression values of Hsp gene expression with respect to the EF1α gene for F. persimilis collected at points in Chelem and Carbonera lagoons that differed in water temperature.  

Discussion

The protein encoded from the sequences belongs to the Hsp70 superfamily chaperones that help to fold many proteins. The variations detected in sequences of Hsp70 gene in Fundulus individuals of the Yucatán Peninsula, are neutral and non-neutral polymorphisms. Neutral polymorphisms occurred in the third position of the codon without affect protein structure, and non-neutral polymorphisms occurred in the first position of the codon affecting amino acid sequences. Non-neutral mutations are probably changing the protein secondary structure, since the changes in polarity and space could be affecting protein interactions. These changes could be induced by environmental conditions; influencing the adaptive capacity of the species to thermal stress.

Both, F. persimilis and F. grandissimus, respond to water temperature by the expression of Hsps, but organisms differ in magnitude expression of Hsp. Fundulus grandissimus exhibit higher expression levels than F. persimilis for both isoforms of Hsp (HSP70-1, Hsp70-2 average fold change 20.54 and 2.89 respectively). Differences in the expression between species could due to natural genetic process, characteristics of each species affecting gene expression (^{Hochachka & Somero, 2002}), however our results show differences in sensitivity to thermal stress of both Fundulus species. According our results, F. grandissimus is more sensitive than F. persimilis, similar result has been reported in Fundulus populations from the north (Fundulus heteroclitus macrolepidotus) and south (Fundulus heteroclitus heteroclitus) of USA. ^{Fangue et al. (2006)} report that organisms of Fundulus from the northern exhibit higher expression levels of Hsp70-2 between 33 to 35 ºC than organisms from southern, they conclude that northern fish are more sensitivity to thermal stress.

The magnitude of the HSPs expression in organisms is proportional to ribosomal protein damage (^{Pelham, 1984}). The Hsp70 protein families are chaperone molecules that promote the folding of newly created polypeptides and reassembly of the ribosomal proteins; they are very important in mechanisms that prevent cell apoptosis (^{Kiang & Tsokos, 1998}; ^{Sreedhar & Csermely, 2004}). Those organisms with elevated expression of HSPs are more sensitive to thermal stress and those with lower expression probably are under adaptation process (^{Dilorio et al., 1996}).

Organisms of F. presimilis exhibit a wide range in magnitude expression of Hsp (between 14.52 and 25.21 fold change) among sites (ChM1 and ChL1), elevated levels of Hsp were recorder in ChM1 point (at 31.94 ºC), 0.6 ºC hotter than point ChL1. These temperatures are reflected in the expression of Hsp, in F. persimilis specimens, when there is an increase of 0.6 ºC organisms expresses 14.52 more HSP70-1 and 25.21 more HSP70-2 (Fig. 3). Similar results were reported in F. heteroclitus, when organisms are exposed to warmer environments, induce expression levels of both Hsp but higher levels of the Hsp70-2 (^{Fangue et al., 2006}).

There are already reports of an increase in water temperature in Chelem Lagoon; the average for 1996 was 25.8 ºC with a maximum of 28.5 ºC (^{Herrera-Silveira, 1997}), in contrast with the average for 2012 -27.5 ºC with a maximum of 31 ºC (^{Tenorio-Fernandez, 2015}). The maximum during the present study was 33.39 ºC; this increase of ~5 ºC since 1996, probably play an important role in the distribution of fish species. The effect of temperature in fish distribution has been reported as a determining factor in the biogeography of the fish species in tropical costal lagoons (^{Harrison & Whitfield, 2006}; ^{Vieira & Musick, 1994}). Whereas the expression of the Hsp genes with respect to the EF1α gene and the non-neutral polymorphisms detected in the sequences of F.grandissimus and F. persimilis indicate that Hsp genes are molecular markers implied in adaptive mechanisms in both Fundulus species and have important role in thermal adaptation.

Our results demonstrate the response mechanisms of F. grandissimus and F. persimilis to thermal stress and the capacity of each species to resist elevated water temperature. Both species are sensitive to thermal stress; however, F. grandissimus seems to have a wide range of response in Hsp genes that affect their distribution in the Chelem Lagoon. Populations of F. grandissimus are generally found in areas where there is freshwater input from inland where water temperature is low with respect to the rest of the lagoon (^{Gallardo et al., 2014}). However, the period of collection of the species was in the dry season, when rainfall is low (^{Conagua, 2017}), resulting in a decrease in the contribution of ground water and higher temperatures in the lagoon. The differential sensitivity to thermal stress of both species due to expression levels of Hsp could be affecting the distribution of the species, and one of the reasons why these species are never found together in the coastal lagoons of Yucatán Peninsula (^{Álvarez, 2016}).

In order to design effective conservation programs, it is crucial to study the fish response mechanisms to variations in climatic conditions and how these mechanisms can influence the distribution of the species. Further studies should seek to detect response to thermal stress under laboratory conditions, and other molecular markers that impart resilience to fish populations. Also, we recommend to perform protein models of interaction to elucidate the effect of non-neutral polymorphisms in the structure and function of Hsp.