Sex-Specific Transcriptomic Differences in the Immune Cells of a Key Atlantic-Mediterranean Sea Urchin

The abundance of the black sea urchin, Arbacia lixula , has been increasing during the last decades likely related to global warming. This thermophilous species has a leading role in maintaining marine barrens in the Mediterranean with the consequent negative impact on coastal rocky ecosystems due to its grazing activity. In this study, we used transcriptomic data from coelomocytes (the cell effectors of the immune system) of females and males of this sea urchin to study potential differences in performance between sexes under laboratory conditions. Differential adaptations, responses to environmental stressors, and resistance against pathogens between sexes may lead to different outcomes in the ongoing expansion of this species in the Mediterranean Sea. Differential expression analyses demonstrated the existence of 120 transcripts, corresponding to 119 genes and two isoforms of the same gene, differentially expressed between coelomocytes of females and males, being 73 up-regulated in males and 47 up-regulated in females. The differential expression patterns were retrieved from a diversity of genes that play different roles related to the immune response due to their antibacterial activity, immune cell activation, cell to cell interaction, intracellular signaling, and detoxification functioning, among others. Our results point out a higher energetic demand of male coelomocytes due to a higher immune activity than females, whereas females have more efficient molecular systems to avoid oxidative stress caused by infections. In conclusion, our study provides evidence of sex-based differences in the expression of genes related to the immune and stress responses in coelomocytes of the sea urchin A. lixula .


