Hydrocarbons occurrence and transcriptomic response of oyster Crassostrea virginica from lagoons of the Southern Gulf of Mexico

The Gulf of Mexico is an important crude oil reserve worldwide, and the oyster Crassostrea virginica is an excellent candidate to study the oil spill impacts on marine invertebrates. In this work, the concentrations of polycyclic aromatic hydrocarbons (PAHs) and aliphatic hydrocarbons (AHs) from eight productive oyster areas in the Gulf of Mexico were measured on sediment, water, and tissues from C. virginica. In water, the highest AHs concentration was detected in Tamiahua (0.50 ng/mL), while for PAHs, the highest concentration was > 0.10 ng/mL in Tampamachoco. In sediment, Tamiahua and Tampamachoco lagoons had the highest AHs concentrations with values near 2.5 μg/g dry weight. Considering the PAHs, Tamiahua, Carmen, and Tampamachoco lagoons registered the highest levels, with values > 60 ng/g dry weight. In tissues from C. virginica, La Pesca, Cármen and Mecoacán presented the highest PAHs concentrations with values between 0.20 and 0.25 μg/g dry weight. Furthermore, from the molecular analysis of genes related with different phases of the xenobiotic detoxification process such as hypoxia inducible factor (hif-1a), cytochrome P450 10 (cyp10), flavin mono-oxygenase (fmo), glutathione S-transferase (gstΩ1), multidrug resistant protein (mdrd1), catalase (cat), among others, the differences between lagoons were significant (P< 0.05) and generally with Las Enramadas showing the highest expression levels. From genes evaluated in this study, hif-1a, cyp10, fmo, mdrd1 and cat presented the highest expression differentials between lagoons. The above allowed us to validate the utility of molecular markers in the assessment of the hydrocarbons effect on oyster under the conditions from the Southern Gulf of Mexico.


Introduction
The study of oil pollution in the global oceans and coastal areas confronts two opposing aspects of human activities: first, the alteration of marine and coastal ecosystems caused by the operations of extraction, refining, transport, storage, and the use of oil as the primary source of energy; and second, the undeniable need to preserve and protect marine resources for our current benefits and future generations (Suprenand et al., 2020;Da Silva et al., 2022). In this sense, the impacts of human activities on coastal ecosystems have been poorly understood and underestimated, with the widespread acceptance of the finite capacity of coastal waters to buffer the pressures of man and his activities. Furthermore, recognizing the complexity of the impact caused by them, so that all living resources or not present in the coastal zone, they cannot be considered available for indiscriminate use (Ponce-Veĺez and Botello, 2005;Magris and Giarrizzo, 2020).
In the case of crude oil, it is composed of a mixture of aliphatic hydrocarbons (AHs) and mono-and polycyclic aromatic hydrocarbons (MAHs and PAHs), with differences in density, chemical structure, and molecular weight (Hayakawa, 2018a). From these, the most toxic components are represented by PAHs because of their mutagenic and carcinogenic effects, such as benzo [a]pyrene, and they are included in most environmental monitoring programs (Fernandes et al., 2021;Wnorowski et al., 2022). In coastal areas, the oil comes mainly from runoff, sewage, oil spills, and industry. Other significant oil inputs include urban runoff in rivers and estuaries discharged into the ocean, wastewater effluents, oil spills from cargo ships, and operational discharges from commercial and recreational vessels (Fingas, 2013;Blackburn et al., 2014). Biogenic and pyrogenic ways represent other sources of AHs and PAHs components. For instance, the waxes of terrestrial plants generate AHs, and it also has been reported that alkanes and alkenes (C 15 to C 19 ) are produced by photosynthesis in the ocean (Li et al., 2010;Coates et al., 2014;White et al., 2019). For PAHs, they also produced by the incomplete production of wood (Lim et al., 2022). In the specific case of the water column and sediment (benthos), oil and gas production and extraction processes contribute significantly to the pollution around production facilities (Hayakawa, 2018b;Carpenter, 2019).
Currently, the increased energy demands of anthropogenic activities have elevated the concentrations of AHs and PAHs in the marine environment, with undesirable effects on different marine ecosystem components (El Nemr et al., 2016). The different sources of n-alkanes (AHs), including the isoprenoids pristane (Pr) and phytane (Ph), and PAHs, made challenging their characterization and identification. Furthermore, hydrocarbon components undergo various biological, chemical and physical weathering processes, adding more variables and affecting the correct hydrocarbon traceability (Leider et al., 2013). However, current methods to identify property sources of AHs and PAHs include multivariate analysis, stable isotope measurement, and diagnostic ratio assessment. In the case of diagnostic ratio, it is based on the relative abundance of the individual or specific groups of AHs or PAHs and is applied to identify hydrocarbon sources in the marine environment (Dashtbozorg et al., 2019;Wu et al., 2021). Unlike exposure tests in the laboratory where the origin of hydrocarbons is known, in the environment, the origin of hydrocarbons in the different matrixes (water, soil and tissue) is an incognita, so the application of diagnostic ratios is usual as a frame of reference (Montuori et al., 2022).
Marine invertebrates, such as oysters, often receive little attention after oil spills unless high mortality of this species or other commercial species is observed. Reports of wildlife impacts are usually limited to the most visible "charismatic megafauna," including seabirds, sea turtles, and mammals. However, oil has multiple adverse effects on marine invertebrates due to habitat degradation, accumulation of petroleum derivatives externally and suffocation of individuals, toxicity, as well as disruption of the food web (Fingas, 2013;Arnberg et al., 2018;Keesing et al., 2018;Hook et al., 2022). In this regard, different groups of marine invertebrates are exposed and respond to oil differently depending on their habitat, feeding mode, and ability to process ingested pollutants (Albers, 2003;Asif et al., 2022) and spills that are not immediately lethal may have short-or long-term impacts on behavior, reproduction, growth and development, immune response, and respiration (Takeshita et al., 2021;Wang et al., 2022). On the other hand, the extent and duration of damage caused by oil spills are further affected by the type of oil spilled and its rate of release; prevailing meteorological and oceanographic conditions; the life stage of the exposed fauna and its mobility within the habitat; substrate characteristics; spill mitigation and response activities (Hook, 2020). Our understanding of the effects of oil pollution on the structure and function of marine communities and ecosystems remains limited because they are complex, variable and infrequently studied, and can be complicated by the presence of other pollutants, such as heavy metals whose synergistic capabilities are unknown (Soares et al., 2022;Zerebecki et al., 2022).
C. virginica is a sentinel specie for the study of the pollutant's effect on aquatic invertebrates. Also, it is a conspicuous member of the invertebrates from the Gulf of Mexico. Its ecological important is based on the fact that they provide habitat for many fish and shellfish species, especially in otherwise soft-bottom environments. Their biogenic reefs supply a hard substrate for sessile invertebrates to settle on, provide refugia from predation within their complex structures, enrich the benthos with biodeposits, and provide structure and forage species for predatory fishes, many of which have value to recreational or commercial fisheries (Kennedy, 1996;Ozbay et al., 2017). Also, oysters, like all filter-feeders, may increase photosynthetically active radiation penetration by filtering phytoplankton and other particles from the water to the point where submerged aquatic vegetation can become reestablished and may help to remove excess nutrients in anthropogenically eutrophied estuaries and bays (Newell, 2004). In addition to its ecological significance, it is also an economically important species in whole the Gulf of Mexico. In the United States, the oyster's production from the coastal regions generated an income of $74 million in 2012 (Vignier et al., 2019). In Mexico, according to CONAPESCA (2014), its fishery production was of 84.5% of the total oyster production in the country (45 392 tm), and is mainly contributed by the states of Veracruz (21 719 tm), Tabasco (20 936 tm) and Tamaulipas (2 737 tm). However, the production levels of this species could be affected in the medium and long term because there are currently a large number of anthropogenic impacts that threaten the health of most ecosystems (Bravo, 2007). Within this context, and based on its physical and chemical characteristics, petroleum is a compound that could negatively affect oyster beds in the Gulf of Mexico.
The Gulf of Mexico is an ocean basin with a significant extension (~1.5 million km 2 ) and represents an extensive oil exploitation and transport area where oil spills may occur incidentally, impacting the marine ecosystem (Hernández-Guzmán et al., 2021). In previous studies with Crassostrea virginica, individuals from the south of Laguna Madre (La Pesca, Tamaulipas, Mexico) were exposed under laboratory conditions to a mixture of super-light and light oil dissolved in a solvent organic , and to light oil only dissolved in seawater to get the Water-Accommodated Fraction (López-Landavery et al., 2022). The maximum oil concentration evaluated was 200 μg/L with 14 and 21 days of exposure, respectively. Both studies focused on the modulation of biological processes through differential gene expression and functional enrichment, integrated with histology and bioaccumulation analyses. Some genes and metabolic pathways differentially expressed were associated with the different categories of the xenobiotic detoxification processes and the defense strategies against such contaminants, as was proposed by Sheehan et al. (1995) and Amiard-Triquet & Amiard (2013). As a consequence of crude oil exposure, we identified differentially expressed genes (DEGs) associated with the master regulation of the stress response, such as the hypoxia-inducible factor (hif-1a), hydrocarbon receptor (ahr) and AHR nuclear translocator (arnt). Other groups of DEGs were the cyp450 family proteins, flavin-containing monooxygenase protein (fmo), glutathione S-transferases (gstW1), multidrug-resistant protein (mdrd1) involved in the biotransformation and excretion of xenobiotics; while catalase (cat), glutathione hydrolase (ggt1) and thioredoxin reductase are involved in response to oxidative and environmental stress, respectively. Also, the crude oil exposure modulated important metabolic pathways involved in energy generation with the down-regulation of cytochrome c oxidase and the increase of metabolic rate at 200 μg/L. Specifically, in a previous study of the hydrocarbon exposure in C. virginica under laboratory conditions , 200 μg/L had a negative impact on oyster's response. For this reason, in this study, we decided to include individuals exposed at 200 μg/L for comparison purposes to have an idea of how other factors in the environment not measured here (heavy metals, plankton, pathogens, currents) modify the oyster's response to hydrocarbons presence from different sources. On the other hand, the assembly of the oyster transcriptome from the digestive gland allows us to identify and choose candidate housekeeping genes (HKGs) to execute the stability analysis necessary for the relative expression with qPCR. Based on the above context, the present work aims to evaluate, through the gene expression, the oyster's response from different lagoons of the Gulf of Mexico to the presence of hydrocarbons, using specific potential markers selected from the previous oil exposure under laboratory conditions in individuals of Crassostrea virginica. In this sense, our working hypothesis was that because hydrocarbons harm the physiological response of the oyster, higher concentrations in its soft tissue will mainly induce a greater expression of the biomarkers.

