RNA-seq analysis of hypoxia-treated and control zebrafish larvae detected a total of 2411 differentially expressed genes (DEGs), including 703 up-regulated and 1708 down-regulated DEGs (Figure 1A). The two groups are distinctly clustered in hierarchical clustering (Figure 1B) and principal component analysis (Supplementary Figure S1A). The RNA-seq results were confirmed by qRT-PCR assays using a different set of biological samples (Figure 1C) and similar changes were detected in 18 out of 20 genes tested (Supplementary Figure S1B).

KEGG analysis indicated that the up-regulated DEGs are highly enriched in the glycolysis/gluconeogenesis pathway, followed by carbon metabolism, amino acid synthesis etc. (Figure 1D). The topmost enriched pathways, however, are “neuroactive ligand signaling” and “hormonal signaling” (Figure 1D). In addition, “calcium signaling”, “FoxO signaling”, and “insulin signaling” pathways are also enriched (Figure 1D). In agreement, GO (Molecular Function) analysis identified “signaling receptor regulator activity” as the top enriched term in the up-regulated DEGs (Figure 1E). “Hormone activity”, largely overlapping with “signaling receptor regulator activity “, is also enriched (Figure 1E). Other enriched GO terms are oxidoreductase activity, iron ion binding, dioxygenase activity, FAD binding etc. (Figure 1E). The enriched down-regulated KEGG pathways include biosynthesis of co-factor, cornified envelope formation, cytoskeleton in muscle cells, focal adhesion etc. (Figure 1F). The top enriched down-regulated GO terms are endopeptidase activity, glycosyltransferase activity, peptidase regulator activity, iron ion binding, enzyme inhibitor activity etc. (Figure 1G). GO CNET plot of the up-regulated DEGs revealed the top five most significantly enriched molecular function terms are dioxygenase activity, 2-oxoglutarate-dependent dioxygenase activity, hormone activity, oxidoreductase activity acting on paired donors with incorporation or reduction of oxygen, and FAD binding. Correspondingly, the top five enriched terms among down-regulated DEGs were endopeptidase activity, peptidase regulator activity, endopeptidase regulator activity, peptidase inhibitor activity, and monooxygenase activity (Supplementary Figures S1C, D).

We further analyzed the RNA-seq dataset by ranking the up-regulated DEGs in hormone activity according to their mRNA abundance (Figure 2A). Among the top on the list is stc2a, which has been implicated in human body height regulation and IGF signaling (22–24) During early development, stc2a mRNA levels increased gradually and reached the plateau at 4 days post fertilization (dpf) (Figure 2B). Data mining from the Zebrahub scRNA-seq database (25), indicated that stc2a mRNA is detected in multiple tissues from 0 somite stages to 10 dpf, and is particularly abundant in the central nervous system, periderm, neural crest, paraxial mesoderm (Figure 2C). In the adult stage, the highest stc2a mRNA levels were detected in the brain, followed by kidney, eyes, and fin. (Figure 2D). Importantly, hypoxia increased stc2a mRNA levels at multiple developmental stages (Figure 2E). The hypoxia-induced stc2a mRNA expression was attenuated in hif2a^-/- deficient fish (Figure 2E), suggesting that stc2a expression is induced by hypoxia, likely through the action of Hif2.