INTRODUCTION
Sea urchins are among the most relevant keystone species in marine ecosystems, with a plethora of examples of the structuring effects of their grazing activity in sublittoral habitats, from tropical seas to kelp forests (e.g., Palacín et al., 1998;Edmunds and Carpenter, 2001;Hernández et al., 2008;Eklöf et al., 2008;Ling et al., 2009). Their overgrazing activity in coastal rocky assemblages causes regime shifts from productive macroalgal beds to impoverished urchin barrens, with anthropogenic Frontiers in Marine Science | www.frontiersin.org July 2022 | Volume 9 | Article 908387 stressors such as overfishing and climate change strengthening the resilience of urchin barrens (Hernández et al., 2008;Ling et al., 2009;Ling et al., 2015;Bernal-Ibáñez et al., 2021). Due to the invaluable ecosystem services provided by macroalgal and kelp beds (e.g., Smale et al., 2013;Rocha et al., 2015;Krause-Jensen and Duarte, 2016), it is of utmost importance to understand all details of the biology of the habitat-modifying echinoid species.
In the Mediterranean Sea, the two sea urchin species producing barrens are Paracentrotus lividus and Arbacia lixula (Palacín et al., 1998;Micheli et al., 2005), having A. lixula a leading role in maintaining the barren state once phase shift occurred . Arbacia lixula is a highly mobile species on open rocky surfaces that "scratches" on the substrate to feed on encrusting coralline algae and benthic invertebrates. This feeding behavior plays an important role in avoiding the macroalgae re-colonization, which finally promotes the persistence of large bare areas Agnetta et al., 2013). As a relatively new colonizer in expansion in the Mediterranean Sea (Wangensteen et al., 2012;Pérez-Portela et al., 2019), the thermophilous A. lixula has increased its abundance during the last century in the Mediterranean (Petit et al., 1950;Francour et al., 1994). Following this trend, and due to its superior physiological performance in warmer conditions (Privitera et al., 2011;Wangensteen et al., 2013a), global warming will likely favor A. lixula expansion and dominance in the whole Mediterranean basin , with its associated negative effects due to its role on the maintaining barrens areas, already evident in the Southern-Central Mediterranean .
Given this projection, an enhancing number of studies are using A. lixula as a model to investigate its genetic variation and adaptation patterns (Wangensteen et al., 2012;Pérez-Portela et al., 2019;Carreras et al., 2021), trophic ecology (Wangensteen et al., 2011), biological cycles (Wangensteen et al., 2013a) and responses to future climatic conditions and stressors (e.g. Pérez-Portela et al., 2020a). During the last decade, a wide array of techniques and approaches has been implemented, for instance, to measure its biological response to underwater highfrequency noise (Vazzana et al., 2020), neurotoxicity (Maisano et al., 2015), transcriptomic responses to temperature changes (Pérez-Portela et al., 2020a), and the effects of environmental changes on larval development, settlement and juvenile survival (Wangensteen et al., 2013b;Gianguzza et al., 2014). For these studies, understanding all factors that may affect the individual response to different treatments is essential, but in most of them, the variable sex is not considered despite its potential relevance in modulating the individual responses.
Genomic, transcriptomic, and proteomic studies of echinoderm coelomocytes have demonstrated the large complexity and flexibility of their innate immune system (e.g. Dheilly et al., 2013;Smith et al., 2018;Zhang et al., 2019). In this work we use transcriptomic data from coelomocytes of females and males of A. lixula to study potential differences in the immune system between sexes. The coelomocytes, the cell effectors of the immune system, consist of heterogeneous cell populations contained in the coelomic fluid, which are also involved in other physiological functions such as the stress response (e.g. Matranga et al., 2000;Pérez-Portela et al., 2016;Smith et al., 2018). In turn, these functions are considered fundamental elements of adaptation to new environments and face environmental shifts associated with climate change (e.g., Palumbi et al., 2014;Cunning et al., 2018;Liu et al., 2018). Previous studies in the sea urchin P. lividus, regarding the differential immune response according to sex, demonstrated that the immune system is more active in females than in males, with higher levels of cytotoxicity and hemolysis activity in females, likely related to a different cell composition within the coelomic fluid and a greater number of phagocytes and uncolored spherulocytes in females (Arizza et al., 2013). Additionally, other studies with marine invertebrates have also evidenced stronger immune responses in females than males (Nguyen et al., 2018, and references therein), although exceptions also exist (Jiang et al., 2017). Indeed, focusing on gene expression responses, more recent studies have identified patterns of sex-specific differentially expressed genes in different marine invertebrates (e.g., Li et al., 2019;Artal et al., 2020;Kadiene et al., 2020), which are also present in somatic tissues (e.g., Rotllant et al., 2017). Hence, it is crucial to understand the sex-specific expression of A. lixula coelomocytes, as different adaptation performance, responses to environmental stressors, and resistance against pathogens of males and females may lead to different outcomes in the ongoing expansion of this species in the Mediterranean Sea. Arbacia lixula lacks sexual dimorphism and has populations with a sex ratio close to 1:1, although slightly deviations in favor of males (Elakkermi et al., 2021) or females (Klaoudatos et al., 2022) have been detected in some populations.
Significant differences in transcriptomic patterns and gene expression have been described between gonads of males and females in A. lixula based on differential functioning of these organs, ovaries and testis (Pérez-Portela et al., 2016), but in the present study, we compare differences between sexes using the same somatic tissue type. Our initial hypothesis is that, following that observed in P. lividus, females of A. lixula would immunologically be more active than males, and therefore both sexes would display differential patterns of transcriptomic expression in coelomocytes.

MATERIAL AND METHODS
Raw reads and a reference transcriptome used in this study were obtained from Pérez-Portela et al, (2020a), available at Mendeley Data https://doi.org/10.17632/5673n552yj.1 and NCBI BioProject number PRJNA642021 (BioSamples SAMN15375803-SAMN15375808). We used transcriptomic information of three females and three males maintained under control conditions in laboratory.