Locations of study
The study locations ( Figure 1) were selected from scientific literature reporting the presence of hydrocarbons in water, sediment, or soft tissue of the organisms in the Gulf of Mexico (Botello et al., 2014;Hernańdez-Martıńez and Ferrari, 2017). For a reliable result of the evaluation of potential marker genes in response to hydrocarbon exposure, at least 30 individuals of Crassostrea virginica were selected per sampled lagoon and 15 from a sample of individuals exposed to light hydrocarbon under laboratory conditions at the facilities of the Universidad Tecnoloǵica del Mar de Tamaulipas Bicentenario (La Pesca, Soto La Marina, Tamaulipas).

FIGURE 1
The geographical location (circles) of the lagoons from the Gulf of Mexico from where the individuals of Crassostrea virginica were obtained to evaluate the presence of hydrocarbons in water, sediment, and organisms; and gene expression of the potential marker genes in the digestive gland of C. virginica.

Assessment of aliphatic hydrocarbons (AHs) and polycyclic aromatic hydrocarbons (PAHs)
The procedures used to get the hydrocarbon extracts in soft tissue from individuals of C. virginica sampled in the different lagoons of the Gulf of Mexico, as well as from water and sediment, were based on Murphy et al. (2012) and Glaser et al. (1981). In the case of tissue samples, each lagoon consisted of three independent pools (n = 3), with each pool consisting of 15 organisms. In the case of sediment, three replicates (n = 3) were used, and three volumes (n = 3) of 1.8 L in the case of water. All samples were stored in precleaned amber glass jugs. Water samples were saturated with a mercuric chloride solution and refrigerated at 4°C until analysis. Loṕez-Landavery et al. (2019) reported more details about the processing of the samples. Organisms from all lagoons were included in the assessment of hydrocarbons, except those from Las Enramadas, which were used only in the expression analysis.

Diagnostic ratios of n-alkanes (AHs) and PAHs
Diagnostic ratios were calculated to identify possible sources of hydrocarbons in the different matrices (water, sediment, and soft tissue) collected from the Gulf of Mexico's coastal lagoons. Carbon preference index (CPI), the terrigenous/aquatic ratio (TAR), the ratio between lower-molecular weight and higher-molecular weight compounds (LMW/HMW), and the natural n-alkane ratio (NAR)) were calculated to n-alkanes. LMW indicates the sum of n-alkanes from n-C 9 to n-C 32 , while HMW indicates the sum of n-alkanes from n-C 24 to n-C 36 . Equations described to continuation were reported previously by Hernańdez-Guzmań et al. (2021). CPI = o (C 23 − C 31 ) odd n−alkanes + o (C 25 − C 33 ) odd n−alkanes 2 o (C 24 − C 32 ) even n−alkanes TAR = (C 27 + C 29 + C 31 ) (C 15 + C 17 + C 19 ) In the case of PAHs, the most used diagnostic ratios are phenanthrene/anthracene (Phe/Ant), fluoranthene/pyrene (Fla/Pyr), and the ratio between lower-molecular weight and high-molecular weight compounds (LMW/HMW) (Nassar et al., 2011;He et al., 2014). LMW/HMW indicates the concentration ratio of 2-3 ring PAHs to 4-6 ring PAHs.

