Octopus maya white body show sex-specific transcriptomic profiles during the reproductive phase, with high differentiation in signaling pathways

White bodies (WB), multilobulated soft tissue that wraps the optic tracts and optic lobes, have been considered the hematopoietic organ of the cephalopods. Its glandular appearance and its lobular morphology suggest that different parts of the WB may perform different functions, but a detailed functional analysis of the octopus WB is lacking. The aim of this study is to describe the transcriptomic profile of WB to better understand its functions, with emphasis on the difference between sexes during reproductive events. Then, validation via qPCR was performed using different tissues to find out tissue-specific transcripts. High differentiation in signaling pathways was observed in the comparison of female and male transcriptomic profiles. For instance, the expression of genes involved in the androgen receptor-signaling pathway were detected only in males, whereas estrogen receptor showed higher expression in females. Highly expressed genes in males enriched oxidation-reduction and apoptotic processes, which are related to the immune response. On the other hand, expression of genes involved in replicative senescence and the response to cortisol were only detected in females. Moreover, the transcripts with higher expression in females enriched a wide variety of signaling pathways mediated by molecules like neuropeptides, integrins, MAPKs and receptors like TNF and Toll-like. In addition, these putative neuropeptide transcripts, showed higher expression in females’ WB and were not detected in other analyzed tissues. These results suggest that the differentiation in signaling pathways in white bodies of O. maya influences the physiological dimorphism between females and males during the reproductive phase.


