A Mouse Systems Genetics Approach Reveals Common and Uncommon Genetic Modifiers of Hepatic Lysosomal Enzyme Activities and Glycosphingolipids

Identification of genetic modulators of lysosomal enzyme activities and glycosphingolipids (GSLs) may facilitate the development of therapeutics for diseases in which they participate, including Lysosomal Storage Disorders (LSDs). To this end, we used a systems genetics approach: we measured 11 hepatic lysosomal enzymes and many of their natural substrates (GSLs), followed by modifier gene mapping by GWAS and transcriptomics associations in a panel of inbred strains. Unexpectedly, most GSLs showed no association between their levels and the enzyme activity that catabolizes them. Genomic mapping identified 30 shared predicted modifier genes between the enzymes and GSLs, which are clustered in three pathways and are associated with other diseases. Surprisingly, they are regulated by ten common transcription factors, and their majority by miRNA-340p. In conclusion, we have identified novel regulators of GSL metabolism, which may serve as therapeutic targets for LSDs and may suggest the involvement of GSL metabolism in other pathologies.


Introduction
Hydrolytic enzymes are abundant in lysosomes [1]. In a healthy cell, the biosynthesis and catabolism of macromolecules are subject to regulatory mechanisms that maintain cellular homeostasis [2]. The degradative processes in lysosomes are controlled by their own enzymes [3,4]. Lysosomes play a central role in several biological processes, including energy metabolism, signaling, plasma membrane repair, secretion, and others [3]. Loss-offunction variants in genes encoding lysosomal proteins cause lysosomal storage disorders (LSDs), a group of diseases characterized by intracellular buildup of partially degraded material [5]. Growing evidence suggests that variants in lysosomal genes increase the risk of developing Parkinson's disease (PD) [6,7].
In the sphingolipidoses, a subset of LSDs, glycosphingolipids (GSLs) accumulate in late endocytic organelles (late endosomes/lysosomes) and participate in their pathological cascades [8]. Current treatments for LSDs include substrate reduction therapy (SRT), which aims to reduce the rate of biosynthesis of stored substrates [5,9,10], and enzyme replacement therapies (ERT) aimed at replacing a deficient enzyme [11,12]. Emerging treatments include gene and cell therapies [13][14][15] and chaperones for improving enzyme folding and trafficking [16]. Although there is a range of therapeutic options for LSDs, they have limitations, such as tissue accessibility [17], antibody-mediated reaction [18], cost [19], and others. So far, therapies aimed at increasing enzyme activity or reducing lipid levels by modulating a second (modifier) gene have not been studied. In this context, a deeper understanding of the regulatory mechanisms that govern GSLs metabolism must be uncovered to fully develop this approximation.
Genome-wide association studies (GWAS) in humans and systems genetics strategies, which include gene mapping in model organisms, have identified genetic regulators of physiological and pathophysiological processes [20][21][22]. The Hybrid Mouse Diversity Panel (HMDP) has been a useful tool because genomes and tissue transcriptomes are freely available, allowing the combination of modifier gene mapping by GWAS and pathway analysis [23,24]. In this study, we have analyzed the activities of 11 lysosomal enzymes and several of their natural substrates in 25 strains of the HMDP panel followed by gene mapping and transcript integration. We identified a lack of correlation between most enzyme activities and their mRNA levels. Similarly, most substrates had no association between their levels and the enzyme activity that catabolizes them. Finally, we mapped putative modifier genes of each lysosomal enzyme and GSL by GWAS. We found associations between the mRNA levels of many modifier genes and enzyme activities or GSL levels. We clustered the putative modifiers in pathways and identified common and uncommon genetic regulators between GSLs and lysosomal enzymes, including transcription factors that regulate them. Our discoveries may help develop novel therapeutics for diseases with altered lysosomal enzyme activities and GSLs.