Gene expression analysis
Total RNA extraction and cDNA synthesis About 30 mg of the digestive gland from the collected organisms were used for RNA extraction. All samples were fixed with RNAlater (Ambion, Austin, TX, USA) and stored at -80C until RNA purification. Following the manufacturer's recommendations, the total RNA was extracted with the RNeasy Mini Plus ® kit (Qiagen, Valencia, CA, USA). The total RNA was eluted with 50 mL of DNasefree water and quantified by absorbance using a Nanodrop ND-2000 spectrophotometer (Thermo Scientific, USA). The RNA integrity was confirmed by electrophoresis using a 2% agarose gel.
RNA samples were treated with DNase (DNase I, RNA Free, Qiagen) to remove the residual genomic DNA and thus prevent it from interfering with the evaluation of the expression of potential marker genes. After DNA elimination, the total RNA was quantified with a Qubit fluorometer ® 3.0 using a Qubit ™ RNA Broad Range assay (Life Technologies, USA), and its integrity was assessed through electrophoresis. Effective DNase treatment was confirmed by verifying no-amplification after PCR when using the housekeeping gene rpl19 (Ribosomal Protein subunit Large 19) ( Table 1).
For cDNA synthesis, we proceeded to synthesize 1 μg of total RNA treated with DNase from each sample. For this, the Impron-II Reverse Transcriptase (Promega) Kit was used. cDNA synthesis was confirmed through PCR using the housekeeping gene rpl19. The correct amplification of the rpl19 gene showed a single amplification band when it was analyzed in 2% agarose electrophoresis.
Real-time PCR analysis (qPCR) for validation of potential marker genes Selection of housekeeping genes (HKGs), genes of interest (GOIs) and efficiency assessment HKGs are necessary for essential functions in the organism. They should not vary significantly regardless of the conditions to which the organism is exposed. The HKGs and GOIs (Table 1) used in this study were selected from RNA-Seq data of C. virginica exposed to a mixture of light and super-light oil and light oil only under laboratory conditions Loṕez-Landavery et al., 2022). For the initial standardization of the selected genes, endpoint PCR was used. After determining the optimum annealing temperature (T A ) of each HKG and GOI primers, we calculated their efficiency through the standard curve using the equation E= (-1 + 10 (-1/slope) ). For this, a sample containing the cDNA of all individuals of C. virginica from the different lagoons of the Gulf of Mexico was prepared and used as a template for qPCR reactions. The qPCR was conducted by triplicate on a CFX96 Real-Time PCR Detection System (Bio-Rad). Each qPCR contained 2.5 mM MgCl 2 , 1× reaction buffer, 0.2 mM dNTP Mix (Promega, WI, USA), 1×EvaGreen fluorescent dye (Biotium, CA, USA), 0.2μM of each forward and reverse primers, 5 U/μL of Accustart Taq DNA Polymerase (Quanta Biosciences), and 3 ml of a cDNA dilution 1:5 in a final volume of 10 ml. Amplification conditions were: denaturation at 95°C for 2 min, followed by 40 cycles for 30 s at 95°C, 30 s at 60°C, 20 s at 72°C acquiring the fluorescence at this step.
In the end, a dissociation ramp from 65°C to 95°C (steps of 0.5°C/s; acquiring the fluorescence at the end of each step) was performed to ensure the absence of artifacts and primer dimers, verifying the specificity of the PCR product.

Stability analysis of HKGs
Once the reaction efficiency of the HKGs was determined, the next step was to select the most stable genes for the conditions of the C. virginica individuals sampled in the Gulf of Mexico. For this, between 3 and 6 biological replicates per lagoon were analyzed individually (Table 2) using the same amplification conditions established in the standard curve. The goal of this analysis was to ensure the less variability possible, in terms of expression, between the individuals sampled from the different lagoons. The average Cq values of the technical replicates were analyzed with the geNorm (Vandesompele et al., 2002) and NormFinder (Andersen et al., 2004) programs to determine the most stable genes between the lagoons' conditions.

Relative expression of GOIs
The relative expression of the potential marker genes was carried out using the most stable HKGs obtained with geNorm and NormFinder and cDNA from the biological replicates of the Crassostrea virginica individuals sampled in the lagoons indicated above. The relative expression values for each of the genes evaluated in this study were performed with the qBase plus ® program (Biogazelle, Belgium), using the geometric mean method (Hellemans et al., 2007).

Statistical analysis
The assumptions of normality and homoscedasticity were evaluated with the Shapiro-Wilk and Levene tests, respectively. A one-way ANOVA was used to determine statistical differences in the presence of hydrocarbons between the lagoons for each matrix (water, sediment, and soft tissue). Also, a one-way ANOVA was used to detect significant differences in the relative expression values between lagoons, with a significance value of 0.05. For both cases (bioaccumulation and gene expression analyses), when there were significant differences between lagoons, a Tukey's honestly significant differences (HSD) test was applied with P< 0.05. The analysis was carried out in R version 3.6.3 (R Core Team, 2019).

Hydrocarbons content in water and sediment from lagoons in the Gulf of Mexico
The chromatography analysis confirmed the presence of aliphatic hydrocarbons (AHs) and aromatic hydrocarbons (PAHs) in the water and sediment of the different lagoons of the Gulf of Mexico. For water, the highest concentrations of AHs (F = 70.14, P*** =8.14e-11) were found in Tamiahua and Laguna de Teŕminos with 0.50 and 0.43 ng/ mL, respectively (Supplementary material (SM), Figure S1A). The lagoons with the lowest levels of AHs in water were Caŕmen, La Pesca, and Mecoacań, with values below 0.05 ng/mL. In the PAHs case (F = 64.97, P*** = 2.86e-08), the highest concentrations were registered in Tampamachoco and Laguna de Teŕminos, with values greater than 0.06 ng/mL ( Figure 2A). La Pesca and Mecoacań registered the PAHs' lowest levels with 0.002 and 0.009 ng/mL, respectively. The most abundant PAH components in water were phenanthrene > pyrene > fluoranthene.

AHs and PAHs in soft tissues of C. virginica from lagoons in the Gulf of Mexico
The levels of AHs (F = 15.34, P*** = 5.23e-06) in the soft tissue of C. virginica were higher in the organisms of the natural environment than those exposed to hydrocarbons under laboratory conditions (data not shown). Within the sampled lagoons, Pajonal, Tampamachoco, and Alvarado had the highest concentrations with values above 10 mg/g dry weight (SM, Figure   S1C). In comparison, La Pesca recorded the lowest value with 3.62 mg/g dry weight, although it did not present significant differences with Tamiahua, Laguna de Teŕminos, Caŕmen, and Mecoacań (P > 0.05). In the case of PAHs (F = 7.45, P*** = 4.56e-04), the highest concentrations were obtained for La Pesca, Caŕmen, and Mecoacań, with values between 200 and 250 ng/g dry weight. A previous study reported a concentration of PAHs of 230 ng/g dry weight in the soft tissue of organisms collected in La Pesca . On the other hand, the lagoon with the lowest PAH concentration was Alvarado in the state of Veracruz, with a value of 65 ng/g dry weight ( Figure 2C). The most abundant PAH component detected in the soft tissue of individuals from the Gulf of Mexico was perylene, except for Alvarado. Other abundant PAH components were phenanthrene and naphthalene. PAHs levels on a volume and dry weight basis detected in water (A), sediment (B) and soft tissues (C) of Crassostrea virginica collected from different lagoons of the Gulf of Mexico. Values represent mean ± standard deviation. Equal letters indicate no statistical difference between locations based on ANOVA one-way (P< 0.05). ND, Non-detected or under limit of detection.

