Global Transcriptomic Analysis of Zebrafish Glucagon Receptor Mutant Reveals Its Regulated Metabolic Network

The glucagon receptor (GCGR) is a G-protein-coupled receptor (GPCR) that mediates the activity of glucagon. Disruption of GCGR results in many metabolic alterations, including increased glucose tolerance, decreased adiposity, hypoglycemia, and pancreatic α-cell hyperplasia. To better understand the global transcriptomic changes resulting from GCGR deficiency, we performed whole-organism RNA sequencing analysis in wild type and gcgr-deficient zebrafish. We found that the expression of 1645 genes changes more than two-fold among mutants. Most of these genes are related to metabolism of carbohydrates, lipids, and amino acids. Genes related to fatty acid β-oxidation, amino acid catabolism, and ureagenesis are often downregulated. Among gcrgr-deficient zebrafish, we experimentally confirmed increases in lipid accumulation in the liver and whole-body glucose uptake, as well as a modest decrease in total amino acid content. These results provide new information about the global metabolic network that GCGR signaling regulates in addition to a better understanding of the receptor’s physiological functions.


Introduction
Glucagon, derived from pancreatic α-cells, activates glucagon receptor (GCGR) signaling to increase hepatic glucose production and maintain normal blood glucose levels during fasting [1]. GCGR is a seven-transmembrane class B G-protein-coupled receptor [2,3]. Stimulation of GCGR activates adenylyl cyclase and increases cAMP levels [4]. Elevated cAMP stimulates several pathways, resulting in an increase of gluconeogenesis, glycogenolysis, and fatty acid oxidation [5]. Inhibition or inactivation of GCGR is shown to lower blood glucose in preclinical models of type 1 and type 2 diabetes [6][7][8]. Consequently, GCGR antagonists and neutralizing antibodies have emerged as a class of novel therapeutics for treating diabetes [9][10][11][12]. Several antagonists and antibodies are undergoing different stages of clinical trials, but none have been approved for clinical use [13].
We previously generated GCGR knockout zebrafish that display hypoglycemia, hyperglucagonemia, and compensatory α-cell hyperplasia, demonstrating that the function of glucagon signaling is at least partially conserved between humans and zebrafish [25]. The amenability of zebrafish to live imaging, chemical screening, and genetic manipulations should facilitate deeper mechanistic understanding of GCGR functions and α-cell biology using this mutant line. Moreover, since it is hard to systemically dissect the genes network of GCGR in rodent models, we took advantage of zebrafish model, and performed RNA sequencing (RNA-seq) analysis of whole fish to provide a comprehensive view of its global transcriptomic changes in the absence of GCGR signaling.

The Transcriptome of Wildtype and Gcgr −/− Larvae
To survey the global transcriptomic changes of gcgr −/− mutants, we extracted total RNA from three biological replicates of both wildtype and gcgr −/− samples and performed high-throughput RNA-seq at seven days post fertilization (dpf). Together, these samples generated 21.94 million (M) pairs of raw reads, with a clean read ratio greater than 99.7%. After the quality filtering, clean reads were mapped to the zebrafish genome (GRCz11) using HISAT (v2.0.4) [26]. Among these reads, 89.55-90.55% clean reads were successfully mapped to the zebrafish genome, and the number of uniquely mapped reads was 63.98-65.53%. Then, clean reads were mapped to a reference transcriptome of zebrafish using Bowtie2 [27], and total 25,955 genes were detected. The average genes mapping ratio was 79.68% (79.28-80.03%), with genes uniquely mapped between 69.0%-69.21% (Table 1). Principal component analysis (PCA) demonstrated that the gcgr-mutant datasets clustered distinctly from wildtype control ( Figure 1A).

Gene Ontology (GO) Enrichment Analysis of Differentially Expressed Genes
A total of 24,831, 24,832, and 24,889 genes were detected with more than one read in WT1, WT2, and WT3, respectively, while 24,666, 24,629, 24,655 genes were detected in gcgr −/− 1, gcgr −/− 2, gcgr −/− 3, respectively. Comparison of RNA-seq data between the two genotypes identified by 1645 DEGs, with 437 upregulated and 1208 downregulated in gcgr −/− mutants ( Figure 1B and Supplemental Table S1). All the DEGs were annotated with gene ontology (GO) terms, and sorted into three major functional categories: "biological process", "cellular component," and "molecular function" (Figure 1C-E). Cellular process (22%), metabolic process (13%), and biological regulation (12%) were the main subcategories of the biological process group. Under the category of "cellular component," "cell" (25%) was the largest class, followed by "membrane" (16%), and "organelle" (16%). As for "molecular function", "binding" (47%) and "catalytic activity" (30%) were the major classes. These results of GO enrichment indicate that multiple biological processes changed in the gcgr-deficient fish. Since many functional categories were affected, we seek to determine whether there was any behavior or development difference between wildtype and gcgr −/− mutant. However, no significant difference in swimming speed and distance between mutant fish and wildtype fish were observed ( Figure 2). Moreover, there were no obvious defects in the development of gcgr −/− mutants compared to wildtype specimens (Supplemental Figure S1)

Kyoto Encyclopedia of Genes and Genomes (KEGG) Analysis Reveals that GCGR Deficiency Alters Multiple Metabolic Pathways
To better understand the GO-annotated DEGs in the zebrafish gcgr −/− mutant larvae, DEGs were subjected to the KEGG database for canonical signaling pathway analysis. The 1645 DEGs were significantly enriched in 44 different signaling pathways, and these pathways were most related to metabolism (12), human disease (11), organismal systems (10), cellular processes (4), genetic information processing (4), and environmental information processing (3).
Strikingly, we found that many genes involved in several metabolic processes were notably affected. Since GCGR is important for the regulation of metabolism, we then further analyzed these pathways. Interestingly, in metabolism pathways, the downregulated genes are far more than the upregulated genes. These data suggest that the GCGR knockout disrupts multiple metabolic pathways in zebrafish ( Figure 3).
To experimentally validate the dysregulation of lipid metabolism in gcgr −/− zebrafish, we performed oil-red O staining in whole mount and frozen sections in the wildtype and gcgr −/− fish in 7 dpf. The results showed a significant accentuation of fat in the liver ( Figure 4D-G). However, the fat accumulation in other tissues were not significantly different ( Figure 4D-G). These results suggest that GCGR is mainly required for lipid metabolism in the liver.
Four DEGs in glycolysis/gluconeogenesis pathway were enriched, two upregulation and two downregulation. The upregulated genes were adh5 (alcohol dehydrogenase 5) (2.01), which converts alcohol to aldehydes or ketones, and acss2l (acyl-CoA synthetase short chain family member 2 like) (1.37), which produces acetyl-CoA from acetate during gluconeogenesis. The two downregulated genes were adpgk2 (ADP-dependent glucokinase 2) (−1.41), which catalyzes ADP and D-glucose to AMP and D-glucose 6-phosphate during glycolysis, and aldh2.2 (aldehyde dehydrogenase 2 family member, tandem duplicate 2) (−1.11) which is responsible for the conversion of acetaldehyde to acetate and participate in glycolysis/gluconeogenesis pathways.
It is known that the blood free amino acid levels are increased in the GCGR knockout mice. This is concomitant with a decrease of intracellular amino acid content in hepatocytes [35]. However, it is unclear whether there is a decrease in the whole-body amino acid content. Here, we measured the total amino acid compositions of the whole organism, and found that Met, Ile, Tyr were slightly decreased in the gcgr −/− mutant zebrafish, while other amino acids did not significantly change ( Figure 6C).

Discussion
Glucagon is a counter-regulatory hormone in glucose homeostasis. It primarily acts on liver GCGR to increase blood glucose by promoting glycogenolysis and gluconeogenesis [1]. Therefore, most studies focus on the GCGR function in the liver. However, it should be appreciated that GCGR is also detectable in many other tissues, including adipose tissue, brain, kidney, heart, and gastrointestinal tract [25,36]. Metabolism is a complex process regulated by the crosstalk of multiple tissues, a systemic dissection of the metabolic network may help better understand the physiological function of GCGR.
In this study, we performed RNA-seq analysis using the whole organism of GCGR knockout zebrafish, and found that it regulates the lipid, carbohydrate, and amino acid metabolism networks.
The lipid metabolism was disrupted in the gcgr −/− knockout zebrafish. Fifty-two genes related to lipid metabolism were altered more than two-fold. Except for an increase of elovl7b in fatty acid elongation pathway, 11 genes enriched in fatty acid related metabolism were all decreased ( Figure 4A). That data suggest that fatty acid anabolism and catabolism are all decreased in the gcgr −/− mutant. For fatty acid oxidation, the hadhaa (hydroxyacyl-CoA dehydrogenase trifunctional multienzyme complex subunit alpha a), whose product catalyzes the last three steps of mitochondrial β-oxidation of long chain fatty acids, was significantly decreased. In addition, the acox1(peroxisomal acyl-coenzyme A oxidase 1), which also plays a crucial role in the fatty acid β-oxidation, was also downregulated ( Figure 4A). Taken together, these data suggest that the gcgr −/− mutants have decreased fatty acid β-oxidation compared to the wildtype fish.
We further found that the gcgr −/− mutants display a significant accumulation of lipids in the liver ( Figure 4D-G). Why lipids accumulate in the gcgr −/− liver is unknown. One possibility is that the gcgr −/− mutant may increase cholesterol metabolism based on our data. The dramatic increase of ldlr mRNA level in the gcgr −/− zebrafish predicts an elevated absorption of LDL-cholesterol (LDL-C) into hepatic cells. The decrease of pcsk9, which degrades LDLR, may further increase hepatic cholesterol. Consistent with our results, the GCGR antagonist-treated mice induces increased liver cholesterol absorption [17]. Additionally, T2D patients treated with glucagon receptor antagonist LY2409021 resulted in a statistically significant increase in hepatic fat fraction in a clinic trial [37]. Taken together, these data suggested that disruption of GCGR causes aberrant expression of lipid metabolism genes, which increased cholesterol absorption and resulted in the accumulation of lipid in the hepatic cells.
Glucose level is lower in GCGR knockout mouse, as well as zebrafish [6,25]. Additionally, the gcgr −/− mice are resistant to develop hyperglycemia in STZ-induced T1D and high-fat diet induced T2D models [7,8]. This phenotype has been attributed the associated increase of plasma GLP-1 levels [6,7,31,38]. Moreover, recent studies revealed that the improvement of glucose control by GCGR blockade required remnant insulin action in the diabetic animals [39,40]. Overall, previous studies have revealed that GCGR plays an important role in glucose metabolism. In our studies, the mRNA levels of most of the key enzymes in the glycolysis and gluconeogenesis pathways were unchanged, such as gck (−0.006), g6pca.1(0.38), g6pca.2 (0.31), g6pcb (0.30), gpia (0.07), pfkla (−0.17), pfklb (−0.09), aldoaa (0.06), aldoab (−0.01), and pgk1 (−0.05). However, the mRNA levels of several regulators of these enzymes were enriched more than two-fold in the glycolysis/gluconeogenesis pathway. The gckr (glucokinase regulator), which forms an inactive complex with glucokinase, were downregulated (−1.27), suggesting improved GCK activity in gcgr −/− animals. By contrast, the adpgk2 (ADP-dependent glucokinase 2) was downregulated (−1.41). Moreover, three other genes, aldh2.2, adh5, and acss2l were also affected. However, limited gene information prevents us from making a conclusive statement on the global glycolysis/gluconeogenesis change in the gcgr −/− mutants. Further information, such as global proteomic data can help us answer the question. In addition to glycolysis/gluconeogenesis pathway, 26 genes in other pathways of carbohydrate metabolism were also altered, suggesting that GCGR signaling may affect these genes, directly or indirectly, and they may contribute to the low blood glucose and the increased glucose tolerance.
We found that the gcgr −/− mutant zebrafish have an increased uptake of exogenous glucose using 2-NBDG as a surrogate. Similar results were found in mice. The GCGR antagonist REMD2.59 improves whole-body insulin sensitivity and glucose uptake in the ob/ob T2D model [41]. Moreover, when we cultured the gcgr −/− mutant larvae in a high glucose medium for three days, the α-cell hyperplasia was suppressed in the mutant ( Figure 5F). These data suggest that low glucose levels facilitate α-cell hyperplasia.
Several studies have confirmed that amino acid catabolism is decreased in mice due to decreased expression of transporters in the liver when GCGR signaling is disrupted, leading to increased plasma free amino acid levels [42][43][44][45]. Glutamine and alanine are dramatically increased and both contribute to α-cell hyperplasia [42,43]. In our study, the transcription level of several key amino acid catabolism genes was significantly decreased. For instance, L-amino acid oxidase gene il4i1, which catalyzes the oxidative deamination of a number of L-amino acids, was 26.4-fold decreased. Another gene, ido1 that regulates the rate-limiting step of oxidative tryptophan degradation, decreased 2.3-fold. Moreover, other amino acid metabolism genes were also dramatically affected. There was a 2.3-fold decrease in alas2 (Delta-aminolevulinate synthase 2), which catalyzes the condensation of glycine and succinyl coenzyme A to form 5-aminolevulinic acid (ALA) during heme biosynthesis. The gene gatm (guanidinoacetate N-methyltransferase), which catalyzes creatine synthesis from amino acids glycine, arginine, and methionine, was downregulated 2.4-fold. Moreover, nags (n-acetylglutamate synthase) that produces N-acetyl glutamate, which regulates the key enzyme of ureagenesis CPS1 (carbomyl phosphate synthetase-1) during ureagenesis, was decreased 4.1-fold. However, the total amino acid content only decreased slightly in gcgr −/− mutants. This reflects the importance of intracellular amino acid homeostasis. Nevertheless, all these data suggested that the zebrafish amino acid catabolism and related metabolism are downregulated. These changes most likely stem from the liver, similar to the mouse.
In summary, we performed whole-organism RNA-seq of wildtype and gcgr −/− mutant larvae. Compared with wildtype, the mutant data showed that many genes in metabolic processes were changed more than two-fold. By further analysis, we found that most of these genes are related to lipid metabolism, carbohydrate metabolism, and amino acid metabolism. Of note, genes related to fatty acid β-oxidation, amino acid catabolism, and ureagenesis were downregulated. Moreover, we experimentally confirmed that the mutants display an increased accumulation of lipids in the liver, an increase in glucose uptake, and modest changes in total amino acids. These results provide new information about the global metabolic network that GCGR signaling regulates in addition to a better understanding of the receptor's physiological functions. Nevertheless, we still do not know that these changes in the gcgr −/− mutants were due to the direct or indirect effects of GCGR deficiency. Further studies need to perform transcriptomic analyses of the gcgr −/− mutants at different stages, as well as in different tissues.

Zebrafish Lines and Maintenance
Zebrafish (Danio rerio) were raised in a recirculating aquaculture system (Shanghai Haisheng Biotech Co., Ltd., Shanghai, China) on a 14:10 h darkness cycle at 28 • C. Embryos were obtained from natural breeding and raised at 28.5 • C in an embryo rearing solution, and staged according to Kimmel et al. (1995) [46]. In this study, the AB strain, gcgra −/− ; gcgrb −/− double mutant fish (referred as gcgr −/− henceforth) [25], and Tg(gcga:GFP) [47] were used. All procedures have been approved by the Xiamen University Institutional Animal Care and Use Committee (Protocol XMULAC20160089, 10 March 2016).

RNA Extraction, cDNA Library Preparation, Sequencing
Total RNA was isolated from pools of 30 larvae using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer's instruction. Three biological repeats were performed for each genotype. To remove any genomic DNA contamination, the RNA were then digested using the RQ1 RNase-Free DNase (Promega, Madison, WI, USA). The concentration and quality of each sample were determined by an Agilent RNA 6000 nano kit in the Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The RNA quality score (RIN/RQN) was 10 for all of the samples.
To construct cDNA libraries, the mRNA of each sample was selected using Oligo (dT) magnetic beads. The total mRNA was then fragmented and reverse-transcribed into double-strand cDNA (ds cDNA) by N6 random primers. The synthesized double-strand cDNA was subjected to end repair, follow by 5 phosphorylation and 3 adenylation. Next, an adapter was ligated to the 3 adenylated cDNA fragments. The ligation products were used for PCR amplification for enrichment of cDNA templates. Following denaturation of these PCR products, the single stranded DNA was then circlized to form the final cDNA library. All the sequencing of these cDNA libraries was performed on a BGISEQ-500 RS platform (BGI, Shenzhen, China) according to the manufacturer's stand protocol [48].

Bioinformatics Analysis of RNA Sequence Data
Raw reads were subjected to a BGISEQ-500 quality control test and filtered into a clean read using SOAPnuke software [49]. This process discarded reads containing adaptor sequences, reads containing unknown bases "N" more than 10%, and low quality reads (the percentage of low quality bases more than 50% in a read). After filtering, the clean reads were aligned with the zebrafish genome (Danio rerio, GRCz11, http://asia.ensembl.org/Danio_rerio/Info/Index) using HISAT (v2.0.4) [26]. Following alignment, the clean reads to zebrafish unigenes using Bowtie2 [27], RSEM (RNA-Seq by expectation maximization) [50] was performed to quantify the gene expression using the FPKM method (fragments per kilobase of transcript per million fragments sequenced) [51]. Moreover, some uncharacterized genes were annotated in the Ensembl (www.ensembl.org/) by BLAST and synteny analysis. These sequence data were deposited in an NCBI database and are accessible via BioProject ID: PRJNA541367.

Analysis of Differentially Expressed Genes (DEGs)
The DEseq2 was used to identify DEGs [52]. Genes with fold change ≥ 2.00 and adjusted p-value ≤0.05 were considered as differentially expressed genes (DEGs) with statistical significance. For duplicate samples, the log2 fold change (Log2FC) and probability for each gene in every comparison was calculated by the NOIseq under the conditions with fold change ≥ 2.00 and probability ≥ 0.8 [53].
To perform the pathway enrichment analysis, the significant DEGs were subjected to GO (gene ontology) and KEGG (Kyoto encyclopedia of genes and genomes) database analysis. The analysis was performed using the hypergeometric test and false discovery rate (FDR) correction methods. Significantly enriched genes were those with a p-values < 0.05 and a FDR ≤0.01.

Quantitative RT-PCR
To validate the DEGs, qRT-PCR was performed. To synthesize the first strand cDNA, the RNA was reverse transcribed by the M-MLV (Promega, Madison, WI, USA) with oligo(dT)16 primers. The qRT-PCR reaction was carried out in an Agilent AriaMx system (Agilent Technologies, Santa Clara, CA, USA) using powerUp SYBR Green Master Mix (Thermo Fisher Scientific, Carlsbad, CA, USA). The amplification program was 95 • C, 10 s; 60 • C, 30 s for 40 cycles. Three biological samples of wildtype or gcgr −/− were used, and each assay for a sample was performed in triplicate. mRNA levels were calculated using the 2 −∆∆Ct method [54] and presented as relative (fold) levels normalized to the level of β-actin. The data were presented as the mean of three biological samples. The primers used were listed in Supplemental Table S6.

Oil Red Staining
Zebrafish at 7 dpf were fixed in 4% paraformaldehyde in PBS for 1 h at room temperature. The fish were washed with PBS subsequently to remove paraformaldehyde, incubated in 60% 2-propanol for 10 min, followed by freshly filtered 0.3% oil red O (ORO) (Xiya Reagent, Chengdu, China) in 60% 2-propanol for 10 min. The staining was terminated by washing twice with 60% 2-propanol. The samples were transferred to 75% glycerol and imaged using Lecia M205 FCA (Lecia Wetzlar, Germany).

2-NBDG Uptake Test
Zebrafish (5 dpf) were incubated in a culture medium containing 600 uM 2-NBDG (Apexbio, B6035, Houston, TX, USA) for 3 h. The larvae were anesthetized for imaging under a M205 FCA microscope (Leica, Wetzlar, Germany), the fluorescence intensity of lens was used as the indicator of glucose uptake according to the reference [32].

Glucose Exposure and α-Cells Counting
The Tg(gcga:GFP) and gcgr −/− ; Tg(gcga:GFP) larvae were incubated with 20 mM glucose or a 0.3× Danieau solution from 4 to 7 dpf for three days. After harvest, larvae were fixed with 4% paraformaldehyde overnight at 4 • C, and then placed on a slide with aqua-mount (Richard-Allan Scientific, Kalamazoo, MI, USA) with the right side of larvae up to expose the islet. The α-cell number was counted according to the GFP under a Zeiss AxioImager A1 microscope (Carl Zeiss, Jena, Germany) with 40× lens.

Behavior Test
Seven dpf zebrafish were transferred individually to a compartment containing 1 mL 0.3X Danieau solution in a 24-well plate. Then, the plate was video-recorded by a top-view camera in the DanioVision tracking system (Noldus IT, Wageningen, Netherlands) for 3 min, the moving tracks were generated and the total moving distance and velocity were determined by Ethovision XT7 software (Noldus IT, Wageningen, Netherlands), 48 animals from each group were analyzed.

Amino Acid Measurement
Seventy embryos of 7 dpf for each group were harvested in 1.5 mL Eppendorf tubes. After lyophilization, samples were hydrolysed in 6 N HCl at 110 • C for 24 h. Then, the volume was adjusted to 10 mL by adding distilled water. One ml of hydrolysate was dried using nitrogen blower and resuspended in 1 mL 0.02 M HCl. After resuspension, samples were filtered through 0.22 µm filters, the amino acid compositions were determined by an automatic amino acid analyzer (L-8900, HITACHI, Chiyoda, Japan). For WT and gcgr −/− , three independent biological repeats were performed. Amino acid in the embryos were expressed as µg/mg wet tissue.

Statistical Analysis
Data were presented as means or means ± S.E.M. The statistics were performed using one-way ANOVA followed by Newman-Keuls post-hoc test or t-test (SPSS). p < 0.05 was considered as significant.