Knockdown of grb10a expression by splice-blocking antisense oligonucleotide

As grb10 has been linked to embryonic growth trajectory, expression was examined over the first 120 hpf. QPCR analysis of grb10a expression at 24-h intervals revealed a strong upregulation at 48 hpf (Fig. 2a). Expression of the grb10 paralogue, grb10b, was not detectable at any time point.

Fig. 2

Grb10a is successfully knocked down in zebrafish injected with splice-blocking antisense oligonucleotides. (A) Grb10a qPCR of WT embryos (24–120 hpf, triplicated, n = 5 embryos per well). Data are shown as gene expression relative to β-actin with maximum levels at 48 hpf. (B) Schematic of the first five exons of the zebrafish grb10a gene. 5’ splice sites are highlighted with the forward and reverse primer triad indicated. (C) Multiplexed PCR amplification of the e3i3 and e4i4 splice site in embryos treated with either Standard Control (SC) morpholino, e3i3, or e4i4. β-actin was used as a positive control. (D) Western blot results of phosphorylated versus total protein ratios for two major signalling molecules of the insulin signalling pathway: AKT and S6. Quantitation using densitometry depicts mean + SEM. Activation of both proteins was found to be significantly elevated in KD zebrafish compared to SC (n = 3, unpaired t-test *** p = 0.0007, * p = 0.04). (E). Expression of grb10 mRNA throughout early life. Grb10 expression is decreased at 5 dpf in the KD but thereafter levels in the KD are not different from those in the SC (Wilcoxon rank sum test. **** p = 0.00003)

To knock down grb10a expression, zygotes were microinjected with splice-blocking antisense oligonucleotides e3i3 and e4i4. Exon 3 and 4 donor splice sites (Fig. 2b) were targeted to confirm the specificity of the phenotype, in accordance with current guidelines [30]. Multiplexed RT-PCR amplification using primers flanking the splice sites (Fig. 2b) showed a single product of the anticipated size for e3i3 and e4i4 embryos, and no product was detected for SC embryos (Fig. 2c), consistent with retention of the corresponding intron.

To confirm grb10a KD induced a quantifiable impact on the downstream insulin signalling pathway, phosphorylation of key proteins was analysed by Western Blot. As shown in Fig. 2d, phosphorylated (active) versus total protein ratios of AKT and S6 were significantly elevated in grb10a KD zebrafish at 96 hpf compared to SC (p = 0.0007 and 0.04 respectively, n = 3), consistent with the expected impact of grb10a KD.

To characterise the transient nature of the morpholino knockdown, grb10 expression was examined across the first thirty days of life. This analysis revealed that grb10 expression was reduced in the morpholino treated group at 5 dpf (p = 0.00003) but returned to control levels by 10 dpf (not significant [NS]) and remained consistent throughout the juvenile phase (Fig. 2e).

Growth trajectory and early life cardiometabolic phenotype is significantly impacted by early life grb10a knockdown

To determine the effect of grb10a KD on growth, total body length was measured at 24-h intervals over the first 5 dpf. As shown in Fig. 3a, total body length was initially comparable between KD and SC zebrafish (2.99 ± 0.06 mm vs 2.83 ± 0.1 mm, not significant (NS), n = 9 and n = 10, respectively). Subsequently, KD zebrafish began to diverge from the SCs at 48 hpf, corresponding to the peak in grb10a expression observed in WT zebrafish (Fig. 3a). KD zebrafish were longer on average than SCs (3.41 ± 0.02 mm vs 3.18 ± 0.02 mm at 72 hpf, p = 0.039, n = 46 and n = 41, respectively). This phenotype was reversed in zebrafish overexpressing grb10a, which were significantly shorter than SC counterparts (3.36 ± 0.02 mm vs 3.51 ± 0.03, p = 0.001, n = 24 and n = 25, respectively), as shown in Fig. 3b. Co-injection of e3i3 and grb10a RNA returned body length to SC levels (3.51 ± 0.03 mm vs 3.56 ± 0.03 mm, NS, n = 25), confirming the validity of e3i3 induced grb10a KD. Moreover, grb10a KD induced by e4i4 (Fig. 3c) also resulted in increased body length, and the ability of grb10a overexpression to suppress growth was shown to be dose dependent (Fig. 3c). Intriguingly, by 120 hpf, body length converged, possibly indicating activation of compensation to regulate growth and return to an “ideal” length post-hatch.