Diagnostic ratios of n-alkanes (AHs) and PAHs
For n-alkanes (Table 3), the CPI values were higher than one (1.06 -7.01), except for samples water in Tampamachoco (0.63) and Pajonal (0.71); and in soft tissues of individuals from Laguna de Teŕminos (0.46). In the case of TAR, sediments from different lagoons predominantly had values higher than one (1.36 -15.30), while NAR had values lower than one in all matrices collected along the Gulf of Mexico. Considering the lower-and higher-molecular weight compounds, water and sediment had values predominantly lowers than one, while soft tissues had values higher than one (1.56 -217.98), except for Tampamachoco (0.48).
Based on diagnostic ratios for PAHs (Table 3), Phe/Ant had values higher than one for all matrices (1.39 -27.16), except for water at Pajonal (0.01). In contrast, the Fla/Pyr ratio values were lower than one in water and soft tissue, except for Tamiahua (matrix = soft tissue, value = 1.77). Interestingly, PAHs of 2-3 rings predominated on PAHs of 4-6 rings in water. However, that relationship was inverse in sediment and soft tissues, with a predominance of PAHs of 4-6 rings, probably due to the sediment's chemical characteristics and the physiology of oysters presenting a low rate of elimination of PAHs.

Selection of HKGs and efficiency
The optimum Tm of the HKGs was generally around 60°C. Also, based on the standard curves (SM, Figure S2), the efficiency values fluctuated between 85 and 100.7% (Table 1). The dissociation curve analysis confirmed the specificity of the primers. From HKGs evaluated (Table 1), N-acetyltransferase domain 1 had (natd1) the lowest efficiency. In contrast, the highest value was obtained with the ribosomal protein S24 (rps24).

Stability analysis of HKGs
The stability values (M) reported by geNorm, indicated that the most stable genes under experimental conditions were tuba, rpe, and ef1a (M< 0.5), while the least stable were rps24, natd1, and ks6a1 ( Figure 3A). In the case of NormFinder, the most stable genes were natd1, ef1a, and rpe with M< 0.017 ( Figure 3B). On the other hand, the paired variation analysis indicated that the optimal number of HKGs for correct expression analysis is a minimum of 3 (M< 0.15). However, M = 0.20 has been reported as the maximum threshold ( Figure 3C).

Amplification of GOIs and efficiency
Like HKGs, potential GOIs had an optimal melting temperature of around 60°C. Also, the dissociation curve confirmed the specificity of each primer, and the standard curves (SM, Figure S3) indicated efficiency values between 82.0 (Metalloproteinase inhibitor) and 104.8% (Flavin monooxygenase) ( Table 1).

Relative expression analysis through qPCR
The relative expression analysis of the GOIs was carried out using the previously selected HKGs (tuba, rpe, and ef1a) and the biological replicates of the Crassostrea virginica individuals sampled in the Gulf of Mexico lagoons (Table 2). In this section, results are described considering three gene categories: master regulators of the stress response (hif-1a, ahr, and arnt), biotransformation and excretion of xenobiotics (cyp302a1, cyp10, gstW1, fmo, sult4a1, and mdrd1), and response to oxidative stress and redox process (cat, ggt1, ttr3, and meti4).

Master regulators of the stress response
For this group of genes (Figure 4), hypoxia-inducible factor (hif-1a), hydrocarbon receptor (ahr), and AHR nuclear translocator (arnt) showed significant differences between lagoons (hif-1a: P***< 0.0001, ahr: P* = 0.037, arnt: P* = 0.038). The lowest levels of hif-1a occurred in La Pesca. In this context, all the evaluated lagoons presented significant differences concerning La Pesca, while the highest levels of expression were shown in Las Enramadas and Mecoacán. For organisms exposed to hydrocarbons under laboratory conditions (200 μg/L during 21 days), differences in expression were detected after 21 days of exposure (T3). In the case of ahr, the lowest level of expression occurred in organisms exposed to hydrocarbons at 14 days (T2), while Las Enramadas presented the highest levels of expression. For arnt, Caŕmen was the lagoon with the lowest expression level, while Las Enramadas again showed the highest expression levels.

Biotransformation and excretion of xenobiotics
The process of biotransformation and excretion of xenobiotics is divided into three phases. During phase I, the cytochrome P450 (CYP) family has been reported to participate in this response. In this study, the expression levels for two CYP members ( Figure 5) were Frontiers in Marine Science frontiersin.org significant between lagoons (cyp302a1: P***< 0.0001, cyp10: P** = 0.0002). The lowest expression values for cyp302a1 were obtained in Laguna de Teŕminos, Mecoacań, Tampamachoco, and Caŕmen, significantly different from those of La Pesca (P*< 0.05). In the case of organisms exposed under laboratory conditions, the differences concerning La Pesca occurred at 14 days of exposure (P*< 0.05). For cyp10, the highest levels of expression occurred in organisms exposed under laboratory conditions, and the values were significantly different from those found in Laguna de Teŕminos (P*< 0.05).

FIGURE 5
Relative expression of genes involved in the phase I (Up) and phase II (low) of the xenobiotic's detoxification processes. Values represent mean ± standard deviation. Equal letters indicate no statistical difference between locations based on ANOVA one-way (P< 0.05).

FIGURE 4
Relative expression of genes related with the master regulation of the stress. Values represent mean ± standard deviation. Equal letters indicate no statistical difference between locations based on ANOVA one-way (P< 0.05). and fmo (P***< 0.0001) had significant differences in expression levels when comparing lagoons. In the case of gstW1, the highest levels of expression occurred in organisms exposed to hydrocarbons for 21 days and in Las Enramadas, whose values were significantly different from those of La Pesca ( Figure 5). Likewise, the lowest expression values were presented in Mecoacań, La Pesca and organisms exposed to hydrocarbons for 14 days. An interesting behavior of expression levels was observed for fmo ( Figure 5), where the lowest values were presented in Las Enramadas and La Pesca. In contrast, the highest values were obtained in Tampamachoco, Caŕmen, and Mecoacań. The values found in Tampamachoco, Caŕmen, and Mecoacań were significantly different from those of La Pesca, Las Enramadas, and Laguna de Teŕminos (P*< 0.05). The expression values for sult4a1 not shown significant differences (P = 0.806, data not shown).
During phase III, multidrug resistance proteins are associated with detoxification processes. In the case of mdrd1 (P***< 0.0001), the highest levels of expression occurred in the population of Las Enramadas and organisms exposed to hydrocarbons for 14 days (T2_200). The expression values of Las Enramadas and T2_200 were significantly different from those found in Laguna de Teŕminos, Mecoacań, Tampamachoco, and Caŕmen (P*< 0.05, Figure 6). So far, it can be seen that most of the genes evaluated are down-regulated concerning the population of Las Enramadas, which could indicate a chronic condition due to the presence of contaminants.