Using CRISPR-Cas9, two mutant zebrafish lines, stc2a(Δ2 + 4)^-/- and stc2a(Δ5)^-/-, were generated for functional analysis. Both are predicted to be null mutations. In both mutant lines, there was a significant reduction in the stc2a mRNA levels, whereas no such changes were found with the mRNA levels of the closely related stc2b (Supplementary Figure S2B). Both mutant fish lines survived to adulthood and reproduced well under standard conditions. The gross morphology of these stc2a^-/- fish was indistinguishable from their siblings (Figure 3A). Compared with wild-type fish, however, the body length of stc2a(Δ2 + 4)^-/- and stc2a(Δ5)^-/- larvae was significantly greater (Figure 3B). Likewise, the head-trunk angle (HTA) and somite number, two parameters of zebrafish developmental speed (73), were significantly greater (Figures 3C, D), suggesting these mutant fish develop more rapidly and grow faster. Next, stc2a(Δ2 + 4)^+/- fish were intercrossed and the offsprings were grew under the same conditions for 6 months. Compared with their wild-type and heterozygous siblings, stc2a(Δ2 + 4)^-/- fish had greater body length and body weight (Figures 3E, F). Likewise, the lengths of the pelvic and dorsal spines in 1 year-old stc2a(Δ2 + 4)^-/- fish were greater than those of the siblings (Figures 3G, H). These results suggest that Stc2a negatively regulates somatic growth, organ size, and body size in adult fish.

Since genetic deletion of stc1a, a structurally related gene, resulted in increased ionocyte cell number, ectopic calcium deposits, kidney stone-like calcium deposits, and reduced bone mineralization (26, 27), we crossed stc2a(Δ5)^-/- fish with Tg (igfbp5a:GFP) fish, a stable transgenic line expressing EGFP in calcium transporting ionocytes or NaR cells (28), and measured the NaR cell number. There were no differences among the different genotypes (Supplementary Figures S3A, B). Alizarin red staining of the juvenile and adult stc2a(Δ2 + 4)^-/- fish did not detect notable differences in calcified tissues (Supplementary Figures S3C, D), suggesting that Stc2a does not play an indispensable role in regulating ionocyte proliferation, calcium balance, bone mineralization nor kidney development or function.

We postulated that Stc2a may regulate somatic growth by inhibiting the pappalysin family metalloproteinase (including PAPP-A and PAPP-A2) and IGF signaling activity. In zebrafish, there are 3 pappalysin family members, including Papp-aa, Papp-ab, and Papp-a2 (29, 30). Currently, there is no specific antibodies against any of these proteins. qRT-PCT analysis showed the levels of papp-aa and papp-a2 mRNA levels were similar between wild-type fish and stc2a^-/- fish (Supplementary Figure S4A). Although there was a markedly increase in papp-ab mRNA levels, this change was not statistically different (Supplementary Figure S4). Hypoxia did not alter their mRNA levels of papp-aa and papp-ab, while it caused a modest but statistically significant decrease in papp-a2 mRNA levels. (Supplementary Figure S4). Commercially antibodies against mammalian phospho-IGF1 receptor did not yield specific signal in zebrafish larvae. Western blotting of the whole body lysates did not detect major differences in phospho-Akt and phospho-Erk levels between stc2a(Δ2 + 4)^-/- and wild-type fish (Figures 4A, B). The low sensitivity of the assay and the lack of tissue/cell resolution may obscure the results.