Fig. 3

Growth and cardiometabolic phenotype are significantly impacted by early life grb10a knockdown. (A) Mean total body length ± SEM of SC and KD zebrafish from 24 to 120 hpf (multiple comparison 2-way ANOVA, * p < 0.05, *** p = 0.0003). (B) Mean total body length ± SEM and individual data points of 96 hpf zebrafish embryos (n = 25). Grb10a KD phenotype was reversed in grb10a overexpression zebrafish. Co-injection resulted in phenotype rescue. One-way ANOVA revealed KD zebrafish were significantly longer, while grb10a overexpression zebrafish were significantly smaller than SC (*** p = 0.0001). Rescue zebrafish were of similar length to SC (p = 0.38). (C) Mean body length measurements of the KD relative to SC ± SEM from 48 to 120hpf. (D) Mean heart rate ± SEM in beats per minute of SC and KD zebrafish (multiple comparison 2-way ANOVA, **** p < 0.0001). (E) Mean yolk area ± SEM of SC and KD embryos over the embryonic life stage (multiple comparisons 2-way ANOVA, *** p = 0.0004. (F) Glucose Uptake-Glo.™ Assay of 96 hpf KD and SC zebrafish, where higher luminescence indicates a greater accumulation of intracellular 2D6P. Luminescence was approximately 30% greater in KD zebrafish compared to SC (unpaired t-test, *** p = 0.0002)

To investigate the impact of grb10a KD on the developing cardiac system, heart rate was measured over the first 5 dpf. As shown in Fig. 3d, while heart rate increased slightly over time in SC zebrafish, average KD heart rate fell. By 120 hpf, mean heart rate was almost 50% lower in the KD compared to SC (61 ± 7 bpm vs 118 ± 3 bpm, p < 0.0001, n = 14 and n = 19, respectively).

To determine whether grb10a KD had an impact on metabolic rate, yolk absorption was measured to indicate nutrient utilisation and energy demand. As shown in Fig. 3e, there was initially no difference in yolk area between the groups (NS, n = 10). SC and KD yolk consumption began to diverge at 48 hpf (0.30 ± 0.01 mm2 vs 0.39 ± 0.01 mm2, p < 0.0001, n = 18 and = 19, respectively). The greater yolk consumption observed in the KD fish suggests an elevated metabolic rate. To support this conclusion, a Glucose Uptake-Glo™ Assay was performed to compare the rate of glucose uptake, measured by 2D6P content. As shown in Fig. 3f, 2D6P accumulation was significantly higher in KD zebrafish compared to SC, an increase of almost 30% (p = 0.0002, n = 5), indicating glucose uptake was elevated. These growth and metabolic findings are consistent with the role of grb10a as a negative regulator of growth and insulin signalling pathway [92].

Early-life grb10a knockdown persistently dysregulates age-associated gene expression

In advance of transcriptomic analyses, tissue types represented in the dataset were quantified by assessing patterns of gene expression from an independent zebrafish scRNA-seq dataset annotated to tissue level [40]. As anticipated, the 5–30 dpf zebrafish samples most closely resembled the transcriptome of mid-hindbrain tissue (> 88%), with additional contribution from other neural tissues as well as cephalic and skeletal muscle (Supplementary Table 2).

To understand whether the changes observed in grb10a KD embryos were coupled with lasting changes in gene expression, the transcriptomic landscape of SC and KD zebrafish was investigated over the first 30 dpf. Unsupervised hierarchical clustering, standard deviation filtering, and maximised projection scores (MPS) were used to define a set of genes with strong age-association in the SC zebrafish (297 genes, MPS = 0.43). These genes fell into four distinct clusters, associating with 5, 10, 15, and 20–30 dpf (163, 15, 32, and 87 genes respectively) (Fig. 1ai and Fig. 4a).

Fig. 4

Transcriptomic analysis of standard control and grb10a knockdown gene expression over the first 30 dpf. Hierarchically clustered heat maps of gene expression generated from an Affymetrix GeneChip™ Zebrafish Genome Array of SC and KD zebrafish RNA, taken at 5, 10, 15, 20, and 30 dpf. (A) Expression of age-related genes in SC segregates into three clusters in SC zebrafish. (B) Clustering of the same age-related genes (identified in A) is disrupted in KD. Notably, the cluster of genes expressed at 5 dpf in SC are expressed at a variety of time points in KD. (C) Analysing the KD dataset independently highlights a set of genes with age related gene expression which fall into two clusters. (D). Venn diagram of age associated genes identified separately in SC and KD demonstrates that only a small proportion of SC age associated genes (30%, 90/297) remained age associated in KD. (E) Gene set enrichment analysis, using the GO Biological Process Ontology gene list, of the age related genes in the SC and KD datasets. The top 20 most enriched pathways with differential expression are included here, with associated normalised enrichment scores and q-values

Upon examining expression of the same genes in KD zebrafish, the clusters of co-expressed genes identified in the SC animals (Fig. 4a) was not present in the KD zebrafish (Fig. 4b). While the clustering was largely conserved in the latter three clusters (10 dpf– 14/15, 15 dpf– 29/32, and 20–30 dpf– 74/87 genes), genes strongly associated with 5 dpf in the SCs were generally expressed at different time points in the KDs (dotted white lines, Fig. 4a and Fig. 4b), with a mapping of only 38/163 genes (23%). Notably, five genes associated with 5 dpf in SC zebrafish mapped to 20–30 dpf in the KD. This suggests that the expression of genes usually associated with early larval development, is instead associated with the late-juvenile stage in KDs. Human orthologues of these dysregulated genes include DGAT2 (fatty acid metabolism), GAMT (energy storage, muscle contraction, and fatty acid oxidation), and PDIA2 (oxidative stress).

As the SC gene clusters were disrupted in the KD dataset, unsupervised analysis of the KD data was performed to identify the subset of age-associated genes in the KD zebrafish. 119 genes were identified (MPS = 0.37) which segregated into 5–15, 15–30, and 20–30 dpf (37, 16, and 66 genes) (Fig. 4c). Although 75% (90/119) of age associated genes in the KD were shown to have been age associated in SC, only a small proportion (30%, 90/297) of genes which were age associated in SC were also age associated in KD (Fig. 4d).

To assess functionality of the identified age-related genes (Fig. 1ai), Gene Set Enrichment Analysis was performed for genes identified in both the SC and KD datasets. Functionality conserved between the SC and KD is described in Supplementary Table 3. Dissimilar pathways (Fig. 4e) included several actin and collagen related pathways, and extracellular structure and RNA processing, which were age-related in KD but not SC zebrafish. Conversely, several metabolic pathways were age-associated in SC zebrafish but not in KD animals. This dysregulation of age-related gene expression in the KD zebrafish implies early life grb10a KD induces remodelling of the transcriptome.

We next assessed transcriptomic organisation using hypergraphs. We discerned a significant increase in the entropy of the transcriptome in KD compared with SC at 5-10dpf and 20-30dpf (FDR < 2.4 × 10–11) (Fig. 5a). Entropy is a measure of information content and serves as an indicator of “disorder”, such that lower entropy (more order) describes a network with little crosstalk and a more discrete function [64, 93], while a network with greater entropy (more disorder) has increased crosstalk and pleiotropic functions. This demonstrates that significant differences in the organisation of the transcriptome occur in the early-larval phase as a result of the grb10a KD with a further reorganisation during later development.

Fig. 5

Hypergraph entropy in the transcriptome reveals differences over time and between functions. (A) hypergraphs were iterated using randomly selected genes and a moving window approach to assess change in transcriptomic entropy over time. SC and KD are significantly different at all time points (FDR < 2.4 × 10–11 ANOVA with Tukey post-hoc) though the largest differences exist at 5–10 dpf (fold change = 0.052) and 20–30 dpf (0.11). (B) hypergraphs were iteratively generated from genes attributed to each pathway and entropy was modelled across the two treatment groups. β values represent the difference in entropy between SC and KD, with a β value of 0 indicating no difference between groups. Positive β values represent higher entropy in SC compared to KD. Ontology classes are considered to be significantly different if the 89% CI of β values does not include 0. Actin filament-based movement and negative regulation of extrinsic apoptotic signalling pathway were the only two pathways with no difference identified between the groups

By investigating pathways enriched in the age-related cluster (Fig. 4e), we were able to demonstrate their contribution to the organisation of the transcriptome by measuring the difference in entropy for each pathway between KD and SC (Fig. 1aii). Most of these pathways had higher entropy in SC than KD (Fig. 5b). The greatest differences in activity between the KD and SC were in Positive Regulation of Lipid Localization, Amino Acid Transport and Cellular Modified Amino Acid Biosynthetic Process. This indicates that entropy and gene coordination at the whole transcriptome level (higher in KD) are different at the level of these pathways (higher in SC).

By employing a random walk, transition between nodes in the network can be modelled probabilistically. Although this does not capture the mechanistic interactions between genes, it reflects flow and coordination in the system, to more closely reflect the causality of a molecular cascade. The structure of the random walk transition matrix was strongly related to the structure of the hypergraph adjacency matrix, as measured by correlation between the two (control r = 0.63, morpholino r = 0.84) (Supplemental Fig. 2). This demonstrates that although the hypergraph adjacency matrix does not consider directionality of association, it nevertheless captures the hierarchical structure of molecular cascades as inferred from random walks.

Early-life grb10a knockdown severely disrupts coordination within larval gene clusters

To investigate the co-ordination of the whole transcriptome with age-associated gene expression, hypergraph models were constructed for the SC and KD datasets based on the gene sets identified in the cluster analysis (297 and 119 genes respectively) (Fig. 1aii). This approach has been used to elucidate clusters of causally associated genes with co-ordinated regulation [44] (Fig. 6a). Figure 6b, c demonstrates this clustering, where colour intensity represents the number of shared correlations between each gene pair. In SCs, genes segregated into three highly connected groups, associated with 5, 10–15, and 20–30 dpf (161/297, 50/297, and 86/297 genes) (Fig. 6b). In KD animals, however, genes segregated into only two groups, either 5–15 dpf or 20–30 dpf (31/119 and 65/119 genes) (Fig. 6c). Notably, the large group of co-ordinated interactions at 5 dpf was absent in the KD dataset, with a combined group of genes correlating with 5–15 dpf identified instead.

Fig. 6

The coordination of age associated genes is altered by grb10a KD, as assessed by hypergraph analysis (A) hypergraphs form clusters of genes with higher order interactions between them; this allows us to refine groups of causally associated elements from a broader set of targets. (B) hypergraph analysis of the identified age-associated genes formed three such clusters of genes in SC zebrafish, corresponding to [1] 5 dpf, [2] 10–15 dpf, and [3] 20–30 dpf. (C) Two clusters were defined in the KD data, corresponding to [1] 5–15 dpf and [2] 20–30 dpf. (D) and (E) Violin plots of connectivity (left) and entropy (right) in the 20–30 dpf cluster. The KD transcriptome was more connected with higher entropy compared to SC, suggesting a more diverse set of connections are made by age associated genes (Wilcoxon Rank sum **** p < 0.0001)

Grb10a knockdown is associated with increased connectivity and crosstalk between age related genes

To quantify these differences in the co-ordination of the transcriptome, two network topology parameters were used: connectivity and entropy (Fig. 1aii). Hypergraph connectivity quantifies the number of higher order interactions within the transcriptome and entropy is a measure of ‘disorder’, as described above, with both related to function [43, 44].

Changes in connectivity (Fig. 6d) and entropy (Fig. 6e) were calculated for genes identified as active at 20–30 dpf by the previously described analyses. The impact of KD on the co-ordination of age-associated genes in the network was identified by comparing the connectivity and entropy of SC and KD datasets. To provide a comparison between age-associated genes and genes associated with other functions, connectivity and entropy were also calculated for randomly selected gene sets (Supplemental Fig. 1).

Genes associated with 20–30 dpf were more highly connected and less entropic (more ordered) than random genes in both the KD and SC datasets (p < 2.2 × 10–16). Comparison of age associated genes in KD and SC revealed a higher connectivity (1.20-fold, p < 2.2 × 10–16) (Fig. 6d) and a higher entropy (1.05-fold, p < 2.2 × 10–16) in the KD (Fig. 6e) This suggests age-associated genes in the KD zebrafish transcriptome have more interactions and more crosstalk (less ordered) than in the SC.

Pathways associated with transcriptome-wide remodelling support an alteration in cardiometabolic phenotype

Having identified a central set of age-related genes, a broader set of genes, which coordinate with these genes, represent upstream or downstream interactions associated with this central core of the hypergraph (Fig. 1aii). This broader set of genes are termed “peripheral genes” as they represent changes in expression consequential but distant from the hypergraph cluster (Fig. 7a). 12775 (KD) and 460 (SC) genes were defined in this set from the 20 to 30 dpf hypergraph clusters described previously, all significantly associated with age (rank-regression: KD q < 1.50 × 10–6, SC q < 4.44 × 10–2; Supplementary Table 4). 28 times more peripheral genes were implicated in the KD. Both datasets showed a skew towards a positive association with age (68% SC vs 59% KD) (Fig. 7b, Fig. 7c). The overlap in gene expression between the two datasets was 244 (77.9% of SC) and 60 (41.1% of SC) in the positively and negatively correlated sets, respectively. Thus, genes in the wider transcriptome with co-ordinated expression at 30 dpf in SC zebrafish demonstrated a similar pattern of expression in KD (Fig. 7d vs. e). However, the inverse is not true, as genes with co-ordinated expression at 15–30 dpf in the KD zebrafish did not show the same pattern in the SC (Fig. 7e vs. d).

Fig. 7

Analysis of the set of genes in the wider transcriptome shows a 27.8-fold increase in the KD ZF. (A) Identification of peripheral elements underpinning the hypergraph implicates a wider set of genes than the initial targets. A summary of those peripheral elements is presented here. (B, C) Venn diagrams of genes positively (B) and negatively (C) correlating with age. (D, E) Hierarchically clustered heat maps of gene expression of the genes identified in the wider transcriptome. Gene expression in the standard control (D) cluster into two age related groups, whereas expression in the knockdown (E) show significant dysregulation. (F, G) Gene set enrichment analysis (GSEA) ranked by R-value of rank age regression in the standard control (F) and knockdown (G). (H) GSEA ranked by R-value of rank age regression of the cluster of genes in the white box in (E)

To assign biological function, GSEA on the peripheral genes was performed as previously described (ranked by R-value, top 100 pathways with weighted set cover) (Fig. 1aiii). Results of this ontology analysis are outlined in Fig. 7f (SC) and Fig. 7g (KD). The two datasets featured distinctly different regulated pathways. RNA processing, a variety of metabolic pathways, and cardiovascular development featured in the KD dataset, whereas growth and immune signalling were featured in the SC dataset.

A distinct pattern of gene expression was identified in the group of genes dysregulated at 20–30 dpf in the KD dataset (white box, Fig. 7e). These genes were upregulated at 15 dpf, downregulated at 20 dpf, and re-upregulated at 30 dpf. This subgroup of 3460 genes (Supplementary Table 5) was associated with an enrichment of gene ontologies including lipid metabolism, muscle development and cardiovascular system development (Fig. 7h) and was synchronous with the spike in growth identified in Fig. 8a. Specific pathways following this pattern of expression included cardiovascular system development, muscle structure development, and developmental maturation.

Fig. 8

Growth in juvenile Zebrafish and adult size and skeletal muscle thickness at 18 months. (A) Mean total body length of SC and KD zebrafish up to 30 dpf. Following the embryonic growth spurt, KD zebrafish experienced an additional period of rapid growth between 15 and 20 dpf (unpaired t-test, **** p < 0.0001). Typical morphology for body length at 18 dpf is included for reference. (B) Individual total body length (left panel), mass (middle panel), and condition factor scores (right panel) for 18-month KD and SC zebrafish (n = 21–24). Length and mass were significantly higher in the KD (unpaired t-test, **** p < 0.0001, ** p = 0.005), while condition factor was significantly lower (* p = 0.02), indicating KD zebrafish have leaner bodies. (C) Thickness of 10 skeletal muscle fibres per fish (n = 5) stained with Masson’s Trichrome. Fibre thickness was marginally higher in KD compared to SC zebrafish (repeated measure ANOVA, p = 0.06). Data are presented as overall group mean and as mean ± SEM for each fish

The impact of early life grb10a KD persists into adulthood

As grb10a KD significantly impacted early-life growth, metabolism, and cardiovascular development and was associated with remodelling of the transcriptome between 5 and 30 dpf, investigation was conducted into the phenotypic differences in older zebrafish.

Total body length measurements up to 30 dpf are shown in Fig. 8a. Body size was significantly elevated in the KDs compared to SCs during larval development between 15 and 20 dpf (30% increase, 7.49 ± 0.27 mm vs 4.32 ± 0.16 mm, p < 0.0001, n = 10). This corresponded to the cluster of dysregulated genes identified by hypergraph modelling (Fig. 7h and white box in Fig. 7e). The growth profile of the KDs was shifted to an earlier age compared to the SC, suggesting a faster rate of maturation. The GSEA highlighted a reduction in the activity of developmental maturation pathways (normalised enrichment score < 1) in SCs vs. KDs, which was reflected in the growth rate at approximately 20 dpf. This time frame is also characterised by rapid morphological changes as the dorsal, anal, and pelvic fin buds develop and mature and the larvae begin to resemble an adult zebrafish [94].

Later-life body length and mass measurements were recorded in 18-month-old zebrafish to investigate the lasting impact of grb10a KD on the phenotype. Final body length (Fig. 8b.) was higher, with KDs 0.17 cm longer than SCs (32 ± 0.03 cm vs 30 ± 0.02 cm, p < 0.0001). KD fish were also approximately 10% heavier (0.54 ± 0.01 mg vs 0.49 ± 0.01 mg, p = 0.005). Fulton’s condition factor (an indicator of body condition) was lower in the KDs (1.64 ± 0.03 vs 1.75 ± 0.03, p = 0.022) despite the increase in mass, suggesting a “leaner” phenotype. To equate this to differences in body composition, skeletal muscle fibres, isolated from the base of the dorsal fin, were sectioned and stained with Masson’s Trichrome (Fig. 8c). The average muscle fibre diameter was approximately 20% greater in the KDs, which was marginally different between groups (p = 0.06).

The impact on the heart was assessed by qPCR and histology. As shown in Fig. 9a, myl7 expression (an index of muscle mass and hypertrophy) in the heart was over 20% greater in KDs (p < 0.0001, n = 3), while nppa expression (activated in response to ventricular stress during hypertrophy and heart failure [95, 96]) was reduced by approximately 40% (p = 0.0012, n = 3). There was no difference in pcna expression (proliferating cell nuclear antigen), suggesting there was no difference in proliferation in the cardiac tissue. These observations were confirmed in the age-related gene expression from the transcriptomic data (Supplementary Table 6). To support these findings, ventricular morphology was also assessed. As shown in Fig. 9b, the ratio of compact myocardial layer to trabeculated was significantly greater in the KD zebrafish (37% vs 19%, p = 0.03). This was coupled with an overall increase in tissue density (76% vs 64%, p = 0.03). The cross-sectional area of cardiac muscle myocyte bundles was also significantly greater in KD zebrafish compared to SC counterparts, with an approximate 50% increase in bundle cross sectional area (214 ± 24 µm2 vs 130 ± 18 µm2, p = 0.0195).

Fig. 9

Cardiac gene expression and Cardiometabolic function in 18-month-old adult Zebrafish. (A) qPCR results of three genes associated with cardiac performance in adult cardiac tissue, relative to β-actin. Myl7 expression was significantly elevated in KD zebrafish (unpaired t-test, **** p < 0.0001) while nppa expression was significantly down regulated (*** p = 0.001) compared to SC zebrafish. There was no significant difference between the expression of pcna (p = 0.3). (B) Ventricular morphometrics obtained by Masson’s Trichrome histology comparing KD and SC compacta thickness, tissue density, and fiber cross sectional area. The compacta layer was significantly thicker (unpaired t-test, * p = 0.03) and overall tissue density was higher (* p = 0.03) in the KD zebrafish, and the fiber cross-sectional area (CSA) was increased (* p = 0.02) (n = 5–6). (C) Maximum (MMR) and Basal (BMR) oxygen uptake rate of adult (18 month) zebrafish, adjusted for body mass. BMR was comparable between the two groups (unpaired t-test, NS), while MMR was greater in the KD zebrafish (**** p < 0.0001), resulting in a greater aerobic scope (dotted line, unpaired t-test, *** p = 0.0007). (D) Fasting Glucose concentrations (left panel), glucose tolerance testing (centre panel) and insulin sensitivity testing (right panel) in adult (18 month) KD and SC zebrafish. Fasting glucose was higher in the KD (n = 16) than in the SC (n = 21) (unpaired t-test, p = 0.01). There were marginally higher glucose concentrations in the KD zebrafish in the glucose tolerance test (repeated measures ANOVA, p = 0.07) but both KD and SC zebrafish responded similarly to insulin (repeated measures ANOVA, NS). Treatment started immediately after the first measurement. Data are presented as mean values ± SEM

To determine whether metabolic rate was elevated later in life, oxygen uptake in adult (18-month) zebrafish was also investigated. Maximum oxygen uptake, achieved following exhaustive activity, was elevated by ~ 25% in KDs, as shown in Fig. 9c (1645 ± 90 mgO2kg−1 h−1 vs 1309 ± 27 mgO2kg−1 h−1, vs p < 0.0001). There was no significant difference in basal metabolic rate (486 ± 47 mgO2kg−1 h−1 vs 474 ± 28 mgO2kg−1 h−1). Consequently, aerobic scope was greater in the KDs (1158 ± 101 mgO2kg−1 h−1 vs 810 ± 43 mgO2kg−1 h−1, p = 0.0007).

As grb10a is involved in regulating the insulin signalling pathway, glucose tolerance and insulin sensitivity tests were performed on adult zebrafish (18 months) to determine whether there was a lasting biological impact. Fasting blood glucose (Fig. 9d) showed a significant difference between SC and KD fish, higher in KDs (10.2 ± 0.9 mmol/L vs 6.0 ± 0.5 mmol/L, p = 0.01). Glucose tolerance testing showed marginally higher glucose levels in the KD zebrafish at p = 0.07, while both groups produced a similar response to insulin administration.

Cardiometabolic genetic risk is functionally linked to the grb10a KD age-related transcriptome

There were 29 shared age-related genes that overlapped between SC (297 genes) and KD (119 genes) hypergraphs (Fig. 6). A total of 19 variants were identified to establish a group of 6 genes we termed ‘cardio-metabolic disease GWAS genes’ (Supplementary Table 7). Recently, 170 serum protein quantitative trait loci have been identified for cardiometabolic trait loci [97] representing 0.89% of the protein coding genome. The observation of 6 cardio-metabolic disease genes out of 29 represents a highly significant enrichment (Fisher’s Exact Test p-value = 0.0006 [2-tailed], odds ratio = 29).

Ranked row sums of the hypergraph incidence matrix (representing the hypergraph connectivity) were compared for each GWAS gene in SC and KD age-related transcriptome hypergraphs (Fig. 10A). Four GWAS genes, APOA1, APOB, BLK and SERPINH1, were ranked lower in the KD hypergraph compared to the control hypergraph. This suggests these cardiometabolic disease associated genes are less functionally important after the transcriptome changes induced by grb10a knockdown in early life. Other cardiovascular and metabolic GWAS genes, MMP3 and MMP12, did not show a difference in rank between the control and knockdown model, suggesting the transient grb10a knockdown does not affect the functional importance of these genes in the transcriptomic network.

Fig. 10

Genes associated with cardio-metabolic diseases. Cardio-metabolic disease GWAS genes ranked by normalized row sum in the age-related transcriptome hypergraphs of control and GRB10 knockdown zebrafish. Six genes out of 29 differentially expressed genes were present in both the SC and KD were shown to be related to cardiometabolic health. A) The impact of these genes on hypergraph connectivity and, B) higher order interactions in the transcriptome hypergraph were measured and ranked using a normalized ranking score, higher value indicates greater number of interactions in the transcriptome and implies greater impact on function by increased cross talk between biological pathways. C) Human Genome Wide Association Study (GWAS) relationships identified by phenome wide association using the fourteen human orthologous genes present in overlap between SC and KD hypergraphs. Identified using PHEWAS catalogue. D) Interactome model of the fourteen human orthologous genes present in overlap between SC and KD hypergraphs, demonstrating coherency and implying co-ordination