High Variability in the Hepatic Activity of Lysosomal Enzymes across Mouse Strains
We measured hepatic enzyme activity of β-hexosaminidase A and B (defective in Tay-Sachs and Sandhoff disease, respectively), α-neuraminidase (defective in Sialidosis/Mucolipidosis Type I), α-galactosidase A and B (defective in Fabry and Schindler disease), β-D-galactosidase (defective in GM1 Gangliosidosis), α-glucosidase (defective in Pompe), chitotriosidase (elevated in Gaucher disease), α-L-fucosidase (defective in fucosidosis), lysosomal acid phosphatase (elevated in patients with Gaucher), and Tartrateresistant acid phosphatase (TRAP; altered in Gaucher disease) by fluorimetry in liver samples derived from 25 inbred mice strains using 4-methylumbelliferone (4-MU) based artificial substrates. We observed significant variability in the average enzymatic activity between the different strains (ANOVA p ≤ 0.05) ( Figure 1). We did not find changes in α-galactosidase A, lysosomal acid phosphatase, and TRAP activities across the tissues analyzed ( Figure 1d,j,k). We observed unique activity distribution patterns across the strains for the other enzymes, suggesting specific modifiers for each enzyme.

Lack of Correlation between the Enzyme Activity and Its mRNA Levels
Advantages of using tissues derived from the HMDP panel of inbred mouse strains include the fact that their genomes are sequenced, and transcriptomic data are available. Thus, we analyzed potential correlations between the genes encoding lysosomal enzymes and their activities. Recently we described the natural variation of hepatic acid β-glucocerebrosidase levels across many different mouse strains and included them in this analysis [20]. We did not identify significant correlations between enzyme activity and its transcript levels ( Figure 2), with the only exception being Glb1, the gene encoding for β-D-

Lack of Correlation between the Enzyme Activity and Its mRNA Levels
Advantages of using tissues derived from the HMDP panel of inbred mouse strains include the fact that their genomes are sequenced, and transcriptomic data are available. Thus, we analyzed potential correlations between the genes encoding lysosomal enzymes and their activities. Recently we described the natural variation of hepatic acid β-glucocerebrosidase levels across many different mouse strains and included them in this analysis [20]. We did not identify significant correlations between enzyme activity and its transcript levels ( Figure 2), with the only exception being Glb1, the gene encoding for β-D-galactosidase (r = 0.5775; p ≤ 0.002) (Figure 2c). These results indicate that mRNA levels are a poor proxy for enzyme activities.
galactosidase (r = 0.5775; p ≤ 0.002) (Figure 2c). These results indicate that mRNA levels are a poor proxy for enzyme activities.
(h) Acid-β-glucosidase. (i) Lysosomal acid phosphatase. (j) Tartrate-resistant acid phosphatase. The Pearson's correlation was performed using 23 strains of mice. Enzyme activities are expressed as nmol/mg/hour and mRNA levels, which were downloaded from repository GSE16780 UCLA Hybrid MDP Liver Affy HT M430A [24] are expressed as log2 transformed. r, correlation; p, p-value.

High Variability in the Hepatic Glycosphingolipid Levels across Mouse Strains
Next, we measured the levels of GSLs in livers of the inbred mice strains in which we had access to enough material for three biological replicates (23/25) by Normal Phase-High-Performance Liquid Chromatography (NP-HPLC). We observed significant variability in GSLs among the strains, especially in total GSLs, GM3-Gc, GM2-Gc, GM1agc, GM3, Gb3, GM1a, GM1b, GD1b, and GD1a ( Figure 3). For example, the levels of GM3-Gc were significantly increased (ANOVA p < 0.0001) in NOD/ShiLtJ compared with the other samples ( Figure 3b). These results indicate that GSLs levels vary across strains.
(h) Acid-β-glucosidase. (i) Lysosomal acid phosphatase. (j) Tartrate-resistant acid phosphatase. The Pearson's correlation was performed using 23 strains of mice. Enzyme activities are expressed as nmol/mg/hour and mRNA levels, which were downloaded from repository GSE16780 UCLA Hybrid MDP Liver Affy HT M430A [24] are expressed as log2 transformed. r, correlation; p, p-value.

High Variability in the Hepatic Glycosphingolipid Levels across Mouse Strains
Next, we measured the levels of GSLs in livers of the inbred mice strains in which we had access to enough material for three biological replicates (23/25) by Normal Phase-High-Performance Liquid Chromatography (NP-HPLC). We observed significant variability in GSLs among the strains, especially in total GSLs, GM3-Gc, GM2-Gc, GM1agc, GM3, Gb3, GM1a, GM1b, GD1b, and GD1a ( Figure 3). For example, the levels of GM3-Gc were significantly increased (ANOVA p < 0.0001) in NOD/ShiLtJ compared with the other samples ( Figure 3b). These results indicate that GSLs levels vary across strains.