We therefore tested the possible role of IGF signaling using BMS-754807, a potent IGF1R inhibitor that has been previously tested in zebrafish (11, 31–33). If Stc2a regulates somatic growth by inhibiting IGF signaling, then inhibition of IGF signaling should reverse these phenotypes. Treatment with BMS-754807 decreased body lengths in both stc2a(Δ2 + 4)^-/- fish and wild-type control fish. The magnitude of decrease, however, was significantly greater in the stc2a((Δ2 + 4)^-/- fish (Figure 4C). Likewise, ZnCl₂, which inhibits pappalysin-mediated Igfbp degradation in human cells and in zebrafish (34), caused a greater percent decrease in body length in the stc2a((Δ2 + 4)^-/- larvae compared to wild-type siblings (Figure 4D). To further test the role of IGF signaling, the mutant fish were treated with wortmannin, a PI3K inhibitor (35), and U0126, a MAPK inhibitor (36). Both inhibitors decreased body length in the wild-type as well as stc2a(Δ2 + 4)^-/- fish. Again, the magnitude of decrease was significantly more pronounced in the stc2a((Δ2 + 4)^-/- fish (Figures 4E, F).

Hypoxia caused growth retardation and developmental delays in wild-type zebrafish embryos (Figure 5A) (8, 11). Although hypoxia reduced the body length and HTA in stc2a(Δ2 + 4)^-/- (Figure 5A), the mutant fish still grew bigger compared to wild-type siblings under hypoxia (Figure 5A). Likewise, while hypoxia lowered HTA values in both genotypes (Figure 6B), the mutant fish developed faster (Figure 5B). Treatment of mutant fish with ZnCl₂ reduced the body length and HTA value to the levels of wild-type fish (Figures 5C, D). We noted that mutant fish were prone to die under hypoxia and quantified their survival rate under hypoxia. They were significantly lower compared to the wild-type fish (Figure 5E, Supplementary Table S1). Addition of ZnCl₂ restored the survival rate to the levels of wild-type control groups (Figure 5F, Supplementary Table S1), suggesting that pappalysin enzyme activity is required in the elevated growth and mortality in stc2a^-/- mutant fish under hypoxic stress.

Because zebrafish larvae are tiny, several dozen were pooled for Western blotting analysis of phospho-Akt and phospho-Erk levels. No difference was detected between stc2a^-/- and wild-type fish. This was likely due to low sensitivity of Western blotting and the lack of resolution. Future studies are needed to develop more sensitive and quantitative assays to detect local IGF signaling activities. In this study, ZnCl₂ was used as a generic pappalysin metalloproteinase inhibitor. The specific pappalysin metalloproteinase isoform(s) involved with the reported Stc2a action is unclear. qRT-PCR analysis did not detect major changes in papp-aa and papp-a2 mRNA levels in stc2a-/- mutant fish. While there was a trend of increase in papp-ab mRNA levels, it was not statistically significant. The actual protein levels and their activities are currently unclear. Among the 3 zebrafish pappalysin members, Papp-aa could cleave human IGFBP5 and IGFBP4 in vitro, whereas it did not cleave other 4 human IGFBPs (34). Zebrafish Papp-ab has been shown to cleave human IGFBP4 and IGFBP5 but not the other four IGFBPs (29). Zebrafish Papp-a2 can cleave human IGFBP3 and IGFBP5, but not other IGFBPs (30). Zebrafish has a total of 9 igfbp genes, including igfbp1a, 1b, 2a, 2b, 3, 5a, 5b, 6a, and 6b (8, 10, 18, 19, 54–56). The specific Igfbp(s) involved in the Stc2a action reported in this study remains unknown. Future studies will be needed to develop reagents and tools to measure these pappalysin metalloproteinase activities, identify their substrates, and develop conditional knockout and isoform-specific inhibitors to dissect their role(s) in somatic growth and survival trade-off under hypoxia stress.

Zebrafish were maintained following standard zebrafish husbandry guidelines (57). All experiments using zebrafish were conducted in line with guidelines approved by the Institutional Animal Care & Use Committee, University of Michigan. Embryos and larvae were raised in standard E3 medium as reported previously (58). 0.003% (w/v) N-phenylthiourea (PTU) was added to the E3 medium to prevent pigmentation when required. Modified low calcium media was prepared following a previously reported protocol (32). In addition to wild-type (WT) fish, Tg(igfbp5a:GFP) fish (28), hif2αbΔ10^-/- (e.g., epas1b.2, (59)), stc1a^-/- (26), stc2a(Δ2 + 4)^-/- and stc2a(Δ5)^-/- (this study) were used.

Zebrafish larvae were subjected to hypoxia (6% O₂) or normoxia (20.9% O₂, atmospheric level) from 81 to 96 hours post-fertilization (hpf). Each group comprised three to four biological replicates, with 30 larvae pooled per replicate. Total RNA was extracted using the PureLink™ RNA Mini Kit (Invitrogen) and treated with DNase using the PureLink™ DNase Set (Invitrogen). Independent RNA sample sets (n = 3~4) were submitted for library preparation and sequencing at the Advanced Genomics Core, University of Michigan. Poly(A) RNA libraries were constructed using the NEBNext Poly(A) mRNA Magnetic Isolation Module and NEBNext UltraExpress RNA Library Prep Kit (New England Biolabs). Sequencing was performed on the Illumina NovaSeq X Plus platform, generating 89.5–104.6 million paired-end reads (151 bp) per sample. Raw sequencing reads were processed to remove low-quality bases and adapter sequences using Trimmomatic (v0.39) (60), and read quality was evaluated before and after trimming using FastQC (v0.12.1) (61). Trimmed reads were then aligned to the zebrafish reference genome (GRCz11, Ensembl release 113) using STAR (v2.7.10) (62), with gene annotations obtained from Ensembl (63). Read counts were obtained using featureCounts (v2.0.7) within the Subread package (64). DEseq2 (v1.46.0) was used to analyze the differential expression between the groups (65). Genes with adjusted p-value (padj) < 0.05 and |log₂ fold change| > 0.6 were considered differentially expressed genes (DEGs). Gene Ontology (GO) and KEGG pathway enrichment analyses of DEGs were performed using the org.Dr.eg.db annotation database (v3.20.0) (Carlson) and the clusterProfiler R package (v4.14.4) (66). Enriched terms with p-value and q-value < 0.05 were considered statistically significant. Top GO terms and KEGG pathways were visualized using the dot plot function, and gene-concept networks (CNETs) were generated using the cnetplot function, both implemented in the clusterProfiler R package. All statistical analyses were conducted using R software (v4.4.1).

RNA was isolated from a pool of 15~30 zebrafish larvae or from adult zebrafish tissue as reported (67). RNA was reverse-transcribed to cDNA using oligo-dT primers and M-MLV reverse transcriptase (Invitrogen). qPCR was performed using SYBR Green (Bio-Rad) on a StepONE PLUS real-time thermocycler (Applied Biosystems) as previously reported (68). The expression level of a target gene transcript was normalized by 18s rRNA or β-actin mRNA levels. PCR primers were designed based on sequences as described in previous studies, the NCBI Gene database, and by NCBI Primer Blast (27, 69) and are listed in Supplementary Table S2. The spatial expression information of stc2a mRNA was extracted from the Single-cell RNA-seq dataset from ZebraHub (https://zebrahub.org) using Scanpy (v1.11.1) (25). Expression of stc2a mRNA was extracted and averaged across tissue categories defined by annotation. Mean expression values were computed per tissue, and relative expression was calculated by normalizing to the sum of expression across all tissues. Visualization and downstream analysis were performed using Seaborn and Matplotlib.

Two sgRNAs targeting the stc2a gene were designed using CHOPCHOP (http://chopchop.cbu.uib.no/). Their sequences are: stc2a-gRNA1-oligo1: 5’- GCTGCTGCTCTCCGTATTGG-3’ and stc2a-gRNA3-oligo3: 5’-GGGTGACTCTCGTGCACATC-3’. sgRNA(30–40 ng/ul) mixed with Cas9 mRNA (200–400 ng/ul) were injected into Tg(igfbp5a:GFP) at the 1-cell stage (70). The injected F0 fish were raised to adulthood and crossed with Tg(igfbp5a:GFP) fish (28). The DNA of the F1 fish were extracted by fin clipping and analyzed by sanger sequencing. Heterozygous F1 male and female fish were crossed to generate the F2 fish.

Genomic DNA, isolated from adult fin or whole larval lysate, were digested with proteinase K (60 μg/mL) in SZL buffer (50 mM KCl, 2.5 mM MgCl₂, 10 mM Tris-HCl (pH 8.3), 0.45% NP-40, 0.45% Tween 20, 0.01% gelatine). Samples were digested at 60°C for 2 hours and at 95°C for 15 minutes. The stc2a^-/- fish genotyping was performed by PCR and by direct DNA sequencing as previously reported (71).

The bright-field images of larvae zebrafish were acquired using a stereomicroscope (Leica MZ16F, Leica, Wetzlar, Germany) equipped with a QImaging QICAM camera (QImaging, Surrey, BC, Canada). Head-trunk-angle, adult zebrafish body length, body weight, and brain weight was measured following published protocols (8, 11). Larval body length measured from inner ear stone to tail end. Image J was used for image analysis and data quantification.

Alizarin red staining was performed as described previously (72). The bright-field images were acquired as described above. The whole-body images were joined together with Adobe Photoshop. Dorsal and pelvic spine lengths measured by ImageJ.

All drugs used in this study, except ZnCl₂, were dissolved in DMSO and further diluted to desired concentrations. ZnCl₂ was dissolved in distilled water. Zebrafish larvae were treated with drugs and vehicle as described previously (28). Drug solutions were changed daily.

Hypoxia treatments were performed using the Invivo2–300 Hypoxia Workstation with an I-CO₂N₂IC advanced gas mixing system (Baker Ruskinn, Sanford, ME). Gas tanks were purchased from Cryogenic Gases (Detroit, MI). CO₂ levels were kept constant at 1.4%. In addition to monitoring the gas values shown by this machine, an oxygen meter was utilized for a secondary reading (Sper Scientific, Scottsdale, AZ). The temperature was maintained at ~ 28°C. Zebrafish larvae were set up in 6-well plates at a density of 15 larvae per well with 3–5 mL of E3 or drug medium. The number of dead larvae was counted hourly based on appearance and motor response to stimulus.

Zebrafish larvae (30–40 per group) were homogenized in RIPA buffer (50 mM Tris-HCl, 150 mM NaCl, 1 mM EGTA, 0.1% Triton X-100, 0.5% sodium deoxycholate, 0.1% SDS, pH 7.4) containing a cocktail of protease inhibitors and phosphatase inhibitors on ice for 30 seconds. Zebrafish homogenates were centrifuged at 13,000 rpm, 4°C, for 20 minutes. The protein concentration of the supernatant was measured by Bradford assay (Bio-rad) and normalized. 2X urea loading buffer (150 mM Tris pH 6.8, 6 M Urea, 6% SDS, 40% glycerol, 100 mM DTT, 0.1% Bromophenol blue). Samples were heated at 95°C for 5 minutes and subjected to 12% SDS-PAGE gels and transferred to a nitrocellulose membrane for western blot analysis. After incubation with primary and secondary antibodies, membranes were scanned using the Odyssey CLx imaging system (LI-COR). The ratios of phosphorylated to total protein were calculated after the images had been analyzed. The primary antibodies were mouse anti-GAPDH (1:4000, Proteintech), rabbit anti-Akt (1:1000, Cell Signaling Technology), rabbit anti-P-Akt (S473) (D9E) XP^® (1:2000, Cell Signaling Technology), rabbit anti-p44/42 MAPK (1:2000, Cell Signaling Technology), rabbit anti-Phospho-p44/42 MAPK (D13.14.4E) XP^® (1:2000, CellSignaling Technology). The secondary antibodies, goat anti-mouse IRDye 680LT and goat anti-rabbit IRDye 800CW were purchased from LI-COR Biosciences and used at 1:10,000 dilution.

Statistical analysis was performed using GraphPad Prism 9 software (GraphPad Software, Inc., San Diego, CA). Values are shown as means ± SEM. Statistical significance between experimental groups was performed using unpaired two-tailed t-test, one-way ANOVA followed by Tukey’s multiple comparison test, or two-way ANOVA with multiple comparisons. Differences between the experimental groups in survival curves were analyzed using the Mantel-Cox log-rank test. Statistical significances were accepted at *p < 0.05, **p<0.01, ***p < 0.001, ****p < 0.0001.