Response to oxidative stress and redox processes
Among the genes associated with oxidative stress, we have catalase (cat) and glutathione hydrolase (ggt1). For both genes, the expression values between lagoons shown significant differences (cat: P***< 0.0001, ggt1: P***< 0.0001). In the case of cat, the highest levels of expression occurred in Tampamachoco, Caŕmen, and Las Enramadas, being, in turn, significantly different concerning La Pesca (P*< 0.05). Of all the lagoons, La Pesca presented the lowest level of expression. For ggt1, Las Enramadas and Tampamachoco again showed the highest expression levels, significantly different from those found in Mecoacań and Laguna de Teŕminos (P*< 0.05). In La Pesca, it only presented significant differences with Las Enramadas (P*< 0.05). For the trr3 gene (P***< 0.0001), which is involved in regulating redox processes and is essential for correct signaling, the highest levels of expression occurred in organisms exposed to hydrocarbons for 21 days (T3_200) and in Las Enramadas. The expression levels of these lagoons were significantly different from the values found in La Pesca, Mecoacań, and Tampamachoco (P*< 0.05). In general, La Pesca presented the lowest level of expression for trr3. Expression levels for meti4, associated with response to environmental stress, did not shown significant differences (P = 0.057, data not shown) Relative expression of genes involved in the phase III of the xenobiotic's detoxification processes (mdrd1), oxidative stress (cat and ggt1), and environmental stress. Values represent mean ± standard deviation. Equal letters indicate no statistical difference between locations based on ANOVA one-way (P< 0.05).

Hydrocarbons bioaccumulation in coastal lagoons of the Gulf of Mexico
The coastal lagoons used in the present work were selected from previous studies that have reported the presence of hydrocarbons in water, sediment, or oyster tissue of the Gulf of Mexico, from Tamaulipas in the north to Campeche in the southeast (Botello, 2005;Montaño-Vera et al., 2017;Botello et al., 2019). In water, aliphatic hydrocarbons (AHs) and aromatic hydrocarbons (PAHs) were found in all the lagoons included in the present study; with the highest levels found in Tamiahua-Laguna de Teŕminos and Tampamachoco-Laguna de Teŕminos, respectively. Celis et al. (1985), in a review of different studies on the concentration of hydrocarbons in water from various bodies of water in the Gulf of Mexico, found that Laguna de Teŕminos and the Tuxpan river presented the highest concentrations (20-48 mg/L). Laguna de Teŕminos, located in the state of Campeche, is a lagoon complex situated in the area of marine oil platforms and receives discharges of particulate and dissolved organic matter from the Grijalva-Usumacinta river system (Paez-Osuna et al., 1987;Botello, 2005). In the case of the Tamiahua and Tampamachoco lagoons, in addition to being in an area of high oil exploitation and having influence from the Tuxpan river and, to a lesser extent, the Panuco river, they have a high annual maritime traffic (Botello and Calva, 1998). The factors mentioned above could explain comparatively the concentrations of hydrocarbons found in water in Laguna de Términos, Tamiahua, and Tampamachoco concerning other coastal lagoons of the Gulf of Mexico.
Sediments serve as depositories of several sources of organic matter and hydrocarbons, and they can be used to assess the historical impact on the environment (Romero et al., 2015). Likewise, the type of organic matter and hydrocarbon that accumulates will depend on the primary productivity of the area, physical and chemical factors, and the nature of the sediment (Waterson and Canuel, 2008). In this study, Tamiahua, Tampamachoco, and Laguna de Teŕminos had the highest levels of AHs, while Tamiahua, Cármen, and Tampamachoco had the highest concentrations of PAHs. In the case of Tamiahua, Tampamachoco, and Laguna de Teŕminos, which also had the highest levels of hydrocarbons in water, the levels of hydrocarbons found in the sediment would reflect the effect of the oil activity carried out over many years in this area of the country and the fluvial influence mentioned above. In addition, in front of the Laguna de Teŕminos (~80 km) are the PEMEX oil exploitation platforms that correspond to the Sonda de Campeche, whose contributions could reach the coast by the action of currents and winds (Botello, 2005). From a quantitative approach, Gogou et al. (2000) reported that areas with low concentrations of AHs had values between 562 and 5,607 ng/g, while, for PAHs, the range fluctuated between 14.6 and 158.5 ng/g. Based on the above ranges and the values found in the present study for both AHs and PAHs, it could be suggested that the coastal lagoons of the Gulf of Mexico have low concentrations of hydrocarbons in the sediments.
The degree of accumulation of pollutants, such as oil, in a marine organism will depend on the biological availability of its different fractions, the time of exposure to which the organism is subjected, and the ability of the organism to metabolize it (Van der Oost et al., 2003). In the present study, the organisms that presented a higher level of bioaccumulation of AHs corresponded to Tampamachoco, Pajonal, and Alvarado. In the case of PAHs, the organisms of La Pesca, Caŕmen, and Mecoacań had the highest concentrations. Pajonal is part of the Caŕmen-Pajonal-Machona lagoon complex located in an area of exploration, drilling, and transport of hydrocarbons (Gold-Bouchot et al., 1997), so it would reflect the concentrations of AHs found in this lagoon. The Alvarado lagoon, in addition to being located in an area of oil exploitation, has the influence of the Coatzacoalcos river, which would be an additive factor in the origin of the AHs. In the specific case of La Pesca, which presented low levels of AHs and PAHs in water and sediment, the values of PAHs in the organisms suggest the idea of the presence of higher concentrations concerning organisms from other lagoons with higher levels of AHs and PAHs found in both water and sediment. The above would seem like a discrepancy since organisms living in areas with higher levels of hydrocarbons would be expected to have higher levels of bioaccumulation in their tissues. However, Gobas et al. (1999) indicated that invertebrates have less ability to metabolize xenobiotics than vertebrates, and it could result in relatively high tissue bioaccumulation even when the contaminant has disappeared from the environment. The scenario described above is one of the main causes of why invertebrates are frequently used as indicators of contamination by organic components (Coates and Söderhäll, 2021). In the case of Mecoacań, the levels of PAHs found could be due to its proximity to the Dos Bocas oil terminal, cataloged within the top three either by traffic of oil vessels or by volume of crude oil shipped (INEGI, 2010).