Sampling and Experimental Design
Details of the sampling and experimental design are described in Pérez-Portela et al, (2020a). Briefly, adult specimens of A. lixula (~48-52 mm in diameter) were collected from a subtidal population of Blanes, in Cala Sant Francesc (41°40′22.47″N, 2°48′10.81″E, North-Western Mediterranean), to minimize potential variation due to genetic divergence among populations.
Once in the laboratory, sea urchins were maintained in an aquarium with airflow and flow-through running seawater at 13˚C. For coelomocytes collection and RNA extraction, sea urchins were pulled out from the aquaria, quickly dissected under RNAase-free conditions, and coelomocyte fluid was collected. Sex determination in sea urchins was done by visual examination of the gonads, and posterior histological analyses. Coelomocytes of each specimen (a total of six: three females and three males) were gathered from the coelomic fluid and RNA was extracted following a protocol optimized elsewhere (Pérez-Portela and Riesgo, 2013;Pérez-Portela et al., 2016). High-quality mRNA was isolated and cDNA libraries were constructed with an average insert size of 300 pb. Libraries were multiplexed and sequenced in an Illumina HiSeq2000 Sequencer (101 base paired-end).
The software FASTQC v. 0.10.0 (www.bioinformatics. babraham.ac.uk) was used to visualize and measure the quality of the raw reads generated. Adapters and bases with low quality (phred scores <33) were trimmed off, and a length filter was applied retaining sequences of >25 bases using TrimGalore v. 0.2.6 (www.bioinformatics.babraham.ac.uk).