Introduction
The white body (WB) was first described by Cuvier [1] as a "corps glanduleux" surrounding the optic lobes and optic traits with the only function to protect these structures during muscle contraction. Later the WB was considered as site for the formation of leucocytes [2,3]. Cazal and Bogoraze [4] proposed a second function for the WB: "fonction nèphrocytaire", and both functions were supported by the work of Bolognari [5]. Then Harrison and Martin [6] demonstrated an important role of WB in urine formation; however Young [7], retained that its function was still unclear.
In cephalopods, it is well known that haemocytes are originated within the WB, in fact, several hematopoiesis genes were found to be expressed in the WB of squid Euprymna tasmanica [8]. In addition, transcripts associated with immune-related signal transduction pathways were found, as well as other genes of the immune response previously identified in E. scolopes [9]. Nevertheless, its general similitude with the mammal lymphoid tissue, its lobular morphology, and its glandular appearance suggest that different WB regions may perform different functions [10,11]; but in octopuses, a detailed functional analysis of the WB is lacking. This organ links the nervous and circulatory systems and may play a role in the cooperation between neuroendocrine and immune responses to environmental stimuli [12]. However, the understanding of how these two systems cooperate, remains unclear [12]. Such interaction may be mediated via secretion/reception of hormones, neuropeptides or signaling peptides. The neuroendocrine crosstalk among different lobes of the nervous system, which regulates physiology and reproduction of O. vulgaris, is well documented [12][13][14][15][16]. However, if the WB takes a place in this crosstalk is poorly understood. Herein, to better understand the WB functions, this study was directed to describe its gene expression, with emphasis in sex-related and reproductive-stage related differences via RNA-Seq. This technology has been successfully implemented for gene discovery and for estimation of gene expression levels in different cephalopod tissues [8,[17][18][19][20][21]. In addition, different tissues were compared in terms of gene expression to identify sex-specific and tissue-specific transcripts.

Ethics statement
In this study, octopuses were anesthetized with ethanol 3% in sea water at experimental temperatures [22,23] to enable humane killing [24] in consideration of ethical protocols [25], and the animals' welfare during manipulations [26,27]. Our protocols were approved by the experimental Animal Ethics Committee of the Faculty of Chemistry at Universidad Nacional Autónoma de México (Permit number: Oficio/FQ/CICUAL/099/15). We encouraged the effort to minimize animals stress and the killing of the minimum necessary number of animals for this study.

Acclimation and experimental design
Female and male octopuses with a body mass that ranged between 400-600 g were captured off the coast of Sisal Yucatán, México, by artisanal fishing fleet and then transported to the Experimental Cephalopod Production Unit at the Multidisciplinary Unit for Teaching and Research (UMDI-UNAM), Sisal, Yucatan, Mexico. Octopuses were acclimated for 10 d with 1:1 sex ratio, in 6 m diameter outdoor ponds provided with aerated natural seawater (26 ± 1˚C). The ponds were covered with black mesh reducing direct sunlight to 70% and connected to seawater recirculation systems coupled to protein skimmers and 50 μm bag filters. PVC 50 mm diameter open tubes were offered as refuges in proportion 2:1 per animal.
Octopuses were fed individually twice a day with a paste made with squid and crab meat at ratio of 8% of its body weight [28]. Food not ingested, and feces were removed daily. During acclimation octopuses paired freely; then, a group of ten males were sampled before copulation and another group of ten were sampled after copulation. On the other hand, twenty fertilized females were individually reared in 80L tanks in a recirculating aquaculture water system. Each tank was provided with a fiberglass box that serves as refuge for the female and for spawn settling. System temperature was maintained at 24˚C, until all females spawned. These conditions were selected because 24˚C is the preferred temperature of this species and has been recommended as the best condition for spawning [29,30]. Water was heated by using a 1200W titanium immersion heater connected to a digital temperature sensor, both placed at system reservoir; and was cooled by using air conditioning according to the temperature required.

White body sampling
Octopus sensibility was minimized before manipulation, by means of an anesthetic procedure that consisted in maintaining animals in 3% alcohol-sea water solution for up to 4 minutes [23]. White body samples were obtained from ten males before copulation (MPRE) and from ten males after copulation (MPOS). In the case of females, ten WB samples were obtained before spawning (FPRE) and another ten samples after spawning (FPOS). All WB samples were taken from the region adjacent to the optic lobe. Additional samples from different tissues (including hearts and gonads) were obtained from the same experimental individuals, for further analyses. All tissue samples were preserved in RNAlater solution until RNA extraction.

RNA sequencing
Total RNA was extracted from 30 mg of WB tissue using the RNeasy extraction kit (Qiagen) following manufacturer instructions. RNA was quantified with NanoDrop 2000 spectrophotometer (Thermo Scientific) and quality was assessed using the Bioanalyzer Instrument 2100 (Agilent Technologies). Samples that presented an RNA integrity number (RIN) equal or higher than 7 were used for sequencing. Using equal amounts of RNA from the three individuals with higher RIN in each experimental condition, paired-end cDNA libraries were prepared using the TruSeq DNA Sample Preparation Kit v2 (Illumina), following the manufacturers protocol. Subsequently, cluster generation and DNA sequencing were performed in MiSeq sequencing system (Illumina) to obtain reads of 250 bp long.

Differential expression analysis
The transcriptome assembled with all female and male reads was used as reference to estimate the abundance of the transcripts for each library: FPRE, FPOS, MPRE and MPOS. The reads from each library were aligned back to the reference transcriptome by using Bowtie2 v2. 3

Functional annotation of transcripts
The reference transcriptome (assembled from the reads of both sexes) was analyzed with BLASTx [37] to find homologs within UniProt Release 2017_12 database with an e-value < E-05 filter. The Gene Ontology (GO) annotations were obtained based on the UniProt IDs using Blast2GO v4.1 with an e-value filter of 1E-08 [38]. The GO annotations for each transcript were analyzed to identify the best-represented biological processes detected in the reference transcriptome, based on the number of sequences included in each GO category. On the other hand, a Fisher exact test (FDR < 0.05) was implemented to identify the enriched biological processes in each library, using the annotated DE transcripts as test-set and the reference transcriptome as background in Blast2GO v4.1. In the same way, an additional enrichment of metabolic pathways from Kyoto Encyclopedia of Genes and Genomes (KEGG) database [39] was performed, using DAVID 6.8 [40] with the annotated DE transcripts as test-set and the assembled reference transcriptome as background. From the enriched categories, transcripts were selected for further analysis to confirm the annotation results. These transcripts were analyzed with TransDecoder v5.5.0 [41] to identify their longest open reading frames (ORF), and to predict their coding sequences and encoded peptides. Then, BLASTp searches against the nonredundant (nr) protein database were performed online (https://blast.ncbi.nlm.nih.gov/Blast. cgi) for the predicted encoded peptides.

Phylogenetic analysis of relevant DE transcripts
A phylogenetic analysis was performed to corroborate the BLASTx results of DE transcripts involved in relevant biological processes. A nucleotide sequence alignment was built, including all the mRNA sequences from mollusks available in the Genbank nucleotide database, corresponding to the putative gene analyzed. Nucleotide and codon alignments were carried out with the ClustalW 2.0 algorithm [42]; then the alignments were analyzed to find the nucleotide substitution model that best described the data, using the Maximum Likelihood method (ML). The model with lowest Bayesian Information Criterion score was considered the best [43]. The phylogenetic relationship among the sequences was inferred by using the ML method based on the General Time Reversible model [43]. Initial tree for the heuristic search was obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood approach, and then selecting the topology with superior log likelihood value. The bootstrap consensus tree inferred from 200 replicates was taken to represent the evolutionary history of the sequences analyzed [44]. Evolutionary analyses were conducted in MEGA7 [45].

Validation of DE transcripts via qPCR
A quantitative real-time PCR analysis was performed to validate the sex-related and tissue specific expression of transcripts, using additional tissues (systemic heart and testis). RNA samples were treated with RQ1 RNase-free DNase (Promega) before cDNA synthesis. The cDNA synthesis was carried out with the Improm II Reverse Transcription System (Promega) following manufacturer's instructions starting with 1μg of RNA from each sample. Primers for qPCR were designed with Primer3 [46] based on the selected transcript sequences (Table 1). To calculate the primer amplification efficiency, a standard curve was built including five stepwise dilutions of the cDNA samples with a constant dilution factor of 1:5, using nuclease-free water. The amplification reaction included 7μL of B-R SYBR Green Super Mix (Quanta Biosciences) and 3 μL of each cDNA dilution. Cq values were obtained in a C1000 Touch Thermal Cycler including the CFX96 Real Time System (Bio-Rad) with an amplification program consisting in: 2 min at 95˚C, 35 cycles of 30s at 95˚C, 30s at 57˚C, 30s at 72˚C and plate read, followed by 5 min at 72˚C, 10s at 95˚C, 30s at 65˚C, 60 cycles of 5s at 65˚C + 0.5˚C/cycle and plate read. The amplification efficiency of the primers (E) was calculated as E = [10(-1/k)-1], with k representing the slope of the standard curve [47]. For relative expression analysis, all reference and target transcripts were amplified by qPCR reactions within 96-well plastic plates (including positive and non-template controls) using the same instrument, reagents and amplification program as performed for amplification efficiency calculation. Nine biological replicates per sex and tissue were analyzed individually and by triplicate. The relative expression of target transcripts in each group, were estimated using the ΔΔ Cq method in Qbase+ 3.0 software [48] using as reference the heterogeneous nuclear ribonucleoprotein D like (HNRNPD) and the V-type proton ATPase subunit D (VATD) putative genes. Reference genes were validated previously using geNorm [49] and NormFinder [50]. For statistical analysis of relative gene expression values from RT-qPCR, one-way ANOVA and Tukey's multiple comparisons test were performed (P < 0.05). Moreover, Kruskal Wallis and Dunn's post hoc tests (P < 0.05) were used when assumptions of normality (Shapiro-Wilk) were not satisfied. To test if the relative expression values obtained by RNA-Seq and qPCR were consistent, a Spearman correlation was performed (P < 0.05). Statistical analyses were performed using Sta-tistica7 software [51].

Transcriptome assembly
For the reference transcriptome 75,265 unigenes and 90,435 isoforms were reconstructed, showing an average contig length of 622 bases, an N50 of 875 bases and a total of 56,280,138 bases assembled.

Differential expression analysis
A total of 3,522 WB transcripts showed significant differential expression between females and males (FDR < 0.01). From these, 2,192 showed higher expression in females and 1,330 in males. By contrast, differential gene expression was much lower comparing before and after the reproductive events. Between pre and post-spawning females, only 60 transcripts showed significant differential expression; and in males, 140 transcripts showed significant differences between pre and post copula condition (Fig 1).

Functional annotation and enrichment analysis
Using BLASTx searches, 17,656 transcripts of the reference transcriptome showed significant hits against the UniProt database, corresponding to 10 339 putative peptides. Then, after separating transcripts by sex, 13,399 were annotated in females, corresponding to 9,175 peptides; on the other hand, 13,265 male transcripts corresponding to 9,429 peptides showed significant hits in this database. From these, 5,926 hits were shared, 3,249 were exclusive for females and 3,503 for males. For the reference transcriptome the best represented gene ontologies in terms of biological processes (BP) are shown in Fig 2 where the response to stress, cell communication, signaling, and developmental processes appeared among the most specific categories. Moreover, in the assembled transcriptome we found at least 74 unigenes from the immune system category (Table 2), as well as 25 additional unigenes coding for different proteins considered as hematopoietic fingerprints [52] (Table 3).
According to the enrichment analysis based on KEGG database, there was high differentiation in signaling pathways between female and male white bodies (Table 4). In the females there was significant enrichment for a wide variety of signaling systems such as Toll-like receptor, PI3K-Atk, TNF, MAPK, mTOR, Rap1, oxytocin, JAK-STAT, VEGF, Ras and insulin among the most important (S1-S13 Figs). By contrast, in males there was enrichment of pathways involving translation initiation factors (eIFs) and apoptotic pathways related to an increment in reactive oxygen species (ROS) (S14-S17 Figs).
On the other hand, there were also significant GO enrichments based on the transcripts with higher expression in each sex (Table 5). In the females, unigenes involved in the neuropeptide signaling, integrin-mediated signaling, and signaling regulation pathways were highly expressed. Among the transcripts included in the category of neuropeptide signaling, we found five sequences putatively coding for FMRF-amide neuropeptide. These sequences were further analyzed to obtain the predicted coding sequence and the encoded peptide. The longest predicted peptide consisted in 143 amino acids with YIPF repeats each 12 amino acids, but not FMRF repeats were detected. This peptide obtained significant BLASTp hits (1E-10) with the Enterin neuropeptide of the clam Mizuhopecten yessoensis (NCBI accession number OWF47724.1). The aligned amino acid sequences presented Y+PF matches each 12 positions with no gaps, and mismatches in the non-repetitive region. This confirms the similarity of these transcripts with other mollusk neuropeptides.
In addition, GO terms like stem cell division, response to starvation, response to cortisol and replicative senescence were enriched in females. By contrast, the males showed significant enrichment for GOs like anatomical structure development, microtubule-based process, cilium morphogenesis, protein modification process and regulation of apoptosis among the bestrepresented categories. The unigenes that best represented the enriched GO categories are shown in Table 6.

Phylogenetic analysis of relevant DE transcripts
Phylogenetic analysis for the putative FMRF-amide sequences was performed, including five assembled transcripts. The analysis involved 90 nucleotide sequences and a total of 255 positions in the final dataset. Different topologies were obtained for the codon-based analysis and the nucleotide-based analysis. In the codon-based analysis, the putative FMRF-amide sequences were clustered with Aplysia californica PRQFVamide precursor protein mRNA (AY231295.1). However, in the nucleotide-based analysis they were clustered with Lymnaea stagnalis pedal peptide preprohormone mRNA (AY297820.1) and Loligo pealei FMRF-amide precursor, mRNA (FJ205479.1). Subtrees for both analyses, including the assembled transcripts and their evolutionarily closer sequences are shown in

Quantitative qPCR analysis
We selected transcripts based on three criteria a) their highly significant differential expression in RNA-seq analysis, like the transcript putatively coding for FMRF-amide neuropeptide, b) their relevance in signaling/hormonal process [53] like estradiol 17-beta-dehydrogenase 8 (HSD17B8), estrogen receptor (ESR1), corticotropin-releasing factor receptor 2 (CRHR2), c) unexpected transcripts that may represent novel functions of the WB like ATP-dependent RNA helicase DDX4 (VASA homolog), Piwi-like protein 1 (PIWIL1) related to germ cell development [54][55][56][57][58], as well as zonadhesin (ZAN) and C-Jun-amino-terminal kinase-interacting protein 4 (SPAG9) transcripts involved in species-specific fertilization [59,60]. All selected transcripts were correctly amplified in real-time qPCR reactions. The expression values of CRHR2, HSD17B and VASA genes failed the Shapiro-Wilk normality test and were analyzed with non-parametric statistics. All transcripts showed significant differential expression among the analyzed groups, considering their parametric or non-parametric distribution (P < 0.01). Their estimated relative expression was represented in histograms with mean values and 95% confidence intervals in each sex-tissue group. Tissue-specific and sex-specific gene expression were detected: the putative FMRF-amide like neuropeptide was detected only in the WB, with higher expression in females than in males (Fig 4); whereas the corticotropinreleasing factor receptor 2, showed higher expression in females' WB and heart (Fig 5). On the other hand, despite the expression of the estrogen receptor in the WB, the highest expression was observed in the testis tissue; nevertheless, its expression was significantly higher in females' WB compared to males' WB (Fig 6). Similarly, HSD17B8 showed relative low expression in WB and the highest expression in testis (Fig 7). In addition, the transcripts related to germ cell development and spermatogenesis detected in the WB like the ATP-dependent RNA helicase DDX4 (VASA) and Piwi-like protein 1 (PIWIL1) were not tissue-specific and showed higher Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase expression in testis (Figs 8 and 9). However, ZAN that is involved in species-specific fertilization, was highly expressed in WB, showing similar expression levels to those in testis (Fig 10). Finally, SPAG9, which is also involved in fertilization showed the highest expression in the females' WB (Fig 11). The Spearman correlation for the expression values obtained by RNA--Seq and qPCR was significant (P < 0.05), suggesting that both methods reported same tendencies in the expression levels across all the analyzed genes and samples. Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase The test-set consisted in the transcripts with significant differential expression. The transcripts from the assembled reference transcriptome were used as reference-set. https://doi.org/10.1371/journal.pone.0216982.t005 Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase

Discussion
One of the most relevant findings in this study was the different transcription pattern between male and female white bodies of O. maya, during the reproductive phase. The unigenes involved in this difference, fall mainly in the category of signaling pathways, and showed higher expression in females. Numerous signaling cascades were highly active and occurring simultaneously in the females. The enriched KEGG pathways suggest that bacterial and viral infections were taking place in females (Table 4), triggering multiple signaling systems. At the same time, the neuropeptide signaling pathway was also enriched in females by the high expression of at least 5 different transcripts coding for a putative FMRF-amide like neuropeptide. FMRF-amide related peptides (FaRPs) constitute an evolutionary conserved and diverse group of neuropeptides in the central nervous system (CNS) of many metazoans [61]. Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase Among the functions of FMRF-amide and FaRPs in mollusks are the following: the modulation of sensory organs, reproduction, motility, osmoregulation, feeding and neurogenesis [62,63]. In O. vulgaris, Di Cosmo and Di Cristo [64], and Di Cristo et al. [65], demonstrated the presence of FMRF-amide both in the CNS and peripheral nervous system (PNS). In the CNS this neuropeptide is involved in the inhibition of the secretory activity of the optic gland, which in turn controls the gonad maturation [64]. FMRF-amide has been detected in several lobes such as optic, subpedunculate and olfactory lobes of O. vulgaris CNS [66]. Significant hits with this neuropeptide using the UniProt database were obtained in this study. However, when the transcripts sequences were phylogenetically analyzed including all reported molluscan neuropeptides within the NCBI nucleotide database, they were clustered in different ways, depending on the sequence alignment method. In the codon-based analysis the putative neuropeptide transcripts were clustered with Aplysia californica PRQFVamide precursor protein mRNA, which is involved in the regulation of the feeding system [67]. By contrast, in the nucleotide-based analysis, the same transcripts were clustered with Lymnaea stagnalis pedal peptide preprohormone mRNA, which was highly expressed during parasitic infections in the great pond snail [68]; this branch also included the Loligo pealei FMRF-amide precursor, mRNA (FJ205479.1). Despite certain similarity at the transcript level with the FMRF-amide sequence, our query sequences do not contain FMRF repeats, instead they present YIPF repeats each 12 aa according to the longest ORF predicted. This 12-aa distance between repeats is also observed in the Enterin neuropeptide of the clam M. yessoensis, which may play a role in non-feeding behaviors [69]. This suggests that these tissue-specific transcripts with higher expression in females, show a neuropeptide-like structure, especially due to the constant and equidistant YIPF repeats. This is the first report of putative neuropeptide-like transcripts detected in the octopus' WB. In this regard, further analyses are needed to clarify the nature and functions of these putative neuropeptide transcripts, to assign it a possible role in the anorexic behavior observed in fertilized females. Furthermore, the response to starvation was another gene ontology enriched in females. This response may be related to the egg protecting Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase behavior during pre and post spawning phases [70,71]. As other incirrate octopods, O. maya females decrease and stop feeding during incubation of the eggs and focus exclusively on their care until the hatching of the juveniles [70]. Maternal care generally includes the protection of the egg mass from potential predators, ventilation by flushing water through the eggs, cleaning the surface of the eggs, and removing dead embryos [72]. These results suggest that WB may contribute to the regulation of this anorexic behavior in adult O. maya females. Similarly, the response to cortisol was another enriched GO in females, which may be linked to their anorexic behavior. Food deprivation has been correlated with high levels of cortisol in fishes [73][74][75]. A similar glucocorticoid analogous to cortisol, the corticosterone has been detected in Enteroctopus dofleini feces after stressful events [76]; so, it is possible that this molecule is released in response to stress conditions [12] also in O. maya with higher expression in females due to its anorexic behavior before spawning. The presence in the WB of the gene encoding for corticosteroid 11-beta-dehydrogenase isozyme 2 (HSD11B2), which Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase participates in the corticosterone inactivation [77,78], together with the corticotropin-releasing factor receptor 2 (CRHR2), which expression was higher in females, supports this hypothesis. These results suggest that the WB plays a role in the glucocorticoid metabolism, and in the response to glucocorticoids. In this regard, the WB can be an important target of glucocorticoids, considering that these molecules can restrain the immune and inflammatory responses [79][80][81][82][83], which are mediated by the WB at least partially [8,84]. But why the immune response should be down-regulated? In fertilized females this regulation could be required to permit the maintenance of foreign cells, specifically the spermatozoa, which are stored in the females' oviductal gland up to four months [70,[85][86][87]. The immune regulation leaded by glucocorticoids with possible action in the WB, may contribute to keep spermatozoa safe from the females' immune system. This down-regulation over females' immune response may be linked to the high expression of signaling genes involved in bacterial and viral infections, as illustrated by the enriched KEEG pathways. On the other hand, despite some evidences of steroids Values represent the fold change (log 10 ) of each target vs the reference genes. Reference genes: HNRNPD and VATD. Samples: FWB = female white body, FH = female heart, MWB = male white body, MH = male heart, MT = male testis. The 95% confidence interval of each group is shown and the significant differences among groups are represented with different letters. https://doi.org/10.1371/journal.pone.0216982.g005 Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase hormones pathways in the nervous system and gonads in O. vulgaris [88,89], the precise gland or nervous lobe were the octopus' glucocorticoids are synthesized remain unclear.
In contrast to females, different processes were enriched by the highly expressed transcripts in males' WB. Transcripts from the androgen signaling pathway, encoding proteins such as nuclear receptor activator 4, prohibitin and protein DJ-1, were significantly more abundant in males. This result together with the higher expression of the estrogen receptor in females, support the idea that steroid hormones are involved in female and male physiological dimorphism during reproduction [89]. Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase Regarding the enriched KEGG pathways, apparently the males suffered a mitochondrial dysfunction, resulting in higher production of reactive oxygen species (ROS), in a similar way to that observed in the Huntington disease pathway (S16 Fig). Moreover, ROS metabolism was also enriched in males by the high expression of transcripts encoding peroxiredoxin-1, peroxiredoxin-4, superoxide dismutase, autophagy protein 5 and NADH dehydrogenase [ubiquinone] iron-sulfur protein 3. Notably, this apparent mitochondrial dysfunction leading to high ROS production and subsequent induction of the antioxidant system, could be an adaptation for the defense against pathogens. This system was also well represented in the transcriptome of O. vulgaris hemocytes [19], which in turn are originated in the WB [10]. Authors have suggested that the antioxidant system enzymes may play a role in the defense against pathogens by the hemocytes [19,90,91]. On the other hand, the apoptotic process was also conspicuous in the O. vulgaris hemocyte transcriptome, with high expression of initiator and effector proteins for apoptosis after parasitic infection [19]. This mechanism is also a major defense against pathogens [92]. In this study, the regulation of apoptosis was an enriched process in males' Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase WB by the high expression of transcripts coding for TNF receptor-associated factor 2, apoptotic chromatin condensation inducer in the nucleus, Bcl-2 homologous antagonist/killer (apoptosis regulator BAK), geranylgeranyl transferase type-2 subunit beta, RNA-binding protein 5, ribosomal L1 domain-containing protein 1, and 40S ribosomal protein S3, which were classified as apoptotic inducers according to the GO. In this regard, the higher expression of antioxidant enzymes and the positive regulation of apoptosis in males' WB could be linked to a higher production of hemocytes compared to females.
In the case of unigenes detected exclusively in males, they were classified mainly in microtubule-based process and cilium morphogenesis, despite there is no evidence so far of ciliated cells within the WB [4,10,11]. However, this organ is tightly connected to nervous tissue, especially to optic tracts and optic lobes, where microtubule-based processes are essential during neurogenesis, playing a role in the organization and dynamics of axons and dendrites [93]. In this regard, the WB may be involved in neurogenesis, possibly due to its high stem cell content, Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase and its proximity to nervous tissue. Recently, Bertapelle et al. [94] detected neurogenesis in the CNS of adult O. vulgaris, which is the first report of adult neurogenesis in lophotrochozoan animals. By contrast, these cilium-related unigenes were not detected in the females' WB. This can be due to a reduction in neurogenesis because of a more advanced senescence in females compared to males.
Finally, in males' WB there were enrichments for spermatogenesis-related GO terms, however this could be a coincidence causing an incorrect GO assignment, considering that the molecular basis of the stem cell system of hematopoiesis and spermatogenesis appears to be very similar [95,96]. This incorrect GO assignment was confirmed by the qPCR results, where the expression levels of transcripts involved in spermatogenesis (VASA homolog and PIWIL1) were compared to those in testis. We observed that despite their presence in the white body, their expression was significantly higher in testis. However, ZAN and SPAG9, which are involved in sperm-egg interaction [59,60], were highly expressed in both female and male white bodies, ZAN showing similar levels to those detected in testis and SPAG9 even higher Sex-specific transcriptomic profiles in the white body of Octopus maya during reproductive phase than testis. These proteins could be synthesized by independent tissues (e.g. WB and gonads) to enhance fertilization, but further research is needed to test this hypothesis.

Conclusion
The results obtained in this study, are evidences of the involvement of O. maya WB in hematopoiesis and in the regulation of immune processes. Notably, there was an important differentiation of signaling pathways between female and male white bodies. Multiple signaling cascades were upregulated in females, some of them related to bacterial and viral infections. At the same time, females showed higher expression of unigenes related to neuropeptide signaling pathways, as well as unigenes involved in the response to glucocorticoids and starvation. By contrast, in males we detected higher expression of unigenes required in the androgen signaling pathway, antioxidant response, and apoptosis. Considering that glucocorticoids can suppress immunity and the antioxidant/apoptotic response can enhance the defense against pathogens; we can infer that immune response was down-regulated in fertilized females and up-regulated in mature males, which is congruent with the higher gene expression related to infection pathways detected in females. Furthermore, our data suggest for the first time, an involvement of this organ in the physiological differences between mature males and females, showing differential gene expression processes in a sex-specific way during the reproductive phase.