Diagnostic ratios of n-alkanes and PAHs
In general, diagnostic ratios were used to identify the possible sources of n-alkanes and PAHs in the different matrices (water, sediment, and soft tissue) from the lagoons of the Gulf of Mexico. The CPI values indicate the possible predominance of odd n-alkanes, even n-alkanes, or not predominance. CPI values near one suggest no preference among odd and even n-alkanes and are associated with petrogenic inputs or anthropogenic sources (Sainakum et al., 2021). In the present study, CPI values from water (0.63-1.72) and soft tissues (0.46-1.92, except Mecoacań) suggest a predominant petrogenic contribution for n-alkanes. In individuals from Carmen, the CPI value (1.92) suggests a contribution from mixed sources, with the predominance of biogenic sources, probably phytoplankton (El Nemr et al., 2016;Sainakum et al., 2021). In contrast, for individuals from Mecoacań, the high CPI value (3.21) suggests the presence of nalkanes from terrigenous sources (Leider et al., 2013), typically from C3 plants, important in the deposition of organic matter in estuaries (Lu et al., 2022).
In petroleum hydrocarbons, the NAR ratio is close to zero and close to one for marine or higher terrestrial plants (El Nemr et al., 2016). In this study, except for La Pesca, the NAR ratio was close to zero in water and soft tissues, ranging from 0.02 to 0.23 in absolute value. The above values should indicate a substantial contribution of petrogenic sources in water and soft tissues collected from the Gulf of Mexico lagoons. In contrast, the NAR value in La Pesca (-1.00) suggests a strong contribution of hydrocarbons from higher terrestrial plants or marine plants. In sediment, NAR ranged from 0.31 to 0.67, indicating a moderate contribution of petrogenic sources or even a mixture of hydrocarbons from high vascular marine or terrestrial plants (Hernańdez-Guzmań et al., 2021). On the other hand, the TAR ratio helps to differentiate the terrestrial contribution from aquatic sources of n-alkanes based on long-and short-chain of odd compounds (Rielley et al., 1991;Jafféet al., 2001). In this study, TAR values from sediment were greater than one (1.36-15.30), suggesting a predominant contribution of n-alkanes from the debris of higher terrestrial plants (Bourbonniere and Meyers, 1996). In contrast, the tendency of TAR values from water and soft tissues was lower than one, suggesting that the main contribution of nalkanes was from aquatic sources (Keshavarzifard et al., 2020). LMW/HMW values lower than one usually indicates n-alkanes produced from sources such as higher plants, marine animals, and sedimentary bacteria (Wang et al., 2006). Also, values close to one and greater than two suggest n-alkanes produced from petroleum and plankton and the presence of fresh oil input, respectively (Commendatore et al., 2000;Vaezzadeh et al., 2015). In this study, the LMW/HMW values in soft tissue (1.56-217.98, except for Tampamachoco) suggest that the predominant sources of n-alkanes are petroleum and fresh oil input. In contrast, LMW/HMW values in sediment were lower than one (0.08-0.46), reflecting the n-alkanes sources from higher plants, marine animals, and sedimentary bacteria. On the other hand, isoprenoids such as pristane (Pr) and phytane (Ph) are derived from the diagenetic processes of phytol, and they are not primary components in the terrestrial biota (Peters et al., 2005;Babcock-Adams et al., 2017). Furthermore, diagnostic ratios based on Pr and Ph have been used as markers of biodegradation for oil spills (Hernańdez-Guzmań et al., 2021). In this study, Pr/Ph values in soft tissues (data not shown) were close to one or slightly greater than one (0.90-1.67), except for Tamiahua (0.42). These values suggest that the predominant source of Ph is crude oil and corroborate the Ph nature because it is not commonly produced by biogenic sources (Caro-Gonzalez et al., 2020).
Diagnostic ratios were also used to assess the sources of PAHs. PAHs are persistent organic pollutants with pyrogenic, petrogenic, and biogenic sources (Liu et al., 2009;Kao et al., 2015;Vaezzadeh et al., 2015). Petroleum often contains PAHs that are thermodynamically more stable than others. For instance, Phenanthrene (Phe) is more stable than Anthracene (Ant), and the Phe/Ant ratio is high. Phe/Ant values from water in this study were greater than ten (12.13-27.16), except for Pajonal (0.01). The above suggests that PAHs source from the water in most of the coastal lagoons of the Gulf of Mexico is mainly petrogenic (Soclo et al., 2000). In contrast, sediment from all lagoons sampled in this study had Phe/ Ant values lower than ten (1.39-6.51), suggesting a pyrogenic source for PAHs. Phe/Ant values for soft tissues in this study focused on bivalves, suggesting petrogenic and pyrogenic sources for PAHs. Oysters from Alvarado, Pajonal, and Mecoacań would be exposed to petroleum hydrocarbons, while oysters from La Pesca, Tamiahua, Tampamachoco, and Caŕmen would be impacted by combustion residues (Ünlü and Alpar, 2006). Another ratio for PAHs source identification is Fluoranthene/Pyrene (Fla/Pyr). Fla and Pyr are usually the most abundant PAHs in pyrolytic sources, with ratios greater than 1 indicating pyrolytic origin and ratios lower than one associated with petrogenic sources (Magi et al., 2002). In this study, Fla/Pyr values of water (0.24-0.97) and soft tissues (0.35-0.95, except for Tamiahua) suggest a petrogenic source for PAHs. More specifically, Fla/(Fla+Pyr) values (data not shown) suggest that the origin of PAHs from water in Mecoacań and PAHs of individuals from La Pesca, Pajonal, and Laguna de Teŕminos is from petroleum combustion, while for other lagoons the PAHs source is petroleum (El Nemr et al., 2016). In sediment from Tamiahua and Caŕmen and individuals from Tamiahua, the Fla/Pyr values were higher than one and suggested a pyrolytic origin of PAHs, probably from coal combustion and wood (Yunker et al., 2002). LMW/HMW, defined as 2-3/4-6 ring PAHs, also indicates sources of PAHs from petrogenic and pyrogenic origins (Lu et al., 2022). PAHs from a petrogenic origin consist mainly of LMW (2-3 rings), while HMW (4-6) are predominant in PAHs with pyrogenic origin (Dashtbozorg et al., 2019). In this study, LMW/HMW values from the water of all lagoons sampled and individuals from Alvarado and Pajonal were greater than one and suggested that the origin of PAHs is petrogenic. Furthermore, in the sediment of all lagoons and individuals of La Pesca, Tamiahua, Tampamachoco, Carmen Mecoacań, and Laguna de Teŕminos, LMW/HMW were lower than one, and it indicates a pyrogenic origin of PAHs.

Hydrocarbons expression profile of potential marker genes
Master regulators of the stress response Three components of the family of master stress regulators were evaluated as possible markers in the response of oysters under the presence of hydrocarbons in the natural environment: HIF-1a, AHR, and ARNT. Although HIF-1a plays an essential role in molecular adaptation to hypoxia (Majmundar et al., 2010), functioning as a transcription factor, it has also been linked to the development of cancer by regulating the expression of oxygen-dependent genes (Fraga et al., 2009). The molecular expression levels of hif-1a found in the present work indicated that Las Enramadas and Mecoacań presented the highest expression levels. In the case of Las Enramadas, initially proposed as a second negative control, it turned out to be a lagoon where several genes had high expression levels. Las Enramadas is a small fishing community located north of La Pesca in the state of Tamaulipas and front of the Cinturoń Plegado de Perdido area. Although there are not many studies carried out in this area of the state of Tamaulipas, it is suggested that, as this area is rich in natural gas, extraction through gas pipelines could be impacting the biota present in this part of the Laguna Madre lagoon system. In the case of Mecoacań, as mentioned above and as has been reported in other studies (Botello, 2005), the presence of the Dos Bocas oil terminal has had an impact on the presence of hydrocarbons in water, sediment, and benthic organisms such as C. virginica.
AHR belongs to the family of PAS transcription factors responsible for sensing alterations in the environment related to circadian rhythm, oxygen gradients, or the presence of xenobiotics (Gu et al., 2000), and together with ARNT, are part of the AHR metabolic pathway. In the specific case of bivalves, most studies have been focused on describing the structure of AHR and ARNT , so there are few expression studies. In the present work, the organisms with the highest expression levels for both genes were those from Las Enramadas, Mecoacań, and direct exposure to light hydrocarbon for 21 days at 200 mg/L. The strong direct correlation between ahr and arnt expression levels (r = 0.5, P-value = 6.52E-05) in these lagoons corresponds to their mode of action within the AHR pathway. ARNT is needed to assemble the active heterodimeric complex that is translocated in the nucleus to promote the expression of genes related to cytochromes P450 or CYP (Mandal, 2005). Furthermore, our results agree with the differences found in the expression of ARNT in the digestive gland of individuals of Crassostrea virginica collected in Fort Morgan (USA), from the beginning of the incident of the Deep Water Horizont and for one year (Jenny et al., 2016); and with the data reported by Chatel et al. (2014) who reported that the expression of AHR was higher in individuals of the zebra mussel Dreissena polymorpha collected from areas impacted by the presence of hydrocarbons. On the other hand, and contrastingly, the results of the present study differ from those recently reported by Wang et al. (2020). They reported inhibition of ahr expression levels in the digestive gland of Ruditapes philippinarum exposed for five days to Benzo[a]pyrene under laboratory conditions. These differences could be associated with the type of compound used (individual vs. mixture), management of organisms (laboratory vs. natural environment), and the type of response assessed (acute vs. chronic).