Differential Expression Analyses and Annotation
Reads from each sample were aligned against a reference transcriptome for coelomocytes of A. lixula (reference_ assembly_CT+T22.fasta, available at https://data.mendeley. com/datasets/5673n552yj/1, Pérez-Portela et al., 2020b). This reference transcriptome was de novo assembled by Pérez-Portela and coauthors (2020b) from 12 different coelomocyte samples of A. lixula, and it was the most complete one among those available for coelomocytes in this species (Pérez-Portela et al. 2020b). Only transcripts over 500 pb from the reference transcriptome were used for mapping our reads, and only paired reads after trimming were mapped using Bowtie2 v. 2.2.1 (Langmead and Salzberg, 2012) as implemented in Trinity (Grabherr et al., 2011). RSEM v. 1.2.11 (Li and Dewey 2011) was then run to generate a table with read counts, and unmapped reads were discarded. We retained information on differential expression of all isoforms detected for a given gene (or component) because different isoforms of the same gene may have different functions.
Differential expression (DE) analyses between females and males (three replicates each) were performed with the package DESeq2 (Love et al., 2014) in R v 4.0 (www.CRAN.R-project.org; R Development Core Team 2008). Before analyses of differential gene expression, read counts were internally normalized by the "DEseq()" function, and then a negative binomial model was fit to accurately estimate differential expression. The significance value for multiple comparisons was adjusted to 0.01 for False Discovery Rate (FDR) with the function "padj" (Benjamini-Hochberg adjustment) as implemented in DESeq2. Visualization of the significant outcomes of isoforms differentially expressed (up-and down-regulated) between females and males was obtained with a heatmap performed with the gplots package of R (Warnes et al., 2009).
Using the Gene Ontology (GO) annotation results for the reference transcriptome, we obtained both the gene annotations and the GO terms associated with the differentially expressed isoforms. GO terms and their associated p-value (adjusted) were then summarized in the REVIGO webserver (Supek et al., 2011), and graphically represented with the treemap R package (Tennekes and Ellis, 2017). The size of the rectangles was adjusted to reflect the adjusted p-value to visually detect the most significant terms. Differentially expressed isoforms without blast hit, unknown function and/or without annotation were isolated, and the InterProScan 5 software (Jones et al., 2014), that predicts protein family membership and the presence of functional domains and sites, was run at the SUPERFAMILY level (De Lima Morais et al., 2011). The InterProScan was run as implemented in the Blast2GO software v 3.31.1 (Conesa et al., 2005).

Differential Expression Between Sexes
A total of 283,880,164 reads from all samples, and 66,985 transcripts from the reference transcriptome, including all genes and isoforms, were used for the differential gene expression analysis in this study. All samples had over 23,000,000 reads (ranging from 23,074,766 to 84,053,170) after filtering and trimming. A total of 120 transcripts, corresponding to 119 genes FIGURE 1 | Heatmap based on the differentially expressed transcripts (DE) between females and males of A. lixula. The clustering on the top of the heatmap represents the similarity among the samples analyzed. F1-F3 and M1-M3 represent the three female and three male replicates, respectively. The term "Value" represents the gene expression value. and two isoforms of the same gene, were differentially expressed (p-adjusted value < 0.01) between coelomocytes of females and males, being 73 up-regulated in males and 47 up-regulated in females. The heatmap that describes the hierarchical clustering of all transcripts related to their expression differences between sexes, including three replicates for each sex, is displayed in Figure 1. The heatmap shows a clear clustering of all replicates per sex.

Functional Annotation
Out of the 120 differentially expressed transcripts, 40 (~33.3%) had blast hit against known genes. The InterProScan analyses Sex-Specific Transcriptomic Differences in Echinoids could only predict protein domains of three uncharacterized transcripts. A complete list of differentially expressed annotated transcripts is presented in Table 1, including transcript identification code (id), log2fold change, p-adjusted value (padj), and gene annotation, and GO-terms. Among the transcripts differentially expressed between females and males, we can find some codifying genes described in other studies as directly related to the immune response, such as: the zinc finger ccch domain-containing protein 14-like (log2fold change= -8.53), the ras-related c3 botulinum toxin substrate 1-like (log2fold change= -8.17), the hepatic lectin (log2fold change= -2.68), and the bactericidal permeability-increasing (log2fold change= -3.3), all of them up-regulated in males. For females, we found up-regulated the echinonectin (log2fold change= 9.07) and ferricchelate reductase 1 like (log2fold change= 6.78).
Other genes codifying for potential proteins involved in the immune activity, due to its role in cell to cell interaction and intracellular signaling, were also DE between sexes, including the activating signal cointegrator 1 complex subunit 3-like (log2fold change= -7.35) and the low-density lipoprotein receptorrelated protein 2-like isoform 2 (log2fold change= -1.00) both up-regulated in males. Transcripts related to the transmembrane transport functions as two isoforms of the major facilitator superfamily domain-containing protein 7-a-like (isoforms 1 and 2, please see c275529_g1_i1 and c275529_g1_i2 in Table 1) were also identified as differentially expressed between sexes (log2fold FIGURE 2 | Gene ontology (GO) treemaps for DE genes annotated in A. lixula. Those up-regulated in males and females are represented separately for Biological processes, Cellular components, and Molecular functions. The size of the rectangles reflects the adjusted p-values associated with the differentially regulated categories. Frontiers in Marine Science | www.frontiersin.org July 2022 | Volume 9 | Article 908387 change= 8.19 and -8.05, respectively). Several genes codifying for integral components of the cellular membranes were also DE, including the transmembrane protein 85-like and alphap integrin precursor, both up-regulated in males (log2fold change= -8.22 and -9.03, respectively), and the ADP, ATP carrier protein-like gene up-regulated in females (log2fold change= 7.51), among others. As expected, the GO terms associated with these DE genes were also different between males and females. Out of the 35 GO categories recovered, 22 were derived from genes up-regulated in males, and 13 from genes up-regulated in females. Figure 2 depicts the up-regulated GO categories associated with the DE genes. The GO term categories found can give us a general idea of the biological processes, cellular components, and molecular functions differentially represented between sexes. Among the most significant GO terms up-and/or down-represented in females and males, we find some terms related to intracellular signaling and transport such as "transmembrane transport", "protein binding", "integral component of membranes", "membrane"; those related with cell to cell interaction such as "cell adhesion"; with the stress-responses such as "ferricchelate reductase activity" and "metal ion binding"; with energy metabolism as "peptide cross-linking", "NADH dehydrogenase complex" and "protein-glutamine gamma-glutamyl transferase activity"; and other general categories such as the biological processes of "RNA-dependent DNA biosynthetic process", and their associated cellular components and molecular functions.

DISCUSSION
Our transcriptomic analysis from coelomocytes of the sea urchin A. lixula demonstrated the existence of significant differences in the expression of some genes, and isoforms, between females and males. These patterns of differential expression were retrieved from a diversity of genes that seem to play different roles during the immune and stress responses in sea urchins. A study in P. lividus identified a different composition of the coelomocyte populations between sexes, with females being immunologically more active (Arizza et al., 2013). In our study, we detect differences in transcriptional levels of expression between females and males of A. lixula, but a higher number of up-regulated genes are associated with the innate immune response in males. Although the relationship between levels of transcription and abundance of proteins, and the functional consequences of differential expression of some particular genes are, on some occasions, still unclear (Feder and Walser, 2005), some studies have shown the correlation between expression levels and protein abundance under experimental conditions (Maier et al., 2009). Therefore, changes in gene expression can be used as an approximation of the physiological performance of a tissue and/or organism under laboratory-controlled conditions (Feder and Walser, 2005;Buckley et al., 2006;Schoville et al., 2012), and so, in this work, the number of immune-related genes up-regulated can be considered as a proxy to infer the immune activity. Hence, A. lixula results suggest that, contrary to that described in P. livudus, males might be immunologically more competent than females, or at least, both sexes use different immune pathways. Immune-competence analyses with sea cucumber Apostichopus japonicus, also evidenced a more powerful immune system in males over females (Jiang et al., 2017), which is in line with our results in A. lixula.
Although many genes are potentially involved in the innate immune response in sea urchins (e.g. Dheilly et al., 2013;Zhang et al., 2019), in this study, we only discuss some of the most relevant ones and/or those playing a well-known role in the immunity, cell adhesion, cell signaling, stress, and detoxification. Regarding the genes directly involved in the immune response and cell adhesion, at least, four of them were up-regulated in males: the bactericidal permeability-increasing protein, the zinc finger ccch domain-containing protein 14-like, ras-related c3 botulinum toxin substrate 1-like, and a hepatic lectin. The bactericidal permeabilityincreasing gene encodes for a lipopolysaccharide-binding protein with potent cytotoxicity against Gram-negative bacteria, which has been isolated in a diversity of echinoderm tissues (Shao et al., 2015), but a DE between sexes has never been described before. The zinc finger ccch domain-containing protein 14-like has a crucial role in the regulation of cytokines production, immune cell activation, and immune homeostasis (Fu and Blackshear, 2017), and it is associated with the molecular function GO term "protein binding". The ras-related c3 botulinum toxin substrate 1-like gene showed significant differences in expression between males and females and its foldchange was among the lowest in the whole dataset. This gene codifies for a plasma membrane protein, involved in cell signaling, a member of the Rho family of small GTPases that regulates processes of inflammation and apoptosis (Hirohama et al., 2020), and stimulates the NADPH oxidase activity and phagocytosis in macrophages (Dheilly et al., 2013). The hepatic lectin gene codifies for a molecule secreted by the coelomocytes to the coelomic fluid as a part of the immune humoral system. Lectins are a large group of proteins and glycoproteins, soluble or integral components of membranes, that act as mediators of the immune system in invertebrates, including echinoderms, due to their function in the cell adhesion, and cell to cell and cell-extracellular matrix interactions (Dheilly et al., 2013). They can specifically recognize and bind to the microbial surface and lead to phagocytosis (Dheilly et al., 2013;Zhang et al., 2019).
In females, two genes potentially related to the immune response were detected as "up regulated": the echinonectin and the ferric-chelate reductase 1-like. The echinonectin gene encodes for an echinoderm adhesion protein, a lectin, that was initially described during the embrionary development, and that also plays an important role in skeleton formation (Mann et al., 2008). Nevertheless, more recent studies might suggest its role in the defense against pathogens as other lectins (Dheilly et al., 2013) due to its important role in cell adhesion. The ferricchelate reductase 1-like gene encodes a molecule linked to the membrane function and iron regulation that reduces Fe +3 to Fe +2 , and is associated to the molecular functioning GO terms "ferric-chelate reductase activity" and "metal ion binding". In sea urchins, iron plays a critical role in the immune defense against pathogens stimulating cell proliferation, and nitric oxide and reactive oxygen elements (Dheilly et al., 2013). In vertebrates, iron sequestration avoids it to be used by microorganisms (Brock and Mulero, 2000;Ganz, 2009), a function that iron-binding proteins could be also playing in sea urchins (Dheilly et al., 2013;Zhang et al., 2019). However, the specific involvement of the ferric-chelate reductase 1 is not mentioned in previous studies on echinoderms. On the other hand, some studies with sea stars have detected expression of this gene in females and signals of selection on it, promoting adaption to reduce oxidative stress (Guerra et al., 2021). This gene is particularly interesting since in other marine invertebrates, such as the mollusk Ruditapes philippinarum, females showed more efficient pathways against oxidative stress caused by infections than males (Matozzo and Marin, 2010), and so our results from A. lixula could also support this idea.
We also highlight the importance of gene isoforms in the coelomocytes' functioning. We detected the existence of two isoforms of the major facilitator superfamily domain-containing protein 7-a-like coding gene, involved in the biological processes of "transmembrane transport" and "integral component of membranes", with an opposite pattern of expression between sexes and almost equivalent foldchange. This gene expression pattern suggests the sex-specificity of these two isoforms. Hence, it does not seem rare to detect the existence of isoforms with opposite patterns of differential expression between sexes, as in other marine invertebrates, such as abalones, the presence of gene isoforms was also sex-dependent (Kim et al., 2017).

Previous transcriptomic analyses in
In general, and independently of the particular genes upor down-regulated per sex, different GO terms inferred from A. lixula were also identified in the immune response against bacteria artificially inoculated in Strongylocentrotus intermedius (Zhang et al., 2019), and coelomocytes' proteomic analyses of Strongylocentrotus purpuratus (Dheilly et al., 2013), including the "transmembrane transport", "membrane", "integral component of membrane", "cell adhesion", "protein binding" and "metal ion binding", among others. These similar results with other sea urchin species support the relevant role of the cell interactions throughout the cellular membranes for the immune response in both sexes of A. lixula. Another interesting molecular pathway potentially linked to immune functions is the lipid metabolism due to its relevance in the development of the immunological response (Miller et al., 2003). Within it, we detected the low-density lipoprotein receptorrelated protein 2-like gene differentially expressed between sexes of A. lixula. These cell receptors are involved in capturing lowdensity lipoproteins for intra-and extra-cellular signaling and are also identified in coelomocytes of S. purpuratus (Dheilly et al., 2013). Finally, we identified the up-regulation of the transglutaminase-like protein and nadh dehydrogenase codifying genes in males, both associated with the energy metabolism and the GO terms "peptide-cross linking", "NADH dehydrogenase complex" and "protein-glutamine gamma-glutamyl-transferase activity". This finding, as in other studies on echinoderms (e.g. Dheilly et al., 2013;Zhang et al., 2019), points out a higher energetic demand in coelomocytes of males due to a higher immune activity.
In conclusion, our study provides evidence of sex-based differences in the expression of some genes directly and indirectly related to the immune and stress responses. The potential higher immune activity in males than in females may lead to an advantage against pathogens but under the cost of a higher energetic demand. On the other hand, females seem to be better adapted to the oxidative stress associated with infections. Both sexes did not display significant differences in transcriptomic responses to thermal stress under controlled conditions in an aquarium (Pérez-Portela et al., 2020a) but did it under control conditions. Whether the absence of transcriptomic divergence between sexes under stressful thermal conditions is due to the low replication of one of the sexes (three replicates per sex) or the dilution of the potential of those differences under stress, is a question that we cannot answer with the current data. Future studies should focus on the differential immune response against pathogens between sexes under variable environmental conditions, to provide a complete picture of the relevance of these sex-associated differences in transcriptional patterns between females and males of A. lixula.

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: NCBI [accession: PRJNA642021].