biology

Article
De novo Assembly and Analysis of Tissue-Specific
Transcriptomes of the Edible Red Sea Urchin Loxechinus albus
Using RNA-Seq
Paulette Antiqueo 1,2 , Rodrigo Zuloaga 1,2 , Macarena Bastias-Molina 1,3 , Claudio Meneses 1,3 ,
Juan Manuel Estrada 4 , Alfredo Molina 1,2,4 and Juan Antonio Valdés 1,2,4, *

1 Departamento Ciencias Biológicas, Facultad de Ciencias de la Vida, Universidad Andres Bello,

Santiago 8370186, Chile; [email protected] (P.A.); [email protected] (R.Z.);
[email protected] (M.B.-M.); [email protected] (C.M.); [email protected] (A.M.)
2 Interdisciplinary Center for Aquaculture Research (INCAR), Concepción 4030000, Chile
3 Centro de Biotecnología Vegetal, FONDAP Center for Genome Regulation, Facultad de Ciencias de la Vida,
Universidad Andres Bello, Santiago 8370146, Chile
4 Centro de Investigación Marina Quintay (CIMARQ), Facultad de Ciencias de la Vida,
Universidad Andres Bello, Valparaíso 2340000, Chile; [email protected]
* Correspondence: [email protected]; Tel.: +56-2661-8363; Fax: +56-2661-8415

Simple Summary: Edible red sea urchin (Loxechinus albus) is an endemic species of echinoderm
distributed along the Chilean coasts. This resource has been overexploited in recent years, depleting
their natural populations. At present, there are few reported gene sequences available in public



databases, restricting the molecular studies associated with aquaculture for this species. The aim of


this study was to present the first annotated reference transcriptome of L. albus using NGS technolo-
Citation: Antiqueo, P.; Zuloaga, R.;
gies and the differential expression transcripts analysis of the evaluated tissues. The transcriptome
Bastias-Molina, M.; Meneses, C.;
data obtained in this study will serve as a reference for future molecular research in the edible red
Estrada, J.M.; Molina, A.; Valdés, J.A.
sea urchin and other sea urchin species.
De novo Assembly and Analysis of
Tissue-Specific Transcriptomes of the
Abstract: Edible red sea urchin (Loxechinus albus) is an endemic echinoderm species of the Chilean
Edible Red Sea Urchin Loxechinus
albus Using RNA-Seq. Biology 2021,
coasts. The worldwide demand for high-quality gonads of this species has addressed the depletion
10, 995. https://doi.org/10.3390/ of its natural populations. Studies on this sea urchin are limited, and genomic information is
biology10100995 almost nonexistent. Hence, generate a transcriptome is crucial information that will considerably
enrich molecular data and promote future findings for the L. albus aquaculture. Here, we obtained
Academic Editor: Patricia Pereiro transcriptomic data of the edible red sea urchin by Illumina platform. Total RNA was extracted
from gonads, intestines, and coelomocytes of juvenile urchins, and samples were sequenced using
Received: 25 August 2021 MiSeq Illumina technology. A total of 91,119,300 paired-end reads were de novo assembled, 185,239
Accepted: 29 September 2021 transcripts produced, and a reference transcriptome created with 38.8% GC content and an N50 of
Published: 2 October 2021
1769 bp. Gene ontology analysis revealed notable differences in the expression profiles between
gonads, intestines, and coelomocytes, allowing the detection of transcripts associated with specific
Publisher’s Note: MDPI stays neutral
biological processes and KEGG pathways. These data were validated using 12 candidate transcripts
with regard to jurisdictional claims in
by real-time qPCR. This dataset will provide a valuable molecular resource for L. albus and other
published maps and institutional affil-
species of sea urchins.
iations.

Keywords: edible red sea urchin; Loxechinus albus; RNA-seq; reference transcriptome

Copyright: © 2021 by the authors.

Licensee MDPI, Basel, Switzerland.
1. Introduction
This article is an open access article
distributed under the terms and The Loxechinus albus (Molina, 1782), or edible red sea urchin, is an echinoderm species
conditions of the Creative Commons of the Chilean and Peruvian coasts, distributed along ca. Cape Horn, Chile (56◦ 700 S) to the
Attribution (CC BY) license (https:// Isla Lobos de Afuera, Peru (6◦ 530 S) [1]. The worldwide demand for high-quality gonads
creativecommons.org/licenses/by/ of this sea urchin has addressed a vast overexploitation of its natural populations [2]. Har-
4.0/). vesting of L. albus represents the major sea urchin fishery among world urchin fisheries [3].

Biology 2021, 10, 995. https://doi.org/10.3390/biology10100995 https://www.mdpi.com/journal/biology

Biology 2021, 10, 995 2 of 19

The aquaculture of this species, involving the rearing tank production of larvae, juvenile,
and later fattening in natural environments, are important approaches to aquaculture
diversification in Chile and to restore the overexploited coastal areas [4].
One of the main difficulties in the study of biological and molecular mechanisms
associated with the farming of this species is the limited genomic information available [5,6].
In this context, transcriptome sequencing is useful to identify genes participating certain
biological processes when genomic data are not available [7]. This analysis allows a broad
comprehension of molecular mechanisms involved in biological processes from data on
predicted function of genes [8]. Progress in the characterization of the transcriptome
in commercial sea urchins is achievable due to advances in next-generation sequencing
(NGS) technologies. NGS has allowed the research of sea urchin transcriptomes and
other non-model species in brief periods of time at a low cost [9–11]. The molecular
information achieved has provided significant value regarding the physiological responses
to adaptation in a variety of commercial sea urchins under fluctuating environmental
conditions [12,13].
At this time, the existing information on L. albus biology is limited and is related to
with oxidative metabolism [14], growth patterns [15], the performance of early juveniles un-
der food type and feeding frequency [16], and cryopreservation of embryos and larvae [17].
However, biological studies with molecular bases carried out in this species are scarce,
mainly due to the low amount of genomic information available [11,18]. Although some
advances have been made in the transcriptome characterization and mitogenome of this
species in recent years, the low coverage of the technology used, as well as the use of gonads
as the only target tissue, has limited the obtainment of a high-quality reference transcrip-
tome [5,6,9,19]. Therefore, we present here the first annotated transcriptome of juvenile
edible red sea urchin using NGS technologies based on three critical tissues for physiologi-
cal homeostasis of echinoderms and the expression analysis of the transcripts present in
each tissue: (i) gonads, involved in reproduction and exportation product for aquaculture,
(ii) intestine, involved in food digestion and nutrient uptake, and (iii) coelomocytes, in-
volved mainly in immune surveillance and inflammatory process. The transcriptome data
obtained here will provide a reference for molecular studies in the farming of L. albus and
other sea urchin species.

2. Materials and Methods

2.1. Experimental Design and Sampling
Loxechinus albus specimens were obtained from the Centro de Investigación Marina
de Quintay (CIMARQ; 33◦ 130 S, 71◦ 380 O, Valparaiso, Chile). Briefly, fertilization was
performed using a pool of gametes from four females and four males stimulated to spawn
by injection of 3 mL of 0.5 M KCl. The embryos generated were cultured in 200 L larval
rearing containers and larvae developed were fed with Chaetoceros gracilis microalgae. The
larvae were grown in 50 L tanks in filtrated and aerated seawater and then preconditioned
to settle in post-larval condition. Juvenile sexually immature sea urchins were maintained
under natural conditions (13 ± 1 ◦ C) in the spring season. The sea urchins were three
years old and weighed 30 ± 5 g. The animals were fed macroalgae ad libitum (Lessonia
sp., Macrocystis sp., Durvillea sp.). A total of 10 sea urchins were selected, dissected, and
three different tissues were collected: intestines, gonads, and coelomocytes. Intestines
were cleaned with phosphate buffer solution (PBS 1×) before storage. In immature gonads,
germ cells were undifferentiated, revealing no sex differentiation. The coelomic fluid
was collected by cutting the peristomal membrane, mixed with anticoagulant (20 mM
Tris–HCl, 0.5 M NaCl, and 30 mM EDTA; pH 7.4), centrifuged for 5 min at 5000× g, and
then coelomocyte pellet was collected. Samples were rapidly frozen in liquid nitrogen and
deposited at −80 ◦ C until use.
Biology 2021, 10, 995 3 of 19

2.2. Isolation of RNA and Sequencing

Total RNA was obtained using columns of the RNeasy Mini Kit (Qiagen, Austin,
TX, USA). The genomic DNA from RNA samples with removed by DNase I treatment.
RNA was quantified by fluorometry using a Qubit 2.0 Fluorometer (Life Technology,
Carlsbad, CA, USA), and the integrity of RNA was measured using the Fragment Analyzer
(Analytical Advanced Technologies, Ames, IA, USA). Total RNA from five sea urchins
were pooled in equal quantities by tissue, in duplicate, and then used to mRNA libraries
construction. These libraries were generated by the TruSeq RNA Sample Preparation Kit v2
(Illumina, San Diego, CA, USA). Finally, libraries were sequenced (2 × 250 bp) utilizing the
MiSeq technology (Illumina) at the Center for Plant Biotechnology (Universidad Andrés
Bello, Santiago, Chile). The raw reads of the present study were uploaded to the NCBI SRA
database under BioProject PRJNA475570, with accession number SRP150640.

2.3. Processing of Raw Data, De novo Assembly, and Validation of Assembly

First, the raw sequence reads were quality checked using FASTQC software. Adapters
were removed, and raw data were trimmed using FlexBar [20] with Phred scores below 38
and 250 bp reads. The de novo transcriptome was assembled using all libraries (two libraries
per tissue) with the Trinity program using default parameters [21]. Transcripts were filtered
based on the minimal number of mapped reads with the Corset program using default
parameters [22]. To evaluate de novo assembly integrity, the assembled transcriptome
by Benchmarking Universal Single-Copy Orthologs (BUSCO) was compared against the
OrthoDBv9 database (Vertebrata and Eukaryota) to identify orthologous genes that were
highly conserved [23].

2.4. Functional Annotation and Analysis of Differentially Expressed Transcripts

For transcriptome annotation, a search in BlastX against the UniProt (https://www.
uniprot.org/blast/; accessed 18 March 2021), Nonredundant (NR; https://blast.ncbi.nlm.
nih.gov/Blast.cgi; accessed 18 March 2021), and Clusters of orthologous groups for eukary-
otic complete genomes (COG; https://www.ncbi.nlm.nih.gov/research/cog; accessed 18
March 2021) databases was performed. The Blast2GO program was employed to acquire
gene ontology (GO) annotation [24], and the WEGO application [25] was used to carry
out GO functional classification for all transcripts. Recognition of differentially expressed
transcripts (DETs) from the gonads, intestines, and coelomocytes were realized using
Bowtie by mapping against the assembled L. albus transcriptome [26]. The RSEM soft-
ware was used to assess expression values of fragments per kilobase million (FPKM) [27].
EdgeR was employed to determine differential expression between intestine vs. gonad,
coelomocytes vs. intestine, and coelomocytes vs. gonads [28]. Transcripts detected with
false discovery rate (FDR)-corrected p values < 0.001 and absolute values of fold-change >
4.0 were incorporated in the GO and Kyoto encyclopedia of genes and genomes (KEGG)
enrichment analyses.

2.5. Gene Ontology and KEGG Enrichment Analysis

The DETs were examined against the DAVID resource [29] and then categorized based
on GO terms for molecular functions, biological processes, cellular components, and KEGG
pathways. To determine a relationship between the DAVID background and L. albus DETs,
a search in BLASTx was performed against Strongylocentrotus purpuratus Ensembl proteins
for major matches with the L. albus transcriptome. Ensembl Gene IDs of S. purpuratus were
acquired from the resultant Ensembl protein entries. Custom IDs set were selected for
DAVID analysis as the “Background” Standard settings for ease (0.1) and gene count (2).
The cut off p value used for molecular functions and cellular components was 1 × 10−3 ,
and for biological processes was 1 × 10−6 .
Biology 2021,10,
Biology2021, 10,995
x 4 4ofof1920

2.6. Validation of RNA-Seq

All quantitative by Real-Time
real-time qPCR chain reaction (qPCR) assays were performed
polymerase
according to MIQE recommendations
All quantitative real-time polymerase [30]. Total
chain RNA (qPCR)
reaction isolation fromwere
assays gonads, intestines,
performed ac-
cording to MIQE recommendations [30]. Total RNA isolation from gonads, intestines,quan-
and coelomocytes was realized using columns of RNeasy Mini Kit (Qiagen). RNA and
tification waswas
coelomocytes measured
realizedby NanoDrop
using columnstechnology
of RNeasy MiniwithKitan(Qiagen).
Epoch Multivolume Spectro-
RNA quantification
photometer
was measured System (BioTek, technology
by NanoDrop Winooski, VT, withUSA). For complementary
an Epoch DNA (cDNA) syn-
Multivolume Spectrophotometer
thesis, only
System RNAWinooski,
(BioTek, with an A260/280
VT, USA). ratio
Forbetween 1.9 and 2.1
complementary DNAwas(cDNA)
selected.synthesis,
This procedure
only
was performed
RNA with an A260/280using 1ratio
μg of RNA by
between 1.9QuantiTect
and 2.1 wasReverse
selected.Transcription
This procedureKit was
(Qiagen),
per-
eliminating first genomic DNA with the wipeout buffer included and
formed using 1 µg of RNA by QuantiTect Reverse Transcription Kit (Qiagen), eliminating then reverse tran-
scribed
first into cDNA
genomic DNA withat 42 the
°C for 30 min.
wipeout qPCR
buffer reactionand
included wasthen
made usingtranscribed
reverse Brilliant II SYBR
into
cDNA
GreenatQPCR 42 ◦ CMaster
for 30 min.
Mix qPCR
(Agilentreaction was made
Technologies, using
Santa Brilliant
Clara, CA,IIUSA),
SYBRaGreen QPCRof
1/dilution
Master
cDNA,Mix and(Agilent Technologies,
5 μM primers Santa Clara,Table
(Supplementary CA, USA), a 1/dilution
S1). Real-time of cDNA,
PCRs were and
run 5on µM an
primers
Mx3000P (Supplementary
qPCR SystemTable S1). Technologies)
(Agilent Real-time PCRsinwere run onThe
triplicate. an Mx3000P qPCR System
PCR amplification pro-
(Agilent Technologies)
gram included 95 °C forin10triplicate. The PCR
min; 40 cycles amplification
of 95 °C for 30 s, Tm program 95 °C for9532◦ C
included
for 30 s, s, for
and
10 min; 40 cycles of 95 ◦ C for 30 s, Tm for 30 s, 95 ◦ C for 32 s, and 72 ◦ C for 30 s. The QGene
72 °C for 30 s. The QGene application was utilized to analyze gene expression [31] and
application was utilizedwith
data were normalized to analyze gene expression
18S ribosomal subunit RNA[31] and data were normalized
as housekeeping gene. with
18S ribosomal subunit RNA as housekeeping gene.
2.7. Statistical Analysis
2.7. Statistical Analysis
All statistical analyses were performed using GraphPad Prism v8.0 (GraphPad Soft-
ware,AllLastatistical
Jolla, CA,analyses
USA). The were performed
generated datausing
wereGraphPad
presented Prism v8.0 (GraphPad
as the mean ± standardSoft-
error
ware,
of the mean (SEM). A one-way analysis of variance (ANOVA) followed bystandard
La Jolla, CA, USA). The generated data were presented as the mean ± error
the Bonferroni
of the mean (SEM). A one-way analysis of variance (ANOVA) followed by the Bonferroni
post hoc test was used to differentiate the means between groups. A value of p < 0.001 was
post hoc test was used to differentiate the means between groups. A value of p < 0.001 was
accepted as significative data. Finally, correlation analysis between real-time qPCR and
accepted as significative data. Finally, correlation analysis between real-time qPCR and
RNA-seq were evaluated by multiple linear regression through p value and coefficients of
RNA-seq were evaluated by multiple linear regression through p value and coefficients of
determination (R2 2).
determination (R ).

3.3.Results
Results
3.1.Raw
3.1. RawData
DataSequencing,
Sequencing,De
Denovo
novoAssembly
AssemblyTranscriptome
Transcriptomeand
andFunctional
FunctionalAnnotation
Annotation
ThecDNA
The cDNAlibraries
librarieswere
weregenerated
generatedfrom
fromintestines,
intestines,gonads,
gonads,and
andcoelomocytes
coelomocytesofof
poolededible
pooled ediblered
redsea
seaurchins
urchins(Figure
(Figure1).
1).

Figure1.1.Juvenile
Figure Juvenileedible
ediblered
redsea
seaurchins
urchins(Loxechinus
(Loxechinusalbus).
albus).
 Biology 2021, 10, 995 5 of 19

The Illumina MiSeq sequencing produced 95,745,640 paired end reads of cDNA library
replicates for each tissue (Table 1). The obtained raw data were trimmed by eliminating
adapters, contaminant sequences and filtering base pairs with low-quality, the high-quality
reads were reduced to 91,119,300 base pairs (Table 1). The Trinity software was used to
de novo assembly using all libraries, resulting in 278,803 transcripts. The high-quality reads
were mapped against the transcriptome generated and reduced by Corset software. The
newly reduced transcriptome had 185,239 transcripts, with N50 = 1769 bp and GC% = 38.81
(Table 1). The length distribution of the obtained transcripts was detailed in Supplementary
Figure S1. The BUSCO database was selected to compare the assembled transcriptome,
which includes information of orthologous genes highly conserved. The assembly found
248 BUSCO genes of the eukaryotic core gene, 222 complete (89.5%), 8 fragmented (3.2%),
and 18 missing genes (7.3%). The annotation of L. albus transcripts was carried out by
BlastX searches against the non-redundant (NR), UniProt, and cluster of orthologous
groups (COG) databases. The statistics of similarity and target species of all transcripts are
shown in Figure 2, revealing that 91.2% of sequences matched S. purpuratus.

Table 1. Summary of transcriptome sequencing for the edible red sea urchin (L. albus) tissue and
y 2021, 10, x assembly statistics. bp: base pair. 5 of 20

Name Number of Reads Number of Reads after Trim

Coelomocyte 20,682,190 19,948,624
The Illumina MiSeq sequencing produced 95,745,640 paired end reads of cDNA li-
Coelomocyte (replicate) 16,865,448 15,699,186
brary replicates for eachIntestine
tissue (Table 1). The obtained raw data
12,145,212 were trimmed by elimi-
11,748,696
nating adapters, contaminant sequences and filtering
Intestine (replicate) base pairs with low-quality,
11,348,164 the
10,620,928
high-quality reads were Gonad
reduced to 91,119,300 base pairs (Table 1). The Trinity software
18,495,858 17,901,954
Gonad (replicate)
was used to de novo assembly 16,208,768
using all libraries, resulting in 278,803 transcripts. The15,199,912
high-
Total 95,745,640
quality reads were mapped against the transcriptome generated and reduced by Corset 91,119,300
software. The newly reduced transcriptome had De
Transcriptome
185,239 transcripts, with N50 = 1769
novo Assembly Afterbp Filter
and GC% = 38.81 (Table 1). The length distribution of the (Trinity) (Corset)
obtained transcripts was detailed
in Supplementary Figure S1.contigs
Total The BUSCO database was278,803 selected to compare the assembled 185,239
Average
transcriptome, which large contig
includes (bp)
information of orthologous 326 genes highly conserved. The 929
assembly found 248 BUSCO genes of the eukaryotic core gene, 222 complete (89.5%),- 8
Coverage contig 708.32
fragmented (3.2%), and 18%GC 38.2
missing genes (7.3%). The annotation of L. albus transcripts 38.81
was
N10 (bp) 5015 5328
carried out by BlastX searches
N30 (bp) against the non-redundant 2645 (NR), UniProt, and cluster of
2945
orthologous groups (COG) databases.
N50 (bp) The statistics of similarity
1418 and target species of all
1769
transcripts are shown inTotal
Figure 2, revealing that 91.2%197,480,887
bases of sequences matched S. purpuratus.
172,122,576

Figure 2. A species-based
Figure BlastX comparative
2. A species-based analysis revealed
BlastX comparative analysisthe major match
revealed Strongylocentrotus
withmatch
the major purpuratus.
with Strongylocen-
trotus purpuratus.

Table 1. Summary of transcriptome sequencing for the edible red sea urchin (L. albus) tissue and
assembly statistics. bp: base pair.

Name Number of Reads Number of Reads after Trim

Biology
Biology2021,
2021,10,
10,x 995 6 6ofof2019

The
Thefunctional
functionalannotation
annotationofoftranscripts
transcriptswas wasrealized
realizedusing
usingthe
theBlast2GO
Blast2GOplatform,
platform,
applying a GO term search through transcripts with BLAST
applying a GO term search through transcripts with BLAST hits matched hits matched against NR da-
against NR
tabase.
database. We detected significant similarity with a total of 57,106 (31%) transcripts. TheGO
We detected significant similarity with a total of 57,106 (31%) transcripts. The GO
analysis
analysisrevealed
revealed38,265
38,265 GO
GO outcomes
outcomes for for biological
biologicalprocesses
processes(20.2%),
(20.2%),36,046
36,046forfor molec-
molecular
ular functions
functions (19.1%)
(19.1%) andand 35,909
35,909 GOcellular
GO for for cellular components
components (19.0%).
(19.0%). A significant
A significant pro-
proportion
portion of the annotated transcripts in biological process were assigned
of the annotated transcripts in biological process were assigned to cellular componentto cellular com-
ponent biogenesis
biogenesis (GO:0044085)
(GO:0044085) and cellular
and cellular localization
localization (GO: 0051641)
(GO:0051641) (Figure(Figure
3). 3).

Figure
Figure3.3.GO
GOfunctional
functional classification
classificationassigned
assignedthe
themost
mostpercentage
percentageofofthe
theannotated
annotatedtranscripts
transcriptstotocellular
cellularcomponent
component
biogenesis
biogenesis term forfor
term biological process;
biological cytosol
process; term for
cytosol cellular
term component;
for cellular and small
component; molecule
and small binding
molecule term for molecular
binding term for
function,
molecular respectively. Analysis wasAnalysis
function, respectively. carried was
out with theout
carried WEGO
with program
the WEGO forprogram
the edible
forsea
theurchin
edible (L.
seaalbus)
urchinreference
(L. albus)
transcriptome.
reference transcriptome.

For
Forcellular
cellularcomponents
componentsand andmolecular
molecularfunctions,
functions,several
severalannotated
annotatedtranscripts
transcriptswere
were
allocated
allocatedtotocytosol
cytosol(GO: 0005829) and
(GO:0005829) andnucleoplasm
nucleoplasm(GO:0005654);
(GO:0005654);and andnucleotide
nucleotidebinding
binding
(GO:0000166),
(GO:0000166),nucleoside
nucleosidephosphate
phosphatebinding
binding(GO:1901265),
(GO:1901265),and andsmall
smallmolecule
moleculebinding
binding
(GO: 0036094) terms,
(GO:0036094) terms,respectively
respectively(Figure
(Figure3).3). Based
Based on sequence
on sequence homology,
homology, 32,23132,231 se-
sequences
quences were classified
were classified into 25 into 25 functional
functional categories
categories (Figure(Figure
4). The4).most
The represented
most represented cat-
categories
egories were General
were General Functional
Functional Prediction
Prediction only, followed
only, followed by Signal
by Signal Transduction.
Transduction. These These
results
results indicate
indicate that wethat we generated
generated a reference
a reference transcriptome
transcriptome for thesea
for the edible edible seabased
urchin urchin
on
based on 91,119,300
91,119,300 high-quality
high-quality reads
reads that that assembled
assembled de novo de novo
into into 185,239
185,239 transcripts
transcripts with anwith
N50
an
ofN50
1769ofbp
1769
andbp and GC
38.8% 38.8% GC content.
content.
 Biology 2021, 10, x 7 of 20
Biology 2021, 10, 995 7 of 19

Figure 4. The General Functional Prediction only was the most represented category for the 25 cluster of orthologous
Figure 4. The General Functional Prediction only was the most represented category for the 25 cluster of orthologous
groups (COG) functional classification in the L. albus reference transcriptome.
groups (COG) functional classification in the L. albus reference transcriptome.
3.2. Assessment of Differentially Expressed Transcripts
3.2. Assessment of Differentially Expressed Transcripts
About ~95.0% of the reads were mapped to the assembled transcriptome by Bowtie,
About
and then RSEM~95.0%
and of the reads
EdgeR werewere
usedmapped to the
to estimate theassembled
expressiontranscriptome
values and toby Bowtie,
calculate
and then RSEM and EdgeR were used to estimate the expression values and
differential expression transcripts (DETs) between tissues, respectively. Transcripts withto calculate
differential
false discoveryexpression transcripts (DETs)
rate (FDR)-corrected p valuesbetween tissues,
< 0.001 and respectively.
absolute Transcripts
fold-change values >with
4.0
falsedefined
were discoveryas rate
DETs.(FDR)-corrected
The mapping pofvalues < 0.001
L. albus and absolute
transcriptome fold-change
indicated valuestran-
that 26,864 > 4.0
were defined
scripts as DETs. The
were differentially mappingamong
expressed of L. albus transcriptome
coelomocytes, indicated
intestine, andthat 26,864
gonad. tran-
These
scripts
DETs were
were differentially
clustered expressed
using hierarchy byamong coelomocytes,
comparisons between intestine,
patterns ofand
genegonad. These
expression
DETs were
(Figure 5). clustered using hierarchy by comparisons between patterns of gene expression
(Figure
The5).
correlation analysis of heat map clustered samples together in the same group,
representing a high reproducibility of RNA-seq data among replicates (Supplementary
Figure S2). To identified DETs, tissues were exposed to a succession of paired comparations
(Figure 6): In coelomocytes vs. gonad, 16,406 DETs were identified (8394 upregulated and
8012 downregulated); in coelomocytes vs. intestine, 14,025 DETs were identified (5558
upregulated and 8467 downregulated); intestine vs. gonad, 10,015 DETs were identified
(6179 upregulated and 3836 downregulated).
 Biology 2021, 10, 995 8 of 19
Biology 2021, 10, x 8 of 20

Biology 2021, 10, x Figure The

Figure5.5.The heat
heat mapmap showed
showed differentiated
differentiated clustering
clustering of expressed
of expressed transcripts
transcripts according
according 9to 20to
of re-
respective
spective expression
expression values
values across
across tissues.
tissues. Parameters:
Parameters: fold change
fold change (absolute
(absolute valuesvalues
> 4.0) > 4.0)FDR
and and
corrected p valuep(p
FDR corrected < 0.001).
value Abbreviations:
(p < 0.001). CL: coelomocyte,
Abbreviations: INT: intestine,
CL: coelomocyte, GN: gonad.
INT: intestine, GN: R1: Repli-
gonad. R1:
cate 1, R2: Replicate 2.
Replicate 1, R2: Replicate 2.

The correlation analysis of heat map clustered samples together in the same group,
representing a high reproducibility of RNA-seq data among replicates (Supplementary
Figure S2). To identified DETs, tissues were exposed to a succession of paired compara-
tions (Figure 6): In coelomocytes vs. gonad, 16,406 DETs were identified (8394 upregu-
lated and 8012 downregulated); in coelomocytes vs. intestine, 14,025 DETs were identified
(5558 upregulated and 8467 downregulated); intestine vs. gonad, 10,015 DETs were iden-
tified (6179 upregulated and 3836 downregulated).

Figure 6. Paired comparisons of differentially expressed transcripts between tissues of L. albus.

Figure 6. Pairedand
Upregulated comparisons of differentially
downregulated transcriptsexpressed
measuredtranscripts between
are indicated tissues
by black andofgray
L. albus. Up-
columns,
regulated and downregulated transcripts measured are indicated by black and gray columns, re-
respectively. Parameters: fold change (absolute values > 4.0) and FDR corrected p value (p < 0.001).
spectively. Parameters: fold change (absolute values > 4.0) and FDR corrected p value (p < 0.001).
Abbreviations: CL: coelomocyte, INT: intestine, GN: gonad.
Abbreviations: CL: coelomocyte, INT: intestine, GN: gonad.
The analysis of Venn diagram revealed that 737 DETs (2.7%) were commonly ex-
The analysis
pressed among theof Venn
three diagram revealed
tissues (Figure 7). that 7377061
Also, DETs (2.7%)
DETs werewere
(26.3%) commonly ex-
specifically
pressed among the three tissues (Figure 7). Also, 7061 DETs (26.3%) were specifically
expressed between coelomocytes and gonad, specific 4280 DETs (15.9%) among coelomo- ex-
pressed between coelomocytes and gonad, specific 4280 DETs (15.9%) among coelomo-
cytes and intestine, and 2678 DETs (10%) between intestine and gonad. A complete list of
the DETs for each tissue was included in Supplementary Table S2. The tissue-specific
DETs were examined in DAVID resource and categorized as biological processes, molec-
 Biology 2021, 10, 995 9 of 19

cytes and intestine, and 2678 DETs (10%) between intestine and gonad. A complete list
of the DETs for each tissue was included in Supplementary Table S2. The tissue-specific
DETs were examined in DAVID resource and categorized as biological processes, molecular
functions, and cellular components. The coelomocyte DETs were considerably enriched
in positive regulation of apoptotic process (GO:0043065) and intracellular signal trans-
duction (GO:0035556) for biological processes (Figure 8a). The most enriched GO terms
for molecular functions and cellular components were assigned to poly(A) RNA bind-
ing (GO:0044822) and ATP binding (GO:0005524), and cytosol (GO:0005829), respectively
(Supplementary Figure S3). In contrast, DNA repair (GO:0006281) and microtubule-based
process (GO:0007017) were the major enriched biological processes in gonads (Figure 8b).
The main transcripts assigned to molecular functions and cellular components were micro-
tubule motor activity (GO:0003777) and ATP binding (GO:0005524) GO terms, and dynein
complex and cilium (GO:0005929) GO terms, respectively (Supplementary Figure S4). For
biological processes, transmembrane transport (GO:0055085) and microtubule-based pro-
cess (GO:0007017) were the most enriched in intestinal DETs (Figure 8c). The majority
21, 10, x
assigned transcripts for molecular functions and cellular components were GTP binding
10 of 20
(GO:0005525) and motor activity (GO:0003774), and extracellular exosome (GO:0070062)
and membrane (GO:0016020), respectively (Supplementary Figure S5).

Figure 7. Differentially expressed transcripts between gonads, intestines, and coelomocytes of the edible sea urchin. Venn
diagram ofFigure 7. Differentially
the total differentially expressed transcripts in
expressed transcripts between
L. albus.gonads, intestines,
Each color and
indicates thecoelomocytes of the tissue and
comparison between
edible
the numbers sea urchin.
of genes Venn
that were diagram ofexpressed.
differentially the total Parameters:
differentiallyfold
expressed transcripts
change (absolute in L.>albus.
values Each
4.0) and FDR corrected
color indicates the comparison between tissue and the numbers of genes that were differentially
p value (p < 0.001). Abbreviations: CL: coelomocyte, INT: intestine, GN: gonad.
expressed. Parameters: fold change (absolute values > 4.0) and FDR corrected p value (p < 0.001).
Abbreviations: CL: coelomocyte, INT: intestine, GN: gonad.
Biology
Biology2021, 10,10,
2021, 995x 1011ofof1920

Figure8.8.The
Figure TheTop-16
Top-16Gene
GeneOntology
Ontologybiological
biologicalprocess
process(BP)
(BP)enrichment
enrichmentofofup-regulated
up-regulatedtranscripts
transcriptsfrom
fromL.L.albus
albustissues.
tissues.
(a) In coelomocytes compared to intestine (CL vs. INT), the most enriched term was positive regulation of apoptotic
(a) In coelomocytes compared to intestine (CL vs. INT), the most enriched term was positive regulation of apoptotic process, pro-
cess, and compared to the gonad (CL vs. GN) was intracellular signal transduction; (b) in the gonad compared
and compared to the gonad (CL vs. GN) was intracellular signal transduction; (b) in the gonad compared to intestine (GN to intestine
(GN vs. INT), the most enriched term was DNA repair, and compared to coelomocytes (GN vs. CL) was microtubule-
vs. INT), the most enriched term was DNA repair, and compared to coelomocytes (GN vs. CL) was microtubule-based
based process; (c) in intestine compared to coelomocytes (INT vs. CL) the most enriched term was transmembrane
process; (c) in intestine compared to coelomocytes (INT vs. CL) the most enriched term was transmembrane transport, and
transport, and compared to the gonad (INT vs. GN) was microtubule-based process.
compared to the gonad (INT vs. GN) was microtubule-based process.
The analysis of KEGG pathway indicated that various transcripts related to cAMP
The analysis of KEGG pathway indicated that various transcripts related to cAMP
signaling pathway,
signaling pathway, Platelet
Platelet activation,
activation, and and Neurotrophin
Neurotrophin signaling
signaling pathwaypathway were
were over-
overrepresented in coelomocytes vs. intestine analysis (Table 2). Other
represented in coelomocytes vs. intestine analysis (Table 2). Other relevant pathways relevant pathways
representedinincoelomocytes
represented coelomocytesvs. vs. gonads
gonads correspond
correspond to to Fc
Fc gamma
gammaR-mediated
R-mediatedphagocyto-
phagocy-
sis, Platelet activation, and Pathogenic Escherichia coli infection
tosis, Platelet activation, and Pathogenic Escherichia coli infection (Table 2). (Table 2). In
Ingonad
gonadvs.vs.
intestine, the Purine metabolism was the main pathway overrepresented,
intestine, the Purine metabolism was the main pathway overrepresented, followed by followed by the
Spliceosome
the Spliceosome andandHuntington’s
Huntington’sdisease (Table
disease 3). We
(Table alsoalso
3). We found thatthat
found Biosynthesis
Biosynthesisof anti-
of
biotics, Gap junction, and Fatty acid degradation were highly represented
antibiotics, Gap junction, and Fatty acid degradation were highly represented in gonads vs. in gonads vs.
coelomocytes(Table
coelomocytes (Table3).
3).Regarding
Regardingthethemain
mainsignaling
signalingpathways
pathwaysoverrepresented
overrepresentedininintes-
intes-
tine(intestine
tine (intestinevs.
vs.coelomocytes),
coelomocytes),we wehighlight
highlightadherens
adherensjunction,
junction,Pathogenic
PathogenicEscherichia
Escherichia
coliinfection,
coli infection,
andand
ABC ABC transporters
transporters (Table
(Table 4). Finally,
4). Finally, the comparison
the comparison of intestine
of intestine vs. gonadvs.
gonad confirmed that pathways such ABC transporters and Pathogenic
confirmed that pathways such ABC transporters and Pathogenic Escherichia coli infection Escherichia coli
infection
were were
relevant relevant overrepresented
overrepresented pathways
pathways (Table 4). In(Table
summary,4). Inwesummary, we detected
detected 26,864 tran-
scripts were differentially expressed among the three tissues, which GO analysis revealed
Biology 2021, 10, 995 11 of 19

several processes and pathways that were expressed in common and, most importantly,
tissue specific.

Table 2. The Top-16 KEGG pathways enrichment of up-regulated transcripts of coelomocytes in the edible red sea urchin
(L. albus). In comparison to intestine, cAMP signaling pathway was the most overrepresented and Fc gamma R-mediated
phagocytosis compared to gonad.

Coelomocytes vs. Intestine

KEGG Term (ID) Number of Genes Percentage of Genes Fold Enrichment p Value
cAMP signaling pathway (4024) 13 2.9% 2.4 8.1 × 10−3
Platelet activation (4611) 10 2.2% 2.8 9.3 × 10−3
Neurotrophin signaling pathway (4722) 9 2.0% 2.7 1.7 × 10−2
Regulation of lipolysis in adipocytes (4923) 6 1.4% 3.8 1.8 × 10−2
Proteoglycans in cancer (5205) 12 2.7% 2.2 2.1 × 10−2
Focal adhesion (4510) 12 2.7% 2.1 2.6 × 10−2
Insulin secretion (4911) 7 1.6% 3.0 2.9 × 10−2
cGMP-PKG signaling pathway (4022) 10 2.2% 2.3 3.0 × 10−2
Coelomocytes vs. Gonad
KEGG Term (ID) Number of Genes Percentage of Genes Fold Enrichment p Value
Fc gamma R-mediated phagocytosis (4666) 15 2.5% 4.9 1.7 × 10−6
Platelet activation (4611) 16 2.6% 3.6 7.0 × 10−5
Pathogenic Escherichia coli infection (5130) 10 1.6% 5.4 7.9 × 10−5
cGMP-PKG signaling pathway (4022) 17 2.8% 2.9 1.9 × 10−4
Gap junction (4540) 12 1.9% 3.7 3.3 × 10−4
Proteoglycans in cancer (5205) 19 2.0% 2.6 3.5 × 10−4
Focal adhesion (4510) 19 3.1% 2.5 4.9 × 10−4
Bacterial invasion of epithelial cells (5100) 11 1.8% 3.8 4.9 × 10−4

Table 3. The Top-16 KEGG pathways enrichment of up-regulated transcripts of gonad in the edible red sea urchin (L. albus).
In comparison to intestine, purine metabolism was the most overrepresented and Biosynthesis of antibiotics compared
to coelomocytes.

Gonad vs. Intestine

KEGG Term (ID) Number of Genes Percentage of Genes Fold Enrichment p Value
Purine metabolism (230) 10 2.9% 2.4 6.9 × 10−3
Spliceosome (3040) 8 2.3% 2.8 1.4 × 10−2
Huntington’s disease (5016) 9 2.6% 2.7 3.2 × 10−2
p53 signaling pathway (4115) 5 1.5% 3.8 4.0 × 10−2
Focal adhesion (4510) 9 2.6% 2.1 4.6 × 10−2
Metabolic pathways (1100) 32 9.2% 3.0 5.9 × 10−2
Pathogenic Escherichia coli infection (5130) 4 1.5% 2.3 7.5 × 10−2
Gap junction (4540) 5 1.7% 2.4 9.0 × 10−2
Gonad vs. Coelomocytes
KEGG Term (ID) Number of Genes Percentage of Genes Fold Enrichment p Value
Biosynthesis of antibiotics (1130) 20 4.5% 3.4 4.5 × 10−6
Gap junction (4540) 10 2.3% 4.1 6.2 × 10−4
Fatty acid degradation (71) 7 1.6% 6.1 8.9 × 10−4
Valine, leucine, and isoleucine degradation
7 1.6% 5.4 1.6 × 10−3
(280)
Fatty acid metabolism (1212) 7 1.6% 5.3 1.8 × 10−3
Phagosome (4145) 12 2.7% 2.9 2.6 × 10−3
Metabolic pathways (1100) 47 10.6% 1.4 1.1 × 10−2
DNA replication (3030) 5 1.1% 5.1 1.6 × 10−2
Biology 2021, 10, 995 12 of 19

Table 4. The Top-16 KEGG pathways enrichment of up-regulated transcripts of intestine in the edible red sea urchin
(L. albus). In comparison to coelomocytes, adherens junction was the most overrepresented and ABC transporters compared
to gonad.

Intestine vs. Coelomocytes

KEGG Term (ID) Number of Genes Percentage of Genes Fold Enrichment p Value
Adherens junction (4520) 13 2.1% 5.0 7.7 × 10−6
Pathogenic Escherichia coli infection (5130) 11 1.7% 5.9 1.2 × 10−5
ABC transporters (2010) 10 1.6% 6.2 2.2 × 10−5
Lysosome (4142) 16 2.5% 3.6 2.9 × 10−5
Metabolic pathways (1100) 68 10.7% 1.5 1.8 × 10−4
Gap junction (4540) 11 1.7% 3.4 1.2 × 10−3
Chemical carcinogenesis (5204) 9 1.4% 3.1 8.2 × 10−3
cAMP signaling pathway (4024) 15 2.3% 2.1 1.3 × 10−2
Intestine vs. Gonad
KEGG Term (ID) Number of Genes Percentage of Genes Fold Enrichment p Value
ABC transporters (2010) 11 1.9% 7.9 7.1 × 10−7
Chemical carcinogenesis (5204) 12 2.2% 4.8 3.4 × 10−5
Pathogenic Escherichia coli infection (5130) 9 1.6% 5.6 1.6 × 10−4
Lysosome (4142) 12 2.2% 3.2 1.3 × 10−3
Gap junction (4540) 10 1.8% 3.6 1.6 × 10−3
Galactose metabolism (52) 6 1.1% 6.4 2.1 × 10−3
Starch and sucrose metabolism (500) 6 1.1% 5.8 3.3 × 10−3
Amino sugar and nucleotide sugar
7 1.3% 4.6 3.6 × 10−3
metabolism (520)

3.3. Transcriptomic Data Validation

We used real-time qPCR to validate RNA-seq and assay upregulated transcripts in
each tissue compared to the others. According to p values and fold enrichments indicated
in GO analysis, we selected four transcripts from coelomocytes: heat shock protein 70 kDa
1 A (HSP70), lysosomal trafficking regulator (LYST), B-cell lymphoma 2 (BCL2), and
ubiquitin A-52 residue ribosomal protein fusion product 1 (UBA52) (Figure 9a). Four
transcripts from gonad were identified: testis-specific serine/threonine-protein kinase 3
(TSSK3), centrin 2 (CETN2), cation channel sperm associated 3 (CATSPER3), and sperm
surface protein 17 (SPA17) (Figure 9b). Four transcripts from intestine were identified:
notch homolog 1 (NOTCH1), toll-like receptor 3 (TLR3), glutathione s-transferase theta
1 (GSTT1), and caspase 3 (CASP3) (Figure 9c). The transcript expression revealed a high
statistical correlation between RNA-seq and qPCR measuring (R2 = 0.7259, p < 0.001)
(Figure 9d). The data indicate that high correlation between RNA-seq and qPCR was
obtained, validating the selected transcripts.
 Biology
Biology 2021,
2021, 10,10,
995x 1314
ofof
1920

Figure
Figure The
9. 9. quantitative
The quantitative real-time PCR
real-time validation
PCR validationof differentially expressed
of differentially transcripts
expressed was was
transcripts highly correlated
highly to RNA-seq.
correlated to RNA-
(a)seq.
coelomocyte; (b) gonad; (c) intestine; and (d) statistical correlation analysis. Expression fold changes
(a) coelomocyte; (b) gonad; (c) intestine; and (d) statistical correlation analysis. Expression fold changes measured measured byby
RT-qPCR
RT-qPCR and RNA-seq
and RNA-seq areare
indicated
indicated byby
black and
black gray
and columns,
gray respectively.
columns, Abbreviations:
respectively. heat
Abbreviations: shock
heat protein
shock 70 kDa
protein 1
70 kDa
A 1(HSP70),
A (HSP70), lysosomal
lysosomal trafficking
trafficking regulator
regulator (LYST),
(LYST), B-cell
B-cell lymphoma
lymphoma 2 (BCL2),
2 (BCL2), ubiquitin
ubiquitin A-52
A-52 residue
residue ribosomal
ribosomal protein
protein
fusionproduct
fusion product 1 (UBA52),
1 (UBA52), testis-specific
testis-specific serine/threonine-proteinkinase
serine/threonine-protein kinase3 3(TSSK3),
(TSSK3),centrin
centrin2 2(CETN2),
(CETN2),cation
cation channel
channel
sperm
sperm associated
associated 3 (CATSPER3),
3 (CATSPER3), sperm
sperm surface
surface protein
protein 17 (SPA17),
17 (SPA17), notch homolog
notch homolog 1 (NOTCH1),
1 (NOTCH1), toll-like 3receptor
toll-like receptor (TLR3), 3
(TLR3), glutathione
glutathione s-transferase
s-transferase theta 1and
theta 1 (GSTT1), (GSTT1),
caspaseand caspase 3 *(CASP3).
3 (CASP3). p < 0.05, ** p < 0.01.

4.4.Discussion
Discussion
InInthis
thisstudy,
study,the
thetranscriptome
transcriptomeofofthe theedible
edibleredredsea
seaurchin
urchin(L.(L.albus)
albus)was wassequenced
sequenced
andannotated
and annotatedby byusing
usingNGS NGS technology.
technology. This This isis the
thefirst
firstreport
reporton onthetheRNARNAsequencing,
sequenc-
transcriptome
ing, transcriptome assembly,
assembly, andandfunctional
functional annotation
annotation from juvenile
from juvenileL. albus based
L. albus on three
based on
different tissues: intestines, gonads, and
three different tissues: intestines, gonads, and coelomocytes. coelomocytes. This transcriptome
transcriptome contains contains
185,239
185,239 transcripts
transcripts with
with similar
similar features
features to existing
to existing transcriptomes
transcriptomes of otherof other sea urchins.
sea urchins. For
For instance,
instance, the GC the GC content
content of L.(38.8%)
of L. albus albus (38.8%) has a similar
has a similar value invalue
relation in relation
to the seatourchin
the sea
urchin transcriptomes
transcriptomes of Evechinusof Evechinus chloroticus
chloroticus (39.0%) (39.0%) [32], Sterechinus
[32], Sterechinus neumayeri neumayeri
(38.6%)(38.6%)
[33],
[33], Strongylocentrotus
Strongylocentrotus intermediusintermedius (39.6%)
(39.6%) [34], and [34], and Mesocentrotus
Mesocentrotus nudus (39.9%) nudus[35].(39.9%)
However, [35].
However,
our our results
results exhibit a GCexhibit
contenta GC content
slightly lower slightly lower than
than values values
reported forreported
the testisfor the testis
transcrip-
transcriptome
tome of L. albus
of L. albus (40.4%) [6]. (40.4%)
This small[6].difference
This smallmay difference may be
be attributed toattributed
the tissue usedto theintissue
the
used(mature
study in the study (mature
gonad), as wellgonad), as well as the
as the sequencing sequencing
technology (Rochetechnology
454 GS-FLX (Roche 454 GS-
Titanium).
This
FLXgroup obtained
Titanium). This1062,716 raw reads
group obtained with araw
1062,716 mean readslength
withof 309.8 length
a mean bp, generating
of 309.8 bp,a
reference
generating transcriptome of 42,530 transcripts
a reference transcriptome of 42,530 with an N50 with
transcripts of 645anbp N50[6],ofin645
contrast
bp [6],toinour
con-
study, which
trast to presents
our study, 91,119,300
which presents paired end reads
91,119,300 with
paired enda mean
reads length
with a of mean250 length
bp, generating
of 250 bp,
a generating
reference transcriptome of 185,239 transcripts
a reference transcriptome of 185,239with an N50with
transcripts of 1769
an N50 bp. of 1769 bp.
AAdetailed
detailed analysis
analysis of tissue
of tissueexpression
expression reveal that that
reveal 26,864 transcripts
26,864 are differentially
transcripts are differen-
expressed among the three tissues. Coelomocytes are the cells
tially expressed among the three tissues. Coelomocytes are the cells responsible responsible for immunity for in
im-
sea urchins,
munity of which
in sea 50%oftowhich
urchins, 70% are
50% motile
to 70% cells
arewith a high
motile cellsenergy
with arequirement
high energyand are
require-
considered
ment and equivalent
are consideredto human macrophages,
equivalent to humanwhich are cellular
macrophages, immune
which system compo-
are cellular immune
nents
systemwith high catabolic
components withactivity and are part
high catabolic of the
activity andinnate immune
are part of theresponse involvedre-
innate immune
insponse
pathogen digestion
involved and autophagy
in pathogen digestion [36,37]. Besides coelomocytes
and autophagy [36,37]. Besides are coelomocytes
predominantly are
Biology 2021, 10, 995 14 of 19

in the coelomic fluid, they also function as wandering cells and infiltrate all tissues [36].
Consequently, the differentially expressed transcripts in the coelomocytes are mainly as-
sociated with biological processes, such as positive regulation of apoptotic process and
intracellular signal transduction, and KEGG pathways associated with Platelet activation,
Fc gamma R-mediated phagocytosis, and Pathogenic Escherichia coli infection. Studies in
sea urchin coelomocytes transcriptomes describe similar observations. In S. intermedius, the
expression of 546 unique transcripts in coelomocytes is associated with lysozyme, lectin,
pattern recognition receptors (PRRs), and the complement system [38]. In a related study
carried out in coelomocytes of Arbacia lixula, the expression of transcripts is associated to
lipid metabolism and the immune response [39]. In addition, an RNA-seq analysis in coelo-
mocytes reveal key functions of NOD-like receptor pathway and phagosomes in spotting
diseased S. intermedius [40]. Recently, the immune response of L. albus coelomocytes by
poly I:C, bacterial lipopolysaccharides (LPS), and temperature reveal a dynamic expression
of TLR genes (tlr3 including), as well as strongylocin-1 and strongylocin-2 [18].
Among the transcripts identified with a high expression in coelomocytes and validated
by RT-qPCR stand out the heat shock protein 70 kDa 1 A (HSP70), the lysosomal trafficking
regulator (LYST), the B-cell lymphoma 2 (BCL2) and the ubiquitin A-52 residue ribosomal
protein fusion product 1 (UBA52). HSP70 is a chaperone protein responsible for protein
folding to protect cells against stressors or presenting antigens for immune response [41].
Interestingly, a recent study has shown that LPS can induce a stress response by increasing
the protein levels of HSP70 in Paracentrotus lividus coelomocytes, suggesting a relevant role
in the sea urchin immune response [42]. LYST plays a role in the transport of materials into
structures called lysosomes, acting as recycling centers within cells [43]. Although there
are no reports of the importance of this gene in sea urchin coelomocytes, in mammalian
macrophages has been linked as a key regulator of membrane trafficking to inflammatory
responses mediated by TLRs [44]. BCL2 is a member of protein regulators for cell death,
through inhibition of apoptosis [45]. The participation of BCL2 as an important mediator
of the immune response in marine organisms has recently been described in Apostichopus
japonicus challenged with Vibrio splendidus [46]. The UBA52 gene encodes to 60S ribosomal
protein L40 (RPL40) and, together with ubiquitin, has a main function of targeting proteins
for degradation by the 26S proteosome. Additionally, UBA52 can regulate gene expression,
chromatin structure, and the stress response [47]. Although there are no reports of the
relevance of UBA52 in the immune response of sea urchin coelomocytes, its participation
in the immune response of higher vertebrates has been described [48]. These observations
suggest a permanent activity of protein catabolism in sea urchin coelomocytes, as sentinel
organisms of the immune response.
In sea urchins, the gonads are considered a dual organ since they regulate reproduction
and nutrient storage [49]. Both functions are carried out by two specific cell populations, so-
matic and germ cells, the latter also called nutritive phagocytes [50]. Several transcriptomic
analyses have been described in gonads of edible sea urchin species, obtaining results in
agreement with ours. Particularly, a de novo assembly of M. nudus gonad transcriptome, key
genes associated with biological processes such lipid metabolism, and the biosynthesis of
polyunsaturated fatty acids have been identified [35]. Furthermore, a comparative analysis
between reproductive tissues of this specie, reveal upregulated GO categories related to
energy generation in the testis and negative regulation of nucleotide metabolism in the
ovary [39]. Additionally, a transcriptomic analysis indicates important differences in the
expression of genes of M. nudus related to the biosynthesis of polyunsaturated fatty acids
and metabolism in relation to high-quality gonads, similar to our observations [34,51]. In
addition, a gonadal transcriptomic of M. nudus present candidate sex-related genes that
could be involved in significant roles in spermatogenesis, oogenesis, and germ cell devel-
opment [52]. More recently, integrated analyses of metabolomic and transcriptomic reveal
key genes for metabolism and eicosapentaenoic acid biosynthesis in the sea urchin S. inter-
medius, identifying six accumulated metabolites and several differentially expressed genes
associated with polyunsaturated fatty acids in the testis compared with the ovary [53].
Biology 2021, 10, 995 15 of 19

The testis-specific serine/threonine-protein kinase 3 (TSSK3), centrin 2 (CETN2),

cation channel sperm associated 3 (CATSPER3) and sperm surface protein 17 (SPA17) were
upregulated specifically in gonad. TSSK3 is a protein involved in the development and
maturation of male germ cell [54]. TSSK3 role as a key regulator in spermiogenesis of
the Scallop Argopecten irradians has recently been described [55]. CETN2 is a member of
calcium-binding proteins and is also a component of the centrosome [56]. However, there
are no previous reports of its importance in the gonadal function of marine echinoderms.
CATSPER3 is a sperm-specific ion channel that plays a central role for the successful fertil-
ization, which involves sperm hyperactivation, acrosome reaction, and chemotaxis towards
the egg [57]. In mammals, it has been determined that CATSPER3 is exclusively expressed
in the testis, which should be validated in future analyzes in adult sea urchins [58]. SPA17
encodes a sperm surface zona pellucida binding protein. In mammals, SPA17 may be impli-
cated in fertilization through zona pellucida attaching of the oocyte [59]. Members of this
group of genes have shown an important participation in the maturation and fertilization
of echinoderm eggs [60]. However, it is difficult to speculate if differentially expressed
transcripts of gonads are relevant in the function of the testis or ovary, because the present
study was performed with sexually immature juveniles L. albus.
Finally, the intestine is the organ responsible for the digestion of food [61]. In echin-
oderms, it also has a key role in the immune response, where it has been described that
it presents a high number of sequences related to this process [62,63]. Consistently, we
observe an overrepresentation of GO and KEGG categories similar to that observed in coelo-
mocytes, showing high enrichment in the sequences related to SSCR and Cargo Receiver.
This is consistent with the importance of the intestine in the immune response; therefore, it
has been postulated that the control of normal flora or pathogen-associated gut response is
responsible for the increase in the expression of the pattern recognition receptors NLR and
TLR [64]. Similarly, we also detect an overrepresentation of lysosomal pathways, which
was also has been found in the transcriptomic analysis in S. intermedius [38].
The notch homolog 1 (NOTCH1), toll-like receptor 3 (TLR3), glutathione s-transferase
theta 1 (GSTT1) and caspase 3 (CASP3) were upregulated specifically in the intestine.
NOTCH1 is a Type 1 transmembrane protein involves in numerous processes of devel-
opment by regulating cellular fate. The Notch signaling pathway is highly conserved
among species and regulates interactions between adjacent cells [65]. Interestingly, the
importance of intestinal epithelial NOTCH1 as a protector for the development of colorectal
adenocarcinoma in a murine model has been described [66]. TLR3 is a member of pattern
recognition receptors that recognizes specific molecules of pathogens during innate im-
mune response [67]. Recently our group described the upregulation of TLR3 expression in
the coelomocyte response to bacterial LPS, poly I:C and temperature in L. albus [18]. GSTT1
is an enzyme that catalyzes the coupling of reduced glutathione to diverse hydrophobic
and electrophilic compounds [68]. Remarkably, the glutathione S-transferase activity in the
anterior portion of P. lividus intestine has been employed as biomarker of environmental
contamination, revealing an important role in the homeostasis of this tissue [69]. Finally,
CASP3 gene encodes for a cysteine-aspartic acid protease (caspase) that has a crucial role
in the execution-phase of apoptosis [70]. This protein has been described as an essential
element if the innate immunity of sea urchins [64]. These observations suggest the presence
of coelomocytes in the sea urchin intestinal tract, as a fundamental part of the crosstalk
between intestinal microbiota and the inflammatory response.

5. Conclusions
This is the first evidence on the RNA sequencing, de novo assembly, and functional
annotation of the edible red sea urchin (L. albus) transcriptome, as well as the differential
expression in the gonads, intestines, and coelomocytes. The de novo assembly produced
185,239 transcripts, creating a reference transcriptome with an N50 of 1769 bp and 38.8% GC
content. Gene ontology analysis of transcripts revealed notable differences in the expression
profiles between gonads, intestines, and coelomocytes, allowing the detection of transcripts
Biology 2021, 10, 995 16 of 19

associated with specific biological processes. In coelomocytes, DETs were mostly associated
to positive regulation of apoptotic process and intracellular signal transduction. In the
gonad, DETs were associated to DNA repair and microtubule-based process. In Intestine,
DETs were associated to transmembrane transport and microtubule-based process. The
dataset generated in this work contribute to enrich the molecular resources of L. albus,
improvement futures biological studies of this species. The acquired information is also
relevant to discovering novel candidate genes that could be employed to estimate the
physiological condition of edible red sea urchins under aquaculture rearing.

Supplementary Materials: The following are available online at https://www.mdpi.com/article/10

.3390/biology10100995/s1, Table S1: Primer sequences for qPCR assays, including amplicon size
and PCR efficiencies for the assessed genes, Table S2: The complete list of differentially expressed
transcripts in the edible sea urchin (L. albus), Figure S1: Length distribution of assembled transcripts,
Figure S2: The heat map analysis of transcript expression across tissues in L. albus was highly
correlated among replicates, Figure S3: The Top-16 Gene Ontology molecular function (MF) and
cellular component (CC) enrichment of up-regulated transcripts from L. albus coelomocytes, Figure S4:
The Top-16 Gene Ontology molecular function (MF) and cellular component (CC) enrichment of up-
regulated transcripts from L. albus gonad, Figure S5: The Top-16 Gene Ontology molecular function
(MF) and cellular component (CC) enrichment of up-regulated transcripts from L. albus intestine.
Author Contributions: Conceptualization, J.M.E. and J.A.V.; methodology, P.A., J.M.E. and J.A.V.;
software, P.A.; validation, J.M.E., C.M., A.M. and J.A.V.; formal analysis, P.A. and M.B.-M.; inves-
tigation, P.A. and J.A.V.; resources, J.M.E., A.M. and J.A.V.; data curation, P.A.; writing—original
draft preparation, R.Z. and J.A.V.; writing—review and editing, R.Z. and J.A.V.; visualization, J.M.E.,
C.M., A.M. and J.A.V.; supervision, J.M.E., C.M., A.M. and J.A.V.; project administration, A.M. and
J.A.V.; funding acquisition, J.A.V. All authors have read and agreed to the published version of
the manuscript.
Funding: This research was funded by Fondo Nacional de Desarrollo Científico y Tecnológico
(FONDECYT), grant number 1201498 (to Juan Antonio Valdés), Fondo de Financiamiento de Centros
de Investigación en Áreas Prioritarias (FONDAP) grant INCAR 15110027 and Fondo de Innovación
para la Competitividad Regional (FIC), grant BIP 30480912-0.
Institutional Review Board Statement: The study was conducted according to the guidelines of
the Declaration of Helsinki and approved by the Ethics Committee of Universidad Andres Bello
and the National Commission for Scientific and Technological Research (CONICYT) of the Chilean
government (Nº 010P/2021).
Informed Consent Statement: Not applicable.
Data Availability Statement: The raw read sequences obtained from sequencing were deposited
in the Sequence Read Archive (SRA) under BioProject PRJNA475570, with accession number
SRP150640. The datasets generated and/or analyzed during the current study are not publicly
available due to privacy or ethical restrictions but are available from the corresponding author
upon reasonable request.
Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design
of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, or
in the decision to publish the results.

References
1. Oyarzún, S.T.; Marín, S.L.; Valladares, C.; Iriarte, J.L. Reproductive Cycle of Loxechinus albus (Echinodermata: Echinoidea) in Two
Areas of the Magellan Region (53oS, 70-72oW), Chile. Sci. Mar. 1999, 63, 439–449. [CrossRef]
2. Vásquez, J.A.; Donoso, G.A. Loxechinus albus. In Developments in Aquaculture and Fisheries Science; Elsevier: Amsterdam, The
Netherlands, 2013; Volume 38, pp. 285–296, ISBN 978-0-12-396491-5.
3. Vásquez, J.A. Loxechinus albus. In Developments in Aquaculture and Fisheries Science; Elsevier: Amsterdam, The Netherlands, 2020;
Volume 43, pp. 431–445, ISBN 978-0-12-819570-3.
4. Olave, S.; Bustos, E.; Lawrence, J.M.; Carcamo, P. The Effect of Size and Diet on Gonad Production by the Chilean Sea Urchin
Loxechinus albus. JWAS 2001, 32, 210–214. [CrossRef]
Biology 2021, 10, 995 17 of 19

5. Cea, G.; Gaitán-Espitia, J.D.; Cárdenas, L. Complete Mitogenome of the Edible Sea Urchin Loxechinus albus: Genetic Structure and
Comparative Genomics within Echinozoa. Mol. Biol. Rep. 2015, 42, 1081–1089. [CrossRef]
6. Gaitán-Espitia, J.D.; Sánchez, R.; Bruning, P.; Cárdenas, L. Functional Insights into the Testis Transcriptome of the Edible Sea
Urchin Loxechinus albus. Sci. Rep. 2016, 6, 36516. [CrossRef]
7. Chandhini, S.; Rejish Kumar, V.J. Transcriptomics in Aquaculture: Current Status and Applications. Rev. Aquac. 2018, 11,
1379–1397. [CrossRef]
8. Inbakandan, D. Transcriptomics in Aquaculture. In Encyclopedia of Marine Biotechnology; Kim, S., Ed.; Wiley: Hoboken, NJ, USA,
2020; pp. 1919–1936, ISBN 978-1-119-14377-2.
9. Jung, G.; Lee, Y.-H. Complete Mitochondrial Genome of Chilean Sea Urchin: Loxechinus albus (Camarodonta, Parechinidae).
Mitochondrial DNA 2015, 26, 883–884. [CrossRef]
10. Oulas, A.; Pavloudi, C.; Polymenakou, P.; Pavlopoulos, G.A.; Papanikolaou, N.; Kotoulas, G.; Arvanitidis, C.; Iliopoulos, L.
Metagenomics: Tools and Insights for Analyzing Next-Generation Sequencing Data Derived from Biodiversity Studies. Bioinform.
Biol. Insights 2015, 9, 75–88. [CrossRef] [PubMed]
11. Vergara-Amado, J.; Silva, A.X.; Manzi, C.; Nespolo, R.F.; Cárdenas, L. Differential Expression of Stress Candidate Genes for
Thermal Tolerance in the Sea Urchin Loxechinus albus. J. Therm. Biol. 2017, 68, 104–109. [CrossRef] [PubMed]
12. Wong, J.M.; Johnson, K.M.; Kelly, M.W.; Hofmann, G.E. Transcriptomics Reveal Transgenerational Effects in Purple Sea Urchin
Embryos: Adult Acclimation to Upwelling Conditions Alters the Response of Their Progeny to Differential pCO2 Levels. Mol.
Ecol. 2018, 27, 1120–1137. [CrossRef]
13. Zhan, Y.; Li, J.; Sun, J.; Zhang, W.; Li, Y.; Cui, D.; Hu, W.; Chang, Y. The Impact of Chronic Heat Stress on the Growth, Survival,
Feeding, and Differential Gene Expression in the Sea Urchin Strongylocentrotus Intermedius. Front. Genet. 2019, 10, 301.
[CrossRef]
14. Malanga, G.; Perez, A.; Calvo, J.; Puntarulo, S. The Effect of Seasonality on Oxidative Metabolism in the Sea Urchin Loxechinus
albus. Mar. Biol. 2009, 156, 763–770. [CrossRef]
15. Flores, L.; Ernst, B.; Parma, A.M. Growth Pattern of the Sea Urchin, Loxechinus albus (Molina, 1782) in Southern Chile: Evaluation
of Growth Models. Mar. Biol. 2010, 157, 967–977. [CrossRef]
16. Cárcamo, P.F. Effects of Food Type and Feeding Frequency on the Performance of Early Juveniles of the Sea Urchin Loxechinus
albus (Echinodermata: Echinoidea): Implications for Aquaculture and Restocking. Aquaculture 2015, 436, 172–178. [CrossRef]
17. Dupré, E.; Carvajal, J. Cryopreservation of Embryos and Larvae of the Edible Sea Urchin Loxechinus albus (Molina, 1782).
Cryobiology. 2019, 86, 84–88. [CrossRef]
18. Dettleff, P.; Villagra, M.; González, J.; Fuentes, M.; Estrada, J.M.; Valenzuela, C.; Molina, A.; Valdés, J.A. Effect of Bacterial LPS,
Poly I:C and Temperature on the Immune Response of Coelomocytes in Short Term Cultures of Red Sea Urchin (Loxechinus albus).
Fish Shellfish Immunol. 2020, S1050464820306513. [CrossRef]
19. Arévalo, A.; Bruning, P.; Sanchez, R.; Cárdenas, L. Development and Characterization of EST-Microsatellites for the Edible Sea
Urchin Loxechinus albus Using next Generation Sequencing. Conserv. Genet. Resour. 2014, 6, 433–435. [CrossRef]
20. Dodt, M.; Roehr, J.; Ahmed, R.; Dieterich, C. FLEXBAR—Flexible Barcode and Adapter Processing for Next-Generation Sequenc-
ing Platforms. Biology. 2012, 1, 895–905. [CrossRef] [PubMed]
21. Grabherr, M.G.; Haas, B.J.; Yassour, M.; Levin, J.Z.; Thompson, D.A.; Amit, I.; Adiconis, X.; Fan, L.; Raychowdhury, R.; Zeng, Q.;
et al. Full-Length Transcriptome Assembly from RNA-Seq Data without a Reference Genome. Nat. Biotechnol. 2011, 29, 644–652.
[CrossRef] [PubMed]
22. Davidson, N.M.; Oshlack, A. Corset: Enabling Differential Gene Expression Analysis for de Novoassembled Transcriptomes. Genome
Biol. 2014, 15, 410. [CrossRef] [PubMed]
23. Simão, F.A.; Waterhouse, R.M.; Ioannidis, P.; Kriventseva, E.V.; Zdobnov, E.M. BUSCO: Assessing Genome Assembly and
Annotation Completeness with Single-Copy Orthologs. Bioinformatics. 2015, 31, 3210–3212. [CrossRef]
24. Conesa, A.; Gotz, S.; Garcia-Gomez, J.M.; Terol, J.; Talon, M.; Robles, M. Blast2GO: A Universal Tool for Annotation, Visualization
and Analysis in Functional Genomics Research. Bioinformatics. 2005, 21, 3674–3676. [CrossRef]
25. Ye, J.; Zhang, Y.; Cui, H.; Liu, J.; Wu, Y.; Cheng, Y.; Xu, H.; Huang, X.; Li, S.; Zhou, A.; et al. WEGO 2.0: A Web Tool for Analyzing
and Plotting GO Annotations, 2018 Update. Nucleic Acids Res. 2018, 46, W71–W75. [CrossRef]
26. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and Memory-Efficient Alignment of Short DNA Sequences to the
Human Genome. Genome Biol. 2009, 10, R25. [CrossRef] [PubMed]
27. Li, B.; Dewey, C.N. RSEM: Accurate Transcript Quantification from RNA-Seq Data with or without a Reference Genome. BMC
Bioinform. 2011, 12, 323. [CrossRef]
28. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital
Gene Expression Data. Bioinformatics 2010, 26, 139–140. [CrossRef] [PubMed]
29. Huang, D.W.; Sherman, B.T.; Lempicki, R.A. Systematic and Integrative Analysis of Large Gene Lists Using DAVID Bioinformatics
Resources. Nat. Protoc. 2009, 4, 44–57. [CrossRef] [PubMed]
30. Bustin, S.A.; Benes, V.; Garson, J.A.; Hellemans, J.; Huggett, J.; Kubista, M.; Mueller, R.; Nolan, T.; Pfaffl, M.W.; Shipley, G.L.; et al.
The MIQE Guidelines: Minimum Information for Publication of Quantitative Real-Time PCR Experiments. Clin. Chem. 2009, 55,
611–622. [CrossRef] [PubMed]
31. Simon, P. Q-Gene: Processing Quantitative Real-Time RT-PCR Data. Bioinformatics 2003, 19, 1439–1440. [CrossRef] [PubMed]
Biology 2021, 10, 995 18 of 19

32. Gillard, G.B.; Garama, D.J.; Brown, C.M. The Transcriptome of the NZ Endemic Sea Urchin Kina (Evechinus chloroticus). BMC
Genom. 2014, 15, 45. [CrossRef]
33. Dilly, G.F.; Gaitán-Espitia, J.D.; Hofmann, G.E. Characterization of the Antarctic Sea Urchin (Sterechinus neumayeri) Transcriptome
and Mitogenome: A Molecular Resource for Phylogenetics, Ecophysiology and Global Change Biology. Mol. Ecol. Resour. 2015,
15, 425–436. [CrossRef]
34. Wang, H.; Ding, J.; Ding, S.; Chang, Y. Transcriptome Analysis to Characterize the Genes Related to Gonad Growth and Fatty
Acid Metabolism in the Sea Urchin Strongylocentrotus intermedius. Genes Genom. 2019, 41, 1397–1415. [CrossRef] [PubMed]
35. Wei, Z.; Liu, X.; Zhou, Z.; Xu, J. De Novo Transcriptomic Analysis of Gonad of Strongylocentrotus nudus and Gene Discovery for
Biosynthesis of Polyunsaturated Fatty Acids. Genes Genom. 2019, 41, 583–597. [CrossRef] [PubMed]
36. Da Silva, J.R.M.C. Immunology in sea urchins. In Developments in Aquaculture and Fisheries Science; Elsevier: Amsterdam, The
Netherlands, 2020; Volume 43, pp. 227–236. ISBN 978-0-12-819570-3.
37. Smith, L.C.; Hawley, T.S.; Henson, J.H.; Majeske, A.J.; Oren, M.; Rosental, B. Methods for Collection, Handling, and Analysis of
Sea Urchin Coelomocytes. Methods Cell. Biol. 2019, 150, 357–389. [CrossRef] [PubMed]
38. Chen, Y.; Chang, Y.; Wang, X.; Qiu, X.; Liu, Y. De Novo Assembly and Analysis of Tissue-Specific Transcriptomes Revealed
the Tissue-Specific Genes and Profile of Immunity from Strongylocentrotus intermedius. Fish Shellfish Immunol. 2015, 46, 723–736.
[CrossRef]
39. Pérez-Portela, R.; Turon, X.; Riesgo, A. Characterization of the Transcriptome and Gene Expression of Four Different Tissues in
the Ecologically Relevant Sea Urchin Arbacia lixula Using RNA-Seq. Mol. Ecol. Resour. 2016, 16, 794–808. [CrossRef] [PubMed]
40. Zhang, W.; Lv, Z.; Li, C.; Sun, Y.; Jiang, H.; Zhao, M.; Zhao, X.; Shao, Y.; Chang, Y. Transcriptome Profiling Reveals Key Roles of
Phagosome and NOD-like Receptor Pathway in Spotting Diseased Strongylocentrotus intermedius. Fish Shellfish Immunol. 2019, 84,
521–531. [CrossRef]
41. Lubkowska, A.; Pluta, W.; Strońska, A.; Lalko, A. Role of Heat Shock Proteins (HSP70 and HSP90) in Viral Infection. IJMS 2021,
22, 9366. [CrossRef]
42. Chiaramonte, M.; Inguglia, L.; Vazzana, M.; Deidun, A.; Arizza, V. Stress and Immune Response to Bacterial LPS in the Sea
Urchin Paracentrotus lividus (Lamarck, 1816). Fish Shellfish Immunol. 2019, 92, 384–394. [CrossRef]
43. Sepulveda, F.E.; Burgess, A.; Heiligenstein, X.; Goudin, N.; Ménager, M.M.; Romao, M.; Côte, M.; Mahlaoui, N.; Fischer, A.;
Raposo, G.; et al. LYST Controls the Biogenesis of the Endosomal Compartment Required for Secretory Lysosome Function: LYST
Controls the Biogenesis of Secretory Lysosome. Traffic 2015, 16, 191–203. [CrossRef]
44. Westphal, A.; Cheng, W.; Yu, J.; Grassl, G.; Krautkrämer, M.; Holst, O.; Föger, N.; Lee, K.-H. Lysosomal Trafficking Regulator Lyst
Links Membrane Trafficking to Toll-like Receptor–Mediated Inflammatory Responses. J. Exp. Med. 2017, 214, 227–244. [CrossRef]
45. Banjara, S.; Suraweera, C.D.; Hinds, M.G.; Kvansakul, M. The Bcl-2 Family: Ancient Origins, Conserved Structures, and Divergent
Mechanisms. Biomolecules 2020, 10, 128. [CrossRef]
46. Guo, M.; Chen, K.; Lv, Z.; Shao, Y.; Zhang, W.; Zhao, X.; Li, C. Bcl-2 Mediates Coelomocytes Apoptosis by Suppressing Cytochrome
c Release in Vibrio splendidus Challenged Apostichopus Japonicus. Dev. Comp. Immunol. 2020, 103, 103533. [CrossRef] [PubMed]
47. Kobayashi, M.; Oshima, S.; Maeyashiki, C.; Nibe, Y.; Otsubo, K.; Matsuzawa, Y.; Nemoto, Y.; Nagaishi, T.; Okamoto, R.; Tsuchiya,
K.; et al. The Ubiquitin Hybrid Gene UBA52 Regulates Ubiquitination of Ribosome and Sustains Embryonic Development. Sci.
Rep. 2016, 6, 36780. [CrossRef] [PubMed]
48. Ebner, P.; Versteeg, G.A.; Ikeda, F. Ubiquitin Enzymes in the Regulation of Immune Responses. Crit. Rev. Biochem. Mol. Biol. 2017,
52, 425–460. [CrossRef] [PubMed]
49. Walker, C.W.; Unuma, T.; Lesser, M.P. Chapter 2 Gametogenesis and reproduction of sea urchins. In Edible Sea Urchins: Biology
and Ecology; Lawrence, J.M., Ed.; Developments in Aquaculture and Fisheries Science; Elsevier: Amsterdam, The Netherlands, 2007;
Volume 37, pp. 11–33.
50. Walker, C.W.; Harrington, L.M.; Lesser, M.P.; Fagerberg, W.R. Nutritive Phagocyte Incubation Chambers Provide a Structural and
Nutritive Microenvironment for Germ Cells of Strongylocentrotus droebachiensis, the Green Sea Urchin. Biol. Bull. 2005, 209, 31–48.
[CrossRef]
51. Jia, Z.; Wang, Q.; Wu, K.; Wei, Z.; Zhou, Z.; Liu, X. De Novo Transcriptome Sequencing and Comparative Analysis to Discover
Genes Involved in Ovarian Maturity in Strongylocentrotus nudus. Comp. Biochem. Physiol. D Genom. Proteom. 2017, 23, 27–38.
[CrossRef] [PubMed]
52. Sun, Z.-H.; Zhang, J.; Zhang, W.-J.; Chang, Y.-Q. Gonadal Transcriptomic Analysis and Identification of Candidate Sex-Related
Genes in Mesocentrotus nudus. Gene 2019, 698, 72–81. [CrossRef]
53. Wang, H.; Ding, J.; Ding, S.; Chang, Y. Integrated Metabolomic and Transcriptomic Analyses Identify Critical Genes in Eicos-
apentaenoic Acid Biosynthesis and Metabolism in the Sea Urchin Strongylocentrotus intermedius. Sci. Rep. 2020, 10, 1697.
[CrossRef]
54. Salicioni, A.M.; Gervasi, M.G.; Sosnik, J.; Tourzani, D.A.; Nayyab, S.; Caraballo, D.A.; Visconti, P.E. Testis-Specific Serine Kinase
Protein Family in Male Fertility and as Targets for Non-Hormonal Male Contraception. Biol. Reprod. 2020, 103, 264–274. [CrossRef]
55. Xue, X.; Zhang, L.; Li, Y.; Wei, H.; Wu, S.; Liu, T.; Liu, L.; Xing, Q.; Wang, S.; Bao, Z. Expression of the Testis-Specific
Serine/Threonine Kinases Suggests Their Role in Spermiogenesis of Bay Scallop Argopecten irradians. Front. Physiol. 2021,
12, 657559. [CrossRef]
Biology 2021, 10, 995 19 of 19

56. Mullee, L.I.; Morrison, C.G. Centrosomes in the DNA Damage Response—The Hub Outside the Centre. Chromosome Res. 2016, 24,
35–51. [CrossRef]
57. Singh, A.P.; Rajender, S. CatSper Channel, Sperm Function and Male Fertility. Reprod. Biomed. Online 2015, 30, 28–38. [CrossRef]
58. Jin, J.-L.; O’Doherty, A.M.; Wang, S.; Zheng, H.; Sanders, K.M.; Yan, W. Catsper3 and Catsper4 Encode Two Cation Channel-Like
Proteins Exclusively Expressed in the Testis. Biol. Reprod. 2005, 73, 1235–1242. [CrossRef]
59. Chiriva-Internati, M. Sperm Protein 17: Clinical Relevance of a Cancer/Testis Antigen, from Contraception to Cancer Immunother-
apy, and Beyond. Int. Rev. Immunol. 2011, 30, 138–149. [CrossRef]
60. Santella, L.; Limatola, N.; Vasilev, F.; Chun, J.T. Maturation and Fertilization of Echinoderm Eggs: Role of Actin Cytoskeleton
Dynamics. Biochem. Biophys. Res. Commun. 2018, 506, 361–371. [CrossRef] [PubMed]
61. Holland, N.D. Digestive system in regular sea urchins. In Developments in Aquaculture and Fisheries Science; Elsevier: Amsterdam,
The Netherlands, 2020; Volume 43, pp. 147–163. ISBN 978-0-12-819570-3.
62. Ramírez-Gómez, F.; Ortíz-Pineda, P.A.; Rojas-Cartagena, C.; Suárez-Castillo, E.C.; García-Arrarás, J.E.; García-Ararrás, J.E.
Immune-Related Genes Associated with Intestinal Tissue in the Sea Cucumber Holothuria glaberrima. Immunogenetics 2008, 60,
57–71. [CrossRef]
63. Yang, A.-F.; Zhou, Z.-C.; He, C.-B.; Hu, J.-J.; Chen, Z.; Gao, X.-G.; Dong, Y.; Jiang, B.; Liu, W.-D.; Guan, X.-Y.; et al. Analysis of
Expressed Sequence Tags from Body Wall, Intestine and Respiratory Tree of Sea Cucumber (Apostichopus japonicus). Aquaculture
2009, 296, 193–199. [CrossRef]
64. Hibino, T.; Loza-Coll, M.; Messier, C.; Majeske, A.J.; Cohen, A.H.; Terwilliger, D.P.; Buckley, K.M.; Brockton, V.; Nair, S.V.; Berney,
K.; et al. The Immune Gene Repertoire Encoded in the Purple Sea Urchin Genome. Dev. Biol. 2006, 300, 349–365. [CrossRef]
[PubMed]
65. Garis, M.; Garrett-Sinha, L.A. Notch Signaling in B Cell Immune Responses. Front. Immunol. 2021, 11, 609324. [CrossRef]
[PubMed]
66. Dunkin, D.; Iuga, A.C.; Mimouna, S.; Harris, C.L.; Haure-Mirande, J.-V.; Bozec, D.; Yeretssian, G.; Dahan, S. Intestinal Epithelial
Notch-1 Protects from Colorectal Mucinous Adenocarcinoma. Oncotarget 2018, 9, 33536–33548. [CrossRef] [PubMed]
67. Habib, Y.J.; Zhang, Z. The Involvement of Crustaceans Toll-like Receptors in Pathogen Recognition. Fish Shellfish Immunol. 2020,
102, 169–176. [CrossRef] [PubMed]
68. Vaish, S.; Gupta, D.; Mehrotra, R.; Mehrotra, S.; Basantani, M.K. Glutathione S-Transferase: A Versatile Protein Family. 3 Biotech.
2020, 10, 321. [CrossRef] [PubMed]
69. Cunha, I.; García, L.M.; Guilhermino, L. Sea-Urchin (Paracentrotus lividus) Glutathione S-Transferases and Cholinesterase Activities
as Biomarkers of Environmental Contamination. J. Environ. Monit. 2005, 7, 288–294. [CrossRef] [PubMed]
70. Asadi, M.; Taghizadeh, S.; Kaviani, E.; Vakili, O.; Taheri-Anganeh, M.; Tahamtan, M.; Savardashtaki, A. Caspase-3: Structure,
Function, and Biotechnological Aspects. Biotechnol. Appl. Biochem. 2021, 1–13. [CrossRef]