Biotransformation and excretion of xenobiotics
Cytochromes P450 or CYP are the enzymes involved in the initial process of detoxification of hydrocarbons and xenobiotics in general (Amiard-Triquet et al., 2011). PAHs are metabolized by the CYP family of enzymes generating two types of PAHs, phenolic and diolic, depending on the chemical groups that adhere during this process. If these PAHs do not follow their metabolization process correctly, they can act as endocrine disruptors, increase the production of reactive oxygen species, or have carcinogenic or mutagenic properties (Motoyama et al., 2009;Tsay et al., 2013). The present study evaluated the expression of cyp302a1 and cyp10. La Pesca presented the highest levels of cyp302a1, while lagoons impacted by the presence of hydrocarbons such as Tampamachoco, Caŕmen, Mecoacań, and Laguna de Teŕminos had low levels of expression. Although CYP302A1 has been related to a lesser degree to detoxification processes, it has other functions, such as the development of the central nervous system, hormone metabolism, and participation in steroid activity (Petryk et al., 2003). This variety of functions, especially that related to the action of steroid hormones, would explain the expression levels of cyp302a1 found in La Pesca. In previous exposure results with La Pesca organisms and under laboratory conditions, we found that the presence of hydrocarbons decreased the production of steroid hormones related to reproduction compared with organisms from the same lagoon but without exposure to hydrocarbons (Tapia-Morales et al., 2019). Within the members of the superfamily of CYPs, those corresponding to groups 1 and 2 have been recognized for having a central role in the detoxification of hydrocarbons (Schlezinger et al., 2006;Zhang et al., 2012). In the case of cyp10, a cyp1-like, the highest expression levels were obtained in organisms exposed to 200 mg/L for 14 and 21 days, respectively. Although the differences were significant only with Laguna de Teŕminos, mainly due to the variance as a result of the low number of samples, our results agree not only with the primary function of CYP1 in the detoxification of hydrocarbons but also corroborate the results reported in Crassostrea brasiliana and C. virginica that showed high expression levels of cyp1 and cyp1-like in organisms exposed to PAHs or from areas impacted by hydrocarbons (Lüchmann et al., 2015;Jenny et al., 2016). The above suggests that cyp10 would be a good candidate as a molecular marker to evaluate the presence of hydrocarbons in the different coastal lagoons of the Gulf of Mexico.
In phase II, PAH derivatives are conjugated with enzymes such as glutathione S-transferases (GST) and monooxygenases (FMO). Conjugation not only makes the derivatives more polar but facilitates their detoxification in bivalves through pathways such as sulfotransferases (SULT) and UDP-glucuronosyltransferases (UGT) (Wu et al., 2013;Vidal-Liñań et al., 2014). In the case of gstW1, the organisms collected towards the end of the hydrocarbon exposure period and those from Las Enramadas presented the highest expression levels. In exposed organisms, there was a positive correlation with cyp10 expression. Studies reporting the same expression pattern in bivalves collected from hydrocarboninfluenced areas or exposed under laboratory conditions include Perna perna and R. philippinarum (De Luca-Abbott et al., 2005), C. virginica (Jenny et al., 2016), Venerupis philippinarum  and C. brasiliana (Lüchmann et al., 2011). In the case of fmo, the organisms that presented a more significant induction were those of Tampamachoco, Caŕmen, and Mecoacań. FMO belongs to a family of monooxygenases involved in the oxidation of sulfur-containing compounds. Previously, we have described that areas such as Tampamachoco, Caŕmen, and Mecoacań are influenced by oil activity or have a synergistic influence on river bodies whose effect is reflected in any of the three components of coastal lagoons such as water, sediment, or organism. Although few studies had reported the expression of fmo in bivalves exposed to hydrocarbons, those carried out in the digestive gland of Crassostrea gigas (Schlenk and Buhler, 1989;Boutet et al., 2004) had found higher expression levels of fmo, not only confirming the expression pattern found in the present study but suggesting its usefulness as a good marker to assess the presence of hydrocarbons present in water bodies.
Phase III of the xenobiotic detoxification process releases metabolites generated from PAHs through urine and substances such as bile or its equivalents (Bekki, 2018). Within this phase, multidrug resistance proteins play an essential role as transporters, and MRP1 or P-glycoprotein has been partially characterized in C. virginica (Ivanina and Sokolova, 2008). Of the organisms collected in the different coastal lagoons of the Gulf of Mexico, those from Las Enramadas and the exposure test under laboratory conditions (14 days at 200 g/L) presented the highest expression levels for mdrd1. In Las Enramadas, and as mentioned above, it is very likely that gas reserves and their exploitation through gas pipelines are having an impact on organisms, so this area needs to be studied more frequently to determine the actual impact not only on organisms but also on the water column and sediments. In organisms exposed under laboratory conditions, previous studies have found that the highest expression levels for GOIs in detoxification processes and oxidative stress occur more consistently at 14 and 21 days of exposure (Tapia-Morales et al., 2018;Loṕez-Landavery et al., 2019). Similar results associated with the up-regulation of a multidrug resistance protein (mrp3) were reported in the digestive gland of C. virginica collected from Fort Morgan and Orange Beach three months after the incident in the Deep Water Horizon (Jenny et al., 2016).