We assessed the impact of higher order interactions in the transcriptomic hypergraph independently using the reduced adjacency matrix row sums. This was performed for each of the identified cardio-metabolic disease GWAS genes (Fig. 10B). This analysis demonstrated persistent differences in higher order interactions between KD and SC with four of the genes more highly ranked in KD (BLK, MMP12, MMP3 and SERPINH1) and two of the genes more highly ranked in SC (APOA1 and APOB).

Permutation of both the hypergraph incidence matrix and hypergraph reduced adjacency matrix was performed 1000 times and used to show robustness of the network ranks associated with the six human orthologous genes (Supplementary Fig. 3).

To examine wider genomic relationships, we examined phenome wide association for the 14 human orthologous genes that mapped to the 29 shared age-related genes that overlapped between SC and KD hypergraphs (Supplementary Table 8). This analysis highlighted a range of related genome wide associations studies (Fig. 10C) that mapped mainly to pathways related to cholesterol metabolism and lung function (7.4 × 10–10 < p < 0.001, Supplementary Table 8). An inferred interactome model of these 14 human genes demonstrated coherence of transcriptomic action (Fig. 10D) and highlighted overall impact on extracellular matrix, protein-lipid complex and collagen pathways (Supplementary Table 9).

We demonstrated by measuring phosphorylated protein that AKT and S6 were significantly elevated in grb10a KD zebrafish indicative of modulation of IGF and Insulin pathways (Fig. 2C). We screened the SC and the KD hypergraphs (Fig. 6) for evidence of phospho-ERK activity by mapping the associated genes to ERK targets [