Correlations between the GSLs and the mRNA Levels of the Biosynthetic Genes
A possibility is that GSL levels could correlate with their biosynthesis rate. Since we started from frozen tissues, we could not test this directly. Instead, we utilized the transcriptomic data available from the repository GSE16780 UCLA Hybrid MDP Liver Affy HT M430A [24]. We found transcript probes for 21 mRNA of the 21 anabolic enzymes of the GSLs pathway and four GSL transfer proteins. The analyzed gene list of the biosynthetic pathway is presented in the Supplementary Table S1. The expression values were organized according to GSLs levels from lowest to highest and presented as a heatmap.

Correlations between the GSLs and the mRNA Levels of the Biosynthetic Genes
A possibility is that GSL levels could correlate with their biosynthesis rate. Since we started from frozen tissues, we could not test this directly. Instead, we utilized the transcriptomic data available from the repository GSE16780 UCLA Hybrid MDP Liver Affy HT M430A [24]. We found transcript probes for 21 mRNA of the 21 anabolic en-zymes of the GSLs pathway and four GSL transfer proteins. The analyzed gene list of the biosynthetic pathway is presented in the Supplementary Table S1. The expression values were organized according to GSLs levels from lowest to highest and presented as a heatmap. The analysis showed significant correlations for Cgt (r = −0.4263; p = 0.042) with total GSLs (Figure 4a Table S2); thus, we analyzed potential correlations between GSL levels and the enzyme activity that catabolizes them across the mouse panel.  Table S2); thus, we analyzed potential correlations between GSL levels and the enzyme activity that catabolizes them across the mouse panel.

Lack of Correlation between Hepatic Lysosomal Enzyme Activity and Their Natural Substrates across Mouse Strains
It is possible to speculate that the strains that present high activity of a particular enzyme should have reduced levels of its natural substrate because the enzyme catabolizes it. Unexpectedly, for most enzymes, we did not find significant correlations between the GSL levels and the enzyme activity that degrades it ( Figure 5), except for neuraminidase and GM3-Gc (r = −0.4706; p = 0.0234) (Figure 5g). These results suggest that for most strains, the rate of biosynthesis and/or uptake of GSLs varies along with the catabolic rates which most likely are genetically regulated.

Identification of Putative Modifier Genes of Lysosomal Enzyme Activity and Sphingolipids Levels
To identify genetic regulators, we conducted genome-wide association studies with a quality control analysis that considered the population structure of the HMDP panel

Lack of Correlation between Hepatic Lysosomal Enzyme Activity and Their Natural Substrates across Mouse Strains
It is possible to speculate that the strains that present high activity of a particular enzyme should have reduced levels of its natural substrate because the enzyme catabolizes it. Unexpectedly, for most enzymes, we did not find significant correlations between the GSL levels and the enzyme activity that degrades it ( Figure 5), except for neuraminidase and GM3-Gc (r = −0.4706; p = 0.0234) (Figure 5g). These results suggest that for most strains, the rate of biosynthesis and/or uptake of GSLs varies along with the catabolic rates which most likely are genetically regulated.

Identification of Putative Modifier Genes of Lysosomal Enzyme Activity and Sphingolipids Levels
To identify genetic regulators, we conducted genome-wide association studies with a quality control analysis that considered the population structure of the HMDP panel strains to reduce false associations [25,26]. We used enzyme activity levels as a trait and included the β-glucosidase activity, which we reported previously in the same and a few other strains [20]. For all the enzymes together, we identified 211 significant Single Nucleotide Variants (SNVs) that passed the empiric threshold of significance p ≤ 4.1 × 10 −6 (−log10P = 5.39), previously calculated by permutations [21,22,26], while the Bonferroni threshold was p ≤ 3.9 × 10 −7 [26]. These SNVs were located in different genomic regions (exonic, intronic, UTR3, downstream, and intergenic) ( Table 1, Supplementary Table S3) in a total of 137 non-redundant genes. Similarly, we identified 3215 SNVs associated with GSLs levels (1744 non-redundant genes) whose variants are located in different genomic regions ( Table 1, Supplementary Table S3). These analyses indicated that our strategy has sufficient power to map putative modifier genes.

Non-Redundant Genes
lysosomal enzymes

Correlations between the Traits and the mRNA Levels of Putative Modifiers
To prioritize the putative modifier genes that could regulate each enzyme, we searched for correlations between the transcript levels of putative modifier genes and their traits (enzyme activity and GSL levels, respectively) ( Figure 6). We found transcript probes for 67 mRNA of the 137 putative modifiers of the enzymes. The expression values were organized according to enzyme activity from lowest to highest and presented as a heatmap. The analysis showed significant correlations in Fip1l1 (r = −0.4462; p = 0.0254) with α-Lfucosidase (Figure 6a). For β-D-galactosidase with Lyplal1 (r = −0.702; p = <0.0001), Arrdc4 (r = 0.627; p = 0.0008), Pde2a (r = 0.5306; p = 0.0064), Glb1 (r = 0.5753; p = 0.0026), Bptf (r = 0.5135; p = 0.0087), Oxr1 (r = −0.447; p = 0.0251) (Figure 6b). No significant correlations were found for the other enzymes analyzed. We used SIFT to explore the impact of genetic variants on the genes identified by GWAS (benign or deleterious changes) associated with changes in enzyme activity [27], because the full genomes of the strains are known [28]. This strategy identified 308 predicted deleterious variants (Supplementary Table S4) in 43 of the 67 genes whose functions are related to organelle biogenesis (Chchd6) [29], intracellular signaling (Pde4dip) [30], and tissue development (Fam181b) [31], among others. These results suggest that amino acid substitution could affect protein function and signaling pathways leading to changes in enzyme activity.

Enrichment Analysis and Common Modifier Genes between Glycosphingolipids Levels and Lysosomal Enzyme Activities
If there is an orchestrated regulation of GSL levels and the enzymes that degrade them, it would be expected to observe enrichment in common pathways [34]. We therefore utilized gProfiler [35] to perform enrichment analysis using the putative modifier genes lists. For the modifier of enzyme activities, we found significantly associated pathways such as cell periphery (p = 5.9 × 10 −4 ), plasma membrane (p = 2.4 × 10 −3 ), and integral components of the plasma membrane (p = 2.6 × 10 −2 ) (Figure 7b), which could be related to endocytic processes necessary to deliver key molecules to the lysosome, including the lysosomal enzymes that can be recycled from the extracellular space. Significant biological processes analysis included regulation of cellular processes (p = 3.9 × 10 −2 ) (Figure 7d) (Supplementary Table S5). We did not find significant enrichment for the molecular function category. For GSLs, we observed enrichment in terms like cytoplasm (p = 3.  Table S5). The same analysis was performed to identify putative modifiers of GSL levels (Figure 6c-f). For 1744 non-redundant SNVs, we found expression values for 994 genes. The analysis identified 45 significant correlations, of which 33 were correlated with GM3-Gc levels, 10 genes with LacCer, and one gene with GD1b and GA2 (Figure 6c-f). Overall, we recorded 4.9% (52/1061) of significant correlations distributed between the two traits. We also explored the impact of genetic variants associated with changes in GSLs with SIFT [27]. This strategy identified 515 deleterious variants predicted to disrupt the protein structure (Supplementary Table S4) in 132 genes related to DNA methyltransferase activity (Setdb1) [32] and synapse (Slitrk1) [33], among others.

Enrichment Analysis and Common Modifier Genes between Glycosphingolipids Levels and Lysosomal Enzyme Activities
If there is an orchestrated regulation of GSL levels and the enzymes that degrade them, it would be expected to observe enrichment in common pathways [34]. We therefore utilized gProfiler [35] to perform enrichment analysis using the putative modifier genes lists. For the modifier of enzyme activities, we found significantly associated pathways such as cell periphery (p = 5.9 × 10 −4 ), plasma membrane (p = 2.4 × 10 −3 ), and integral components of the plasma membrane (p = 2.6 × 10 −2 ) (Figure 7b), which could be related to endocytic processes necessary to deliver key molecules to the lysosome, including the lysosomal enzymes that can be recycled from the extracellular space. Significant biological processes analysis included regulation of cellular processes (p = 3.9 × 10 −2 ) (Figure 7d) (Supplementary Table S5). We did not find significant enrichment for the molecular function category. For GSLs, we observed enrichment in terms like cytoplasm  Table S5). Biological processes terms revealed 328 pathways, including system development (p = 4.3 × 10 −39 ), anatomical structure development (5.6 × 10 −38 ), and multicellular organism development (p = 1.4 × 10 −37 ). We searched for the overlap between the cellular component domains of modifiers of enzyme activity and GSLs, which resulted in three common pathways (GO:0071944-cell periphery, GO:0005886-plasma membrane, and GO:0005887-integral component of plasma membrane) (Figure 7c) and one pathway associated with biological processes (GO:0050794; regulation of cellular process) (Supplementary Table S5).
To better understand the molecular regulation of these 30 genes, we analyzed the
Others participate in neurodegenerative conditions; PD, schizophrenia, and intellectual disability (Tenm4, Pde4dip, Grid2, Arhgap18) [51,52]. These results suggest that lysosomal enzymes and GSLs may play a role in their pathophysiology and should be explored further ( Table 2).
To better understand the molecular regulation of these 30 genes, we analyzed the transcription factors that bind to their promoters and/or enhancers (Figure 7f). We found no information for three of the 30 genes since they are putative (Rik) genes. The following transcription factors can bind to the 27 genes for which we have information: REST, TBP, CEBPB, EP300, POLR2A, FOS, DPF2, CTCF, RAD21, and SP1. Some of these transcription factors are broad regulators of transcription, such as TBP and POLR2A, while others are selective for specific processes, such as CTCF and RAD21. Considering all the promoters/enhancers of the 27 shared genes, we identified a total of 533 transcription factors that can bind them, although some only bind a few genes (Supplementary Table S6). We also searched for potential shared microRNA (miRNA) regulators using miRTarBase, a curated microRNA database [53]. We identified that miR-340-5p can bind to 11 of the 27 known common genes (Tusc1, Fam91a1, Zc3h12c, Adamts5, Tmem135, Tenm4, Grid2, Csnk1g3, Cdh6, Fam181b, and Pde4dip; p = 2.2 × 10 −2 ) (Figure 7g). This result suggests that miRNA-340-5p regulates GSLs metabolism and may be involved in the pathogenesis of LSDs and the disorders described in Table 2.

Discussion
In this study we searched for genetic modulators involved in the regulation of the lysosomal enzyme activities and the levels of substrates related to GSLs, with the idea of finding novel therapeutics targets for disorders in which they participate. By GWASs, we identified common and uncommon genetic regulators, evaluated the associations between modifier gene mRNA levels and each trait, and also clustered them in pathways. We identified 30 shared putative modifiers and described the transcription factors that are predicted to regulate them, and we noted that the miRNA340-5p can bind to 11 of these genes.
Our first unexpected finding was that most lysosomal enzyme activities do not correlate with their mRNA levels, nor with most of their substrate levels. Although enzyme activity can decrease with age [92], we used sex and age-matched samples; thus, the variation observed across strains was shown not to be due to any of these factors.
Another unexpected finding was that GM2-Gc levels correlate with the mRNA levels of the Cgt gene, which encodes for the UDP-galactose ceramide galactosyltransferase (CGT). CGT is a key enzyme for the biosynthesis of galactocerebrosides. Gangliosides, including GM2 derivates, are built from glucosylceramide and not from the galacto series [93]. However, for most of the biosynthetic genes there were no associations between the amount of lipids and the transcript levels of their anabolic pathways. Altogether, our results suggest that the GSL biosynthesis rate and uptake differ across the mouse strains, suggesting the existence of specific modifier genes for each trait.
Our third unexpected finding was that TFEB, the master transcriptional regulator of lysosomal genes [94], did not appear in the list of modifiers of lysosomal enzymes. This may be due to the fact that we screened for enzymatic activity instead of mRNA levels, and we showed a lack of correlation between transcript levels and enzyme activity under physiological conditions, at least for most enzymes. One exception was β-D-galactosidase, for which we found a positive correlation between its transcript levels and activity. Furthermore, the GWAS for this enzyme identified Glb1, the gene encoding for β-D-galactosidase, as a putative modifier of its activity, validating the power of discovery of our population-based strategy [95].
Our study had some limitations: First, we quantified lysosomal traits from liver homogenates that were not in living or isolated organelles, which may have diluted enzyme activity or promoted molecular interactions that might not occur in vivo because of cellular compartmentalization. Second, we could not directly measure GSLs biosynthesis and uptake because we started with mouse liver samples. Third, we used SNV catalogs with imputation, which may lead to false associations, though with increased mapping resolution.
Most of the enzymes we assayed are associated with LSDs [5,8]. For many LSDs, no therapies are available, and the few currently available treatments have severe limitations [5]. In this context, targeting a modifier gene could be a novel therapeutic approach. For example, lack of β-D-galactosidase activity triggers GM1 gangliosidosis, a disease with no approved therapies [96]. Our study identified the druggable Lypla1 and Pkm genes as putative modifiers of β-D-galactosidase activity, which can be pharmacologically modulated [97,98]. We found other druggable genes as well for several traits, and with the current gene editing technologies virtually any gene can be targeted. The potential modifying effects of these genes and compounds can be tested in LSDs disease models.
A hallmark of the sphingolipidoses is the intracellular buildup of GSLs, so strategies aimed at reducing their levels could lead us to novel therapies [5,8]. GSLs comprise a ceramide moiety with one or more sugar residues linked to it [99]. An approved therapy for Gaucher and Niemann-Pick disease type C is Miglustat [100,101], a small molecule inhibitor of GSL biosynthesis, thus reducing their levels. Our GSLs GWAS identified more than 50 genes previously associated with sphingolipid metabolism, which served as a positive control, including B3gnt5, Cln8, Hexb, Pnpla1, St8sia1, and Cgt. B3gnt5 regulates GSLs metabolism and lung tumorigenesis [102]. Our study also identified Lipc as a modifier of GM3-Gc levels, which has been previously associated with elevated serum levels of liver enzymes (alkaline phosphatase and γ-glutamyl transferase) [103], suggesting a new connection between GM3-Gc and liver damage. Variants in LIPC, CPS1, PABPC4, CITED2, TRPS1, and MVK are associated with changes in plasma lipoprotein levels [104], connecting novel traits to GSLs metabolism.
Lysosomal leakage has been associated with Alzheimers' [105], cancer, and inflammation among other conditions [106]. Recently, the phosphoinositide signaling pathway was implicated in lysosomal repair [107]. Many genes of this pathway appear in our discovery list (Osbpl9, Osbpl6, Pde4dip, Pde2a, Pde1a, Pde7a, Pde7b, Pde4d, Pde8b, Pld5, Pik3r1, Pip4k2a, Pip5k1a, Pip5k1b, Pi4kb, Pdpk1, Atg4c, Atg10), suggesting that integrity of the lysosomal compartment is key to the proper functioning of enzymes and/or that these enzymes and lipids participate in lysosomal repair. Furthermore, this novel lysosomal repair pathway may facilitate the development of novel therapeutics for these diseases with lysosomal leakage.
In conclusion, we described putative regulators of hepatic lysosomal enzymes and GSLs, many of them druggable and associated with diseases where alterations in GSL metabolism have not been previously described and should be assessed. We expect our findings may facilitate the development of novel therapeutics for conditions with alterations in these traits.

Mouse Tissues
We used 8 weeks-old mice livers derived from 25 inbred mouse strains, which were kindly donated by Dr. Aldons Lusis (University of California, Los Angeles, CA, USA).  Table S7).

Glycosphingolipids Levels Quantification
The GSLs were extracted and measured by Normal Phase-High-Performance Liquid Chromatography (NP-HPLC) following published methods [145]. Briefly, the aqueous tissue extract was homogenized in chloroform/methanol (C:M) (1:2 v/v) and kept overnight at 4 • C. Then, the extracts mixture was centrifuged at 3000 rpm for 10 min at room temperature. We added 0.5 mL of PBS and 0.5 mL of chloroform to the supernatant followed by a 3000rpm centrifugation for 10 min at room temperature. The lower phase was carefully removed and dried under a stream of nitrogen gas (N2) in a heating block (42 • C), resuspended in 40 µL C:M 1:3 v/v and mixed with the upper phase. Afterwards, glycosphingolipidsderived oligosaccharides were purified from the samples using C18 columns (Telos, Kinesis, UK) previously pre-equilibrated with 1.25 mL methanol (four times) and 1.25 mL deionized water (three times). We loaded the mixed phase (lower/upper) onto a column and rinsed the sample tube with 1 × 1 mL of deionized water. Then, the C18 column was washed with 4 × 1.25 mL deionized water and eluted it with 1 × 1 mL (C:M) (98:2 v/v), 2 × 1 mL (C:M) (1:3 v/v), 1 × 1 mL methanol. The eluates were dried under N2 current and digested with a recombinant Endoglycoceramidase I (rEGCaseI) (GenScript, Oxford, UK) in buffer 50 mM sodium acetate, pH 5.0, 0.6% TritonX-100 (4 µL enzyme + 86 µL buffer) at 37 • C for 16 h. The released glycans were labeled with 310 µL of labelling mix (30 mg/mL anthranilic acid (2AA) and 45 mg/mL sodium cyanoborohydride) in 4% sodium acetate, 2% boric acid in methanol, and heated at 80 • C. Then, we cooled the samples and mixed them with 3 × 1 mL acetonitrile: deionized water (97:3) (v/v) and added them to a Discovery DPA-6S-SPE tube (Supelco, PA, USA), pre-equilibrated with 1 × 1 mL acetonitrile, 2 × 1 mL deionized water, and 3 × 1 mL acetonitrile. The columns were cleaned with 3 × 1 mL acetonitrile: deionized water (95:5) (v/v), and the tubes were washed with 2 × 1 mL acetonitrile: deionized water (95:5) (v/v) and eluted in 0.6 mL deionized water. We took 60 µL from 0.6 mL sample eluted, added 140 µL acetonitrile, and injected 50 µL of this mix (deionized water: acetonitrile) (30:70) (v/v) onto NP-HPLC (Waters Alliance 2695 separations module and multi-fluorescent detector set at Ex 360/Em 425 nm). To calculate molar quantities from peaks in the chromatogram, we included a calibration standard containing 2.5 pmol 2AA-labelled chitotriose (Ludger, Oxford, UK) for each NP-HPLC run [145]. The chromatographic data were processed using Waters Empower software 3 (Waters, Milford, MA, USA). Fluorescence values by sample were normalized to protein content using a BCA Assay kit (Merck KGaA, Darmstadt, Germany).

Genome-Wide Association Studies (GWAS)
We used the genotype of each strain, and the enzymatic activity or substrate as trait, and its kinship matrix to perform the GWAS using The Efficient Mixed Model Association (EMMA) v.1.1.230 in the R package [26,146]. We used PLINK to remove SNVs in linkage disequilibrium to avoid false associations [25], considering an R 2 = 0.25, leaving 127,285 independent variants out of the initial four million variants downloaded from the mouse HapMap reference panel (http://mouse.cs.ucla.edu/mousehapmap/full.html, accessed on 28 September 2020) [147].

Gene Expression Array and Heat Maps
For gene expression correlations, we obtained inbred mouse hepatic transcript data from the repository GSE16780 UCLA Hybrid MDP Liver Affy HTM430A [24]. The mRNA levels in the repository were expressed as log2 transformed and were calculated from the Affimetrix chip with the robust multiarray average (RMA) method. To plot the heatmaps, we used Morpheus software (https://software.broadinstitute.org/morpheus, accessed on 15 February 2022).

Enrichment Analysis
We used gProfiler [35] with the default settings to perform the pathway enrichment analyses.

Identification of Transcription Factors
We consulted the GeneHancer (GH) database, a catalogue of genome-wide enhancerto-gene and promoter-to-gene associations, through GeneCards ® (https://www.genecards. org/Guide/GeneCard, accessed on 6 September 2022) [148]. Only transcription factors with a significative GH Score were considered.

Statistics
We used Student's t-test, ANOVA with Bonferroni correction, and Pearson correlation. All tests were two-tailed. The significance was considered to be p < 0.05. We used an R package [146] and Prism v9.1.0 (GraphPad software, San Diego, CA, USA) for these analyses.

Informed Consent Statement: Not applicable.
Data Availability Statement: We obtained inbred mouse hepatic transcript data from the repository GSE16780 UCLA Hybrid MDP Liver Affy HTM430A reported in [24].