Response to oxidative and environmental stress
The reactive oxygen species (ROS) are compounds generated as part of detoxifying PAHs. When antioxidant defenses do not neutralize ROS, they could induce different types of damage at the cellular level, such as lipid and protein oxidation (Lushchak, 2011). In aquatic organisms, ROS exposition makes DNA susceptible to oxidative damage, and such damage can be measured as micronuclei, DNA breaks, and others using the comet assay test (Abele et al., 2011). Among the antioxidant defense system components, the most studied in marine organisms include catalase (CAT) and superoxide dismutase (SOD). Catalase is an enzyme that metabolizes hydrogen peroxide to molecular oxygen and water (Van der Oost et al., 2005). Considering the lagoons where cat was evaluated, Tampamachoco, Caŕmen, and Las Enramadas had a higher expression level. Studies reporting the up-regulation of cat in organisms exposed to hydrocarbons in the wild or under laboratory conditions include the polychaeta Perinereis aibuhitensis (Zhao et al., 2017) and Mus musculus (Wang et al., 2009). In general, the results reported in previous studies and the present work suggest the usefulness of cat as an environmental bioindicator.
As mentioned above, one of the results of the detoxification process of xenobiotics and PAHs is the generation of ROS and, therefore, oxidative stress conditions at the cellular level. Within this context, compounds such as thioredoxin reductase (TRR), glutathione reductase (GR), and glutathione peroxidase (GPx) are the primary antioxidant enzymes of the thioredoxin and glutathione systems found in plants, animals, and organisms (Meister and Anderson, 1983), so their fluctuation has been used as indicators of oxidative stress and cell damage in mussels (Gonzalez-Rey and Bebianno, 2013). In the case of thioredoxins, they are involved in different biological functions such as the elimination of reactive oxygen species, regeneration of proteins with oxidative damage, and control of apoptosis through the regulation of the transcription factor NF-kB (Fernando et al., 1992;Saitoh et al., 1998;Das and Das, 2000). In the present study, thioredoxin 3 (trr3) expression analysis showed that organisms exposed to hydrocarbons (21 days) and those from Las Enramadas showed the highest expression levels. In exposed organisms, direct and acute exposure compared to those from the natural environment would reflect the organism's rapid response to the pollutant's presence. On the other hand, it has been described that the presence of PAHs with a greater number of benzene rings is more difficult to eliminate and tends to produce stronger oxidative stress in the body (Motoyama et al., 2009). The above supports the results found working with C. virginica under laboratory conditions, where soft tissue bioaccumulation analyses indicated a higher proportion of PAHs of 4 to 6 benzene rings compared with those of 2 to 3 rings.
Genes from different function categories had been used as potential markers to assess the effect of hydrocarbon components on marine organisms (Albarano et al., 2021;Varea et al., 2021). In this sense, the most frequently selected genes as potential markers correspond to a member of the CYP family proteins, specially CYP1 and CYP2 (Ertl et al., 2016;Zacchi et al., 2019). In zebrafish, the increase of retene concentrations, from 0.205 to 50 μM, not only increased the number of differentially expressed genes but also increased the log 2 FC (FC = Fold-change) of cyp1a, proposing it as a marker of the negative effect of hydrocarbons on de development of zebrafish (Wilson et al., 2022). In the case of invertebrates, upregulation of cyp members in response to the presence of PAHs had been reported in Ruditapes philippinarum (Liu et al., 2014), Crassostrea virginica (López-Landavery et al., 2019; polychaeta Perinereis aibuhitensis (Zhao et al., 2022), Crassostrea gigas and Mytilus coruscus exposed to 1 and 5 μg/L during 15 days (Li et al., 2021). Based on the bioaccumulation results of this study and the organisms exposed to hydrocarbons at 200 mg/L under laboratory conditions (Loṕez-Landavery et al., 2019), our working hypothesis is supported in most of the selected genes, except for cat, fmo, and ggt1, where organisms from Tampamachoco had the highest levels of expression. Organisms exposed to 200 mg/L of hydrocarbons had levels of PAHs in soft tissues > 5 mg/g dry weight, and Tampamachoco had the highest level of PAHs in water and was among the top three of the lagoons with high content of PAHs in sediment. In the present study, genes such as cyp10, fmo, mdrd1, cat, and trr3 showed a high induction. From these, cyp10 or fmo gene expression, or a combination of both, would provide a complementary effective tool as biomarkers to study the effect of PAHs in marine invertebrates from the Gulf of Mexico.

CIGOM and other similar monitoring programs
This work is part of the Consorcio de Investigacioń del Golfo de Mexico (CIGOM), the first multidisciplinary research group focused in to study the physical, chemical, and biological characteristics of this important oceanic basin on the Mexican side. CIGOM was founded in 2015 and it specializes in multidisciplinary projects related to possible environmental impacts of the oil and gas industry on the marine ecosystems of the Gulf of Mexico (https://cigom.org). The biological component, includes the study of the effects of hydrocarbons, mainly PAHs, on invertebrates, fish, turtles, and plankton and recently incorporated the problem of sargassum macroalgae (Sargassum spp.). Considering invertebrate monitoring as indicators of pollution in the coastal lagoons of the Gulf of Mexico, a similar and older initiative was established initially by the US Environmental Protection Agency, continued by the National Oceanic and Atmospheric Administration (NOAA), called "The Mussel Watch Program (MWP)" (Kimbrough et al., 2008;Bricker et al., 2014), and then it extended as a model around the world (Ramu et al., 2007;Thebault and Rodriguez y Baena, 2007;Cossa and Tabard, 2020). CIGOM and MWP reflect the need to assess the impact of anthropogenic activities on marine and estuarine ecosystems. In comparative terms and considering only the monitoring of contaminants, the MWP has a broader scope and intensive sampling than CIGOM. MWP monitors nearly 300 sites from the Atlantic, Pacific, and Great Lakes coasts, and measures more than 140 chemical contaminants (Beiras, 2018). In contrast, and on the same base, CIGOM includes more species such as oysters (Crassostrea virginica), benthic and pelagic teleost, Chondrichthyes, cetaceans, turtles, and plankton (https://escenarios.cigom.org/en/). In a complete scope, although CIGOM includes more components than MWP such as oil circulation modeling, biological connectivity, and sceneries of vulnerability, among others, it has the challenge of staying for many years just like MWP.

Conclusion
The present study is one of the first where the AHs and PAHs levels present in organisms of different coastal lagoons of the Southern Gulf of Mexico, from Tamaulipas to Campeche, are evaluated together. It also integrated bioaccumulation data with molecular data to validate potential biomarker genes to assess the presence and impact of hydrocarbons in the different coastal lagoons on the east side of the Southern Gulf of Mexico, which is where the country's oil activity is concentrated. The first thing to highlight here is the perspective to continue studying the area of Las Enramadas, due to the expression patterns shown during the present study. Although bioaccumulation data have not been integrated with those from the water, sediment, and soft tissue of organisms from Las Enramadas, it suggests that the exploitation of gas and the proximity of gas pipelines could impact the southern region of the Laguna Madre lagoon complex due to its proximity to the hydrocarbon-producing basins of Burgos and Tampico-Misantla. The above generates the need to continue studying this zone to support the results found here. In the case of the presence of hydrocarbons, lagoons related to the existence of oil terminals with high maritime traffic or with fluvial influence, such as Tampamachoco, Mecoacań, Laguna de Teŕminos, and Caŕmen, presented the highest levels of hydrocarbons in any of the components evaluated, which confirms the partial impact of the presence of hydrocarbons on oyster response due to that others contaminants not evaluated here such as heavy metals influence the oyster response. Finally, and based on expression, although different biomarkers have been proposed using different species in other parts of the world, our results suggest using cyp10 or fmo as biomarkers for future assessment studies in coastal lagoons in the Gulf of Mexico.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/ genbank/, SRP183024; https://www.ncbi.nlm.nih.gov/genbank/, SRR18047696 to SRR18047710.