Healthy European sea bass (Dicentrarchus labrax), with an average weight of 50 g, were provided by a commercial fish farm in the north of Spain (Sonríonansa S.L., Pesués, Cantabria, Spain) (22, 23). Prior to the experiments, fish were acclimated for 30 days to the fish holding facilities of the Instituto de Ciências Biomédicas Abel Salazar (ICBAS), Porto. Fish were kept in 110 liters recirculating sea water (28 ± 1‰ salinity) tanks at 22 ± 1°C, with a 13/11-hour light/dark cycle and fed daily ad libitum with commercial fish feed with an iron content of approximately 200 mg iron/kg feed. Before each treatment, fish were anaesthetized with ethylene glycol monophenyl ether (2-phenoxyethanol, 3 mL/10L, Merck, Algés, Portugal). All animal experiments were carried out in strict compliance with national and international animal use ethics guidelines (including ARRIVE guidelines), approved by the animal welfare and ethic committees of ICBAS (permits P293/2019/ORBEA and P375/2020/ORBEA) and DGAV – Direcção-Geral da Alimentação e Veterinária, and conducted by experienced and trained FELASA Function A+B+D investigators.

Pairs of oligonucleotide PCR primers were designed according to conserved regions of erfe mRNA sequences from other fish and mammalian species, available in the National Center for Biotechnology Information nucleotide database (http://www.ncbi.nlm.nih.gov/) and Ensembl (http://www.ensembl.org/) and cDNA preparations from spleen and kidney were used in PCR amplifications (22, 23). PCR products were run on 1.2% agarose gels, relevant fragments purified with the NZYGelpure kit (NZYtech, Lisbon, Portugal), cloned into pCR™2.1-TOPO^® vectors, propagated in One Shot^® Mach1™-T1R competent cells (Invitrogen, Life Technologies, Carlsbad, CA) and sent for sequencing (Eurofins Genomics, Ebersberg, Germany). Both strands were sequenced, and chromatograms were analyzed in FinchTV (Geospiza, Seattle WA, USA) and assembled using Multalin (http://bioinfo.genopole-toulouse.prd.fr/multalin/multalin.html).

Genomic DNA was isolated from sea bass red blood cells using the NZY Blood gDNA Isolation kit (NZYTech), according to the manufacturer’s instructions; quantification was performed using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific); quality was checked by agarose gel electrophoresis (22, 23). One microgram of genomic DNA was amplified by RT-PCR with the primers based on the previously obtained cDNA sequences, with the following cycling profile: 94°C for 5 min, 30 cycles of 94°C for 90 s, 59°C for 60 s, and 72°C for 90 s. Several PCR products were purified, cloned, and sent for sequencing as described earlier. Comparisons were made between cDNA and genomic DNA to assess the similarity of the coding regions and to identify intron/exon boundaries. A comparison between the genomic sequences of sea bass erythroferrone with those of other vertebrates was made, with Ensembl accession numbers for sequences used as follows: Homo sapiens (ENST00000546354); Mus musculus (ENSMUST00000086861); Monodelphis domestica (ENSMODT00000020244); Anolis carolinensis (ENSACAT00000006555); Ornithorhynchus anatinus (ENSOANT00000024567); Xenopus tropicalis (ENSXETT00000049409); Sparus aurata (ENSSAUT00010007639); Gasterosteus aculeatus (ENSGACT00000002111); Esox lucius (ENSELUT00000026569); Cyprinodon variegatus (ENSCVAT00000011742); Parus major (ENSPMJT00000028878); Serinus canaria (ENSSCAT00000008685).

The 5’ and 3’ RACE were carried out as previously described (22, 23) using the 5’/3’ RACE Kit, 2nd Generation (Roche Applied Science, Amadora, Portugal) according to the manufacturer’s instructions. Conditions for PCR were as follows: 94°C for 2 min, 94°C for 15 s, 59°C for 30 s, 72°C for 40 s, for 10 cycles; 94°C for 15 s, 59°C for 30 s, 72°C for 40 s (plus 20 s/cycle), for 25 cycles, with a final elongation at 72°C for 7 min. When necessary, a second PCR amplification was performed using the same conditions for an additional 30 cycles. Amplification products were run on agarose gels, relevant fragments purified, cloned, and sequenced as previously described.

Amino acid sequence alignments were performed using CLUSTALW from MEGAX (24). A phylogenetic tree was constructed using the Maximum Likelihood method, with the Jones–Taylor–Thornton (JTT) model, Nearest-Neighbor-Interchange heuristic model, complete deletion of gaps, and 10000 bootstrap replications. Sequences used for comparisons and phylogenetic trees include organisms representative of the various teleost fish and tetrapod classes, and their accession numbers were as follows: Homo sapiens (NP_001278761); Mus musculus (NP_775571); Monodelphis domestica (XP_007502443); Monodon monoceros (XP_029081431); Camelus dromedaries (KAB1279206); Ornithorhynchus anatinus (XP_001513652); Corvus cornix cornix (XP_010403787); Taeniopygia guttata (XP_030136162); Parus major (XP_015492724); Alligator mississippiensis (XP_014454239); Anolis carolinensis (XP_008104647); Xenopus tropicalis (NP_001072387); Carassius auratus (XP_026137994); Cyprinus carpio (KTG44855); Danio rerio (XP_002660750); Electrophorus electricus (XP_026885811); Esox lucius (XP_010890517); Fundulus heteroclitus (XP_021170161); Gadus morhua (XP_030219798); Ictalurus punctatus (XP_017351597); Lates calcarifer (XP_018534513); Larimichthys crocea (XP_019134833); Maylandia zebra (XP_004572440); Oncorhynchus kisutch (XP_020323002); Oncorhynchus mykiss (XP_021468332); Oreochromis niloticus (XP_013121495); Oryzias latipes (XP_004078828); Paralichthys olivaceus (XP_019965497); Salmo salar (XP_013994632); Salmo truta (XP_029559874); Seriola lalandi dorsalis (XP_023282631); Sparus aurata (XP_030259008)

Prediction of the three-dimensional structure of sea bass erythroferrone was performed by protein homology detection/modeling, using the Phyre2 Protein Fold Recognition Server (25). To maximize sequence coverage and confidence, intensive modeling was selected. Obtained models were submitted to the SAVES metaserver (http://services.mbi.ucla.edu/SAVES/) to check and validate protein structures. 96.2% of the residues were found to be in the most favored and allowed regions of the Ramachandran plot. Molecular graphic images were obtained using the Polyview-3D webserver (http://polyview.cchmc.org/polyview3d.html) (26).

Total RNA was isolated from tissues with the NZY Total RNA Isolation kit protocol for tissue samples (NZYtech, Lisboa, Portugal) with the optional on-column DNase treatment, according to the manufacturer’s instructions (22, 23). Total RNA quantification was performed using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific), and quality was assessed by running the samples in an Experion Automated Electrophoresis Station (Bio-Rad, Hercules, CA). For all samples, 2.5 µg of each were converted to cDNA using the NZY First-Strand cDNA Synthesis Kit (NZYTech) according to the manufacturer’s protocol.

Five healthy sea bass were sacrificed and several tissues (namely, liver, spleen, head kidney, intestine, gill, heart and brain) were collected for RNA isolation and cDNA synthesis, as previously described (22, 23). Relative levels of erfe mRNA were quantified by real-time PCR analysis using a CFX384 Touch Real-Time PCR Detection System (Bio-Rad). A total of 1 µL of each cDNA sample was added to a reaction mix containing 7.5 µL iTaq Universal SYBR Green Supermix (Bio-Rad), 5 µL double distilled H₂O, and 250 nM of each primer ( Supplementary Table 1 ), making a total volume of 15 µL per reaction. A non-template control was included for each set of primers. The cycling profile was the following: 95°C for 3.5 min, 40 cycles of 95°C for 20 s and 59°C for 20 s. Samples were prepared in duplicates, a melting curve was generated for every PCR product to confirm the specificity of the assays, and a dilution series was prepared to check the efficiency of the reactions. Beta actin (actb) was used as the housekeeping gene. The comparative CT method (2^-ΔΔCT method) based on cycle threshold values was used to analyze gene expression levels.

Anemia - To induce experimental anemias, fish were individually weighted and phlebotomized from the caudal vessels the equivalent v/w of or 2% body mass (21, 22, 27), while control fish were subjected to the same manipulation (anesthesia, weighting, pinching), but no blood was removed (for a total of 80 animals randomly divided in 2 groups, 40 fish per group). Iron overload - To induce iron overload ( 21, 22, 28 ), fish were intraperitoneally injected with 100 µl iron dextran (Sigma-Aldrich, St. Louis MO, USA) diluted in sterile PBS to a final concentration of 50 mg/mL, while control fish were injected with 100 µL sterile PBS (for a total of 80 animals randomly divided in 2 groups, 40 fish per group). Infection – Fish were intraperitoneally injected with 10⁵ CFU of Photobacterium damselae spp.piscicida, strain PP3, diluted in 100 µL of sterile PBS (for a total of 100 animals randomly divided in 2 groups, 50 fish per group). Control fish were similarly injected with 100 μL of sterile PBS. Hepcidin administration - Fish were intraperitoneally injected with commercially synthesized sea bass hepcidin peptides, either with 100 µL of a 50 µM solution of Hamp1 (QSHLSLCRWCCNCCRGNKGCGFCCKF), or Hamp2 (HSSPGGCRFCCNCCPNMSGCGVCCRF) (Bachem AG, Bubendorf, Switzerland), diluted in sterile PBS. Control fish were similarly injected with 100 µL of sterile PBS (for a total of 120 animals randomly divided in 3 groups, 40 fish per group).

In iron modulation and peptide models, fish were collected at 1, 4, 7, 10, 14 and 21 days after treatment. In the infection model, fish were collected at 1, 2, 3, 4 and 7 days post infection. Five fish from each of the experimental groups were euthanized with an overdose of anesthetic, blood collected for hematocrit determination, dissected, their tissues excised, snap frozen in liquid nitrogen and stored at -80°C until further use. RNA was isolated and converted into cDNA, and gene expression analysis (erfe, epo, hamp1, hamp2 and hbb) was conducted as previously described for the liver, spleen and head kidney. No animals were excluded in any of the experiments.

Statistical analysis was carried out using GraphPad Prism 8 (GraphPad Software Inc, La Jolla CA, USA). The normality of the distribution was evaluated with the Kolmogorov-Smirnov test. Multiple comparisons were performed with One-way ANOVA and post hoc Student Newman-Keuls test. A p value <0.05 was considered statistically significant.

A single sea bass erythroferrone transcript was obtained by PCR amplification and 5’/3’ RACE. Erythroferrone coding DNA (deposited in GenBank under accession number MT043300) consists of an open reading frame of 918 bp, with flanking 5’ and 3’ untranslated regions of 356 bp and 2090 bp, respectively, and encodes a 305 amino acid protein. The peptide has a predicted MW of 32870 Da and an isoelectric point of 5.29. A signal peptide cleavage site was predicted between positions 36 and 37 (MEA-AS) with the SignalP-5.0 server (http://www.cbs.dtu.dk/services/SignalP/), and several protein features were identified using Prosite (https://prosite.expasy.org/scanprosite/): a complement component 1q (C1q) domain profile (between positions 152-305), a proline rich region (105-125) and three N-glycosylation sites (197-200, 246-249 and 284-287) ( Supplementary Figure 1 ).

The genomic structure of sea bass erythroferrone consists of eight exons and seven introns, with a single initiator methionine and stop codon. The genomic interval from the initiator methionine to the stop codon is 6952 bp. Genomic organization of sea bass erythroferrone was analyzed and compared with those of other vertebrates, including mammals, amphibians, birds, reptiles and other fish. Significant variations were observed in terms of exon/intron numbers, as well as sizes, not only between vertebrate classes, but also between different fish species ( Figure 1 ).

Sequence comparison with other vertebrate species erythroferrone proteins showed a high degree of conservation particularly among fish species, with sequence identity varying between 57.48% (D. rerio) and 93.46% (L. crocea), while the lowest identity observed was with birds, of 26.47% (P. major) ( Supplementary Table 2 ). The highest degree of sequence conservation is seen at two of the conserved features, the proline rich region and the C1q domain profile ( Figure 2 ).

Additionally, although we can observe some conservation of the N-terminus among the various fish species, this region is less conserved when compared with tetrapods.

Phylogenetic analysis clusters fish erythroferrones separated from tetrapods ( Figure 3 ) and among fish, there seem to be three distinct clusters. One cluster includes members of the Cypriniformes and Gymnotiformes orders, a second cluster includes Salmoniformes and Esociformes and a third and more heterogeneous cluster, where sea bass erythroferrone is located, includes fish from the Perciformes, Pleuronectiformes, Beloniformes, Cyprinodontiformes, Cichliformes, Carangiformes and Gadiformes orders.

To predict protein structure, a profile-profile alignment webserver (Phyre2) was used to search for sequence similarities between sea bass erythroferrone and other proteins with known 3D structures that could serve as template(s) for homology modeling. Only the region corresponding to the C1q domain (comprised by amino acids 152-305) was successfully modelled, as of yet there are no available templates for a full length erythroferrone protein. Phyre2 selected 6 templates with a confidence level of 100% to model the protein: (1) globular head of the complement system protein C1q, PDB d1pk6c, 19% identity; (2) single chain recombinant globular head of the complement system protein C1q, PDB c5hkjA, 19% identity; (3) structure of the human collagen x nc1 trimer, PDB c1gr3A, 19% identity; (4) structure of the human collagen x nc1 trimer, PDB d1gr3a, 19% identity; (5) crystal structure of a single-chain trimer of human adiponectin globular domain, PDB c4douA, 18% identity; (6) globular domain of zebrafish complement 1qA protein, PDB c5hbaA, 11% identity. From a total of 154 residues, 136 (88%) were modeled at 100% confidence, whereas the remaining 18 (12%) were modeled ab initio. Secondary structure predictions reported by Phyre2 indicate a beta sheet configuration, with several disordered areas connecting the various beta strands ( Figure 4 ).

Constitutive expression of sea bass erythroferrone was evaluated in several healthy sea bass tissues, namely, liver, spleen, head kidney, intestine, gill, heart and brain ( Figure 5 ). Highest expression was observed in the spleen, followed by the head kidney, which is not surprising, considering these are the major hematopoietic organs in teleost fish. Moderate expression was observed in the intestine, gill, heart and brain, with the liver presenting the lowest constitutive expression of all tested organs.

Gene expression was evaluated in the liver, spleen and head kidney of anemic fish, at day 1, 4, 7, 10, 14, and 21 post anemia induction. As expected, a significant decrease in the hematocrit was observed as early as day 1, with a gradual recovery towards day 21 but still not completely recovering to control levels ( Figure 6A ). This anemia was accompanied by significant decreases in the expression of hamp1 in all organs, particularly in the earlier days, with later gradual recoveries to control levels, whereas no significant variations were observed for hamp2 ( Figures 6B-D ). Erythroferrone (erfe) levels in the spleen ( Figure 6E ) and head kidney ( Figure 6H ) mirrored hamp1, with substantial increases in expression, although reaching peak levels at different time points, with the zenith at day 4 in the spleen and 14 in the head kidney, followed in both by decreases and recovery to control levels. Increased expression of erythropoietin was also observed in the spleen ( Figure 6F ) and head kidney ( Figure 6I ), also with different patterns of expression, with a gradual increase up to day 21 in the spleen, while reaching the zenith at day 10 in the head kidney, followed by a gradual decrease to control levels. Additionally, hemoglobin levels mirrored the hematocrit, with gradual increases in both the spleen ( Figure 6G ) and head kidney ( Figure 6J ), up to day 7 and 10 respectively, followed by a return to control levels at the same time the hematocrit returned to near normal levels.

Gene expression was evaluated in the liver, spleen and head kidney of iron overloaded fish, at day 1, 4, 7, 10, 14, and 21 post iron overload. A steady increase in the hematocrit was observed up to day 7, followed by a gradual decrease to control levels ( Figure 7A ). Considerable increases in expression were observed for hamp1 in all tested tissues ( Figures 7B-D ), with similar patterns of early increase in the liver and spleen (both peaking at day 4), but a bit later in the head kidney (with the zenith at day 7), but nevertheless with recoveries to control levels in all organs at day 21 post-overload. For hamp2, no significant variations were observed. No significant variations in erfe expression were observed in the spleen ( Figure 7E ) or head kidney ( Figure 7H ), despite the slight decreases observed in erythropoietin expression in the spleen, at day 4 ( Figure 7F ), and in head kidney, at days 4 and 7 ( Figure 7I ). Furthermore, significant increases in the expression of hemoglobin were observed both in the spleen ( Figure 7G ) and head kidney ( Figure 7J ) at day 7, with a return to normal levels at day 10.

To evaluate a possible involvement of erythroferrone during infection, gene expression was evaluated in the liver, spleen and head kidney of fish infected with the Gram-negative bacterium P. damselae spp. piscicida, 1, 2, 3, 4 and 7 post-infection. Infection had a significant impact on hematocrit, with a gradual decrease up to 4 days, followed by a recovery but still not returning to normal levels after 7 days ( Figure 8A ). Early increases in hamp1 were observed in all tested tissues, followed by a rapid decline towards normal levels. Hamp2 expression on the other hand was highly increased in the liver, up to 7 days, with a similar but more moderate pattern of overexpression in the head kidney, as well as a slight increase in the spleen, at 4 days post-infection ( Figures 8B-D ). Similar patterns of downregulation were observed for erfe, epo and hbb, in both spleen ( Figures 8E-G ) and head kidney ( Figures 8H-J ), with levels kept under normal levels mostly between days 1 and 4, and with a recovery of all genes to normal levels of expression at day 7.

It is known that erythroferrone has an effect on hepcidin, but we also wanted to know if hepcidin could have an impact on erythroferrone. As such, erythroferrone, erythropoietin and hemoglobin expressions were evaluated in the spleen and head kidney of fish administered Hamp1 and Hamp2 synthetic peptides, at day 1, 4, 7, 10, 14, and 21 post-administration. Administration of Hamp1 had a significant impact on the hematocrit, leading to a reduction starting as early as day 1, reaching the nadir at day 4 and then gradually recovering to control levels up to day 21 ( Figure 9A ); Hamp2, on the other hand, had no impact on hematocrit levels. Erfe expression was only influenced by the administration of Hamp1, with Hamp2 having no significant effects. Both the spleen ( Figure 9B ) and kidney ( Figure 9E ) presented similar patterns of under expression up to day 10 post-administration, but whereas in the kidney expression levels slowly recovered to control levels, in the spleen they were surpassed at day 14, with a significant overexpression followed by recovery to control levels at day 21. Similarly, epo and hbb expressions were also only influenced by Hamp1 administration, with both genes being progressively downregulated in the spleen ( Figures 9C, D ) and head kidney ( Figures 9F, G ) up to day 7, followed by gradual recoveries to normal levels.