Involvement of IGFBP5 in the Development of Goose Fatty Liver via p38 Mitogen-Activated Protein Kinase (MAPK)

: Previous studies showed that insulin-like growth factor-binding protein 5 (IGFBP5) plays a role in non-alcoholic fatty liver disease; however, its expression and function in goose fatty liver remain unknown. To address this, we obtained a full-length mRNA sequence of the goose IGFBP5 gene using a 5 (cid:48) -rapid ampliﬁcation of cDNA ends assay and nested polymerase chain reaction (PCR). Additionally, using the newly acquired sequence of 5’-untraslated region, we determined the missing sequence of the ﬁrst intron. Bioinformatics analysis revealed three exons and three introns in the goose IGFBP5 gene. Quantitative PCR analysis indicated that the mRNA abundance of IGFBP5 was signiﬁcantly lower in goose fatty liver than in the normal liver. Comparison of transcriptomes of goose primary hepatocytes transfected with IGFBP5 overexpression vector versus those transfected with empty vector identiﬁed 777 differentially expressed genes (DEGs). The enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways indicated the focal adhesion, ECM-receptor interaction, regulation of the actin cytoskeleton, mitogen-activated protein kinase (MAPK) signaling, and GnRH signaling pathways. Immunoblotting revealed the induction of the p38 MAPK pathway by IGFBP5 overexpression, which is in line with the suppressed expression of IGFBP5 and p38 MAPK in goose fatty liver than in normal liver. These ﬁndings suggest that IGFBP5 is involved in the development of goose fatty liver via the p38 MAPK pathway. to validate the differentially expressed genes identiﬁed by transcriptome analysis.


Introduction
The liver plays an important role in insulin-mediated regulation of metabolism, especially through glucose and lipid homeostasis [1,2]. It is also a major site for the synthesis of endocrine factors, such as insulin-like growth factors 1 and 2 (IGF1 and IGF2) and their binding proteins (IGF-binding proteins, IGFBPs) [3]. Several metabolic disorders, including non-alcoholic fatty liver disease (NAFLD), are associated with liver dysfunction. NAFLD is characterized by hepatic steatosis and is usually accompanied by insulin resistance (IR) [4][5][6]. NAFLD comprises multiple pathological syndromes, ranging from simple steatosis and steatohepatitis to fibrosis and even cirrhosis [6]. In contrast, fatty liver in geese is physiological rather than pathological and is considered a consequence of evolutionary adaptation to the needs for migration [7]. Goose fatty liver does not show any overt pathological syndromes [2,6]. Previous studies identified several protective components, including downregulation of TNF-α and complement genes, as well as upregulation of adiponectin receptors, fatty acid desaturases, and mitochondrial genes [6,8]. Goose liver has a great capacity for lipid deposition; 3-4 weeks of artificial overfeeding can result in goose liver weighing 8 to 10 times its normal weight [6]. Thus, the goose fatty liver is a unique model for NAFLD studies. However, the mechanisms underlying its formation remain unclear.
IGFBPs are a family of proteins that bind to IGFs with high affinity, including six (IGFBP1-6) or more IGFBPs [3,9,10]. IGFBP1-6 have masses of 24-50 kDa (216-289 amino acids) and share a highly conserved structure with three domains of similar size: the conserved N-terminal cysteine-rich and the C-terminal cysteine-rich regions are connected by a less structural and less conserved linker region [9]. The biological functions of IGFBPs can be IGF-dependent or IGF-independent. IGF protein hormones play extensive roles in cell cycle progression, growth, proliferation, metabolism, and overall cell survival [1,3,[9][10][11].
IGFBP5 is one of the most evolutionarily conserved members of the IGFBP family [12]. The human IGFBP5 protein contains 252 amino acids with a globular high-affinity IGF binding region at the N-terminus of IGFBP5 [13]. IGFBP5 has a modest binding preference for IGF2 over IGF1 and can inhibit or potentiate IGF activity in a range of cells in vitro. However, IGFBP5 also has several IGF-independent actions in vitro. IGFBP5 exerts its transactivation activity via its N-terminal domain within the nucleus [13]. IGFBP5 plays an important role in the bone, ovary, mammary gland, kidney, and liver [12]. The main function of IGFBP5 is to regulate cell growth by either suppressing or inducing cell proliferation [12]. IGFBP5 has been implicated in some diseases [13]. For example, a study showed that the level of IGFBP5 in the serum is correlated with NAFLD [14]; thus, serum IGFBP5 may serve as a biomarker for NAFLD [15]. However, the expression and function of IGFBP5 in goose fatty liver remains unknown.
In this study, we performed sequence analysis of goose IGFBP5 gene and compared the IGFBP5 expression levels in goose fatty liver with those in normal liver by quantitative polymerase chain reaction (PCR). We investigated the potential role of IGFBP5 in the formation of goose fatty liver by transcriptome analysis of goose primary hepatocytes transfected with IGFBP5 overexpression vector vs. empty vector and validated a key pathway mediating the role of IGFBP5 in goose fatty liver vs. that in normal liver. This study lays a foundation for addressing the role of IGFBP5 and its related mechanisms in the development of goose fatty liver.

Experimental Animals and Sample Collection
A total of 66-d-old healthy male Landes geese from the same batch raised in Licheng Livestock and Poultry Breeding Co., Ltd. (Huai'an, Jiangsu, China) were used in this study. The geese were randomly divided into two groups: control and over-feeding (n = 8). Diet and dietary schemes were applied according to our previous study [16]. Briefly, all geese were housed in cages and raised under natural light and temperature conditions. For the overfeeding group, there was a 5 day pre-overfeeding period, followed by a 19 day formal overfeeding period, during which the amount of diet and the frequency of feeding were gradually increased to 1200 g of daily intake, five meals per day. Geese in the control group were fed ad libitum. All geese had free access to water and the feed used for each group was the same. On the last day of overfeeding, six geese were randomly selected for overnight fasting followed by sacrifice for sample collection. The liver samples were collected and snap-frozen in liquid nitrogen before being stored at -70 • C. All animal protocols were approved by the Animal Welfare Committee of Yangzhou University (permission number SYXK (Su) IACUC 2016-0020).

Primer Design
The primers used in this study were designed using PRIMER 3 software (http://www. bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi, accessed on 21 December 2020) or the NCBI Primer-Basic Local Alignment Search Tool (BLAST) software (https://www.ncbi. nlm.nih.gov/tools/primer-blast/, accessed on 13 May 2021). The sequences of the primers used to acquire and validate the mRNA and genomic sequences of the goose IGFBP5 gene are listed in Table 1, and the sequences of the primers used for qPCR are listed in Table 2. Genomic and mRNA sequences of the goose IGFBP5 gene (NW_013185861.1, XM_013196711.1, and a new sequence from this study) were used for primer design. Note: 1 The primers were designed to acquire and validate the 5 -end of goose mRNA sequence and the missing intronic sequence between the 5 -UTR and the first exon of goose IGFBP5 gene. 2 The underlined sequence GATTACGCCAAGCTT added to the 5 ends of IGFBP5-GSP1/2/3/4 primer sequences were used to facilitate in-fusion cloning of PCR products. Table 2. Information on the primer sets used in quantitative PCR analysis.

Isolation of Genomic DNA (gDNA)
Total gDNA was extracted from a normal liver sample using a TIANamp Genomic DNA Kit (DP304, TIANGEN Biotech, Co, Ltd., Beijing, China) according to the manufacturer's instructions. The quality and quantity of gDNA samples were determined using a NanoDrop spectrophotometer (NanoDrop Technologies, Inc. Wilmington, DE, USA). The gDNA samples were stored at −20 • C until use.

General PCR and Product Separation by Electrophoresis
Each PCR mixture contained 1 µL (100 ng) of DNA template, 1 µL (10 µM) of for-ward primer, 1 µL (10 µM) of reverse primer, 10 µL of 2× Taq Plus PCR MasterMix (Vazyme Biotech, Co., Ltd., Nanjing, China), and 7 µL water. PCR was performed on a Veriti™ 96-Well Thermal Cycler (Applied Biosystems Instruments, Singapore) under the following conditions: initial denaturation at 95 • C for 3 min, 35 cycles of denaturation at 95 • C for 15 s, annealing at 60 • C for 15 s, extension at 72 • C for 1 min, and a final extension at 72 • C for 5 min. Electrophoresis of the PCR product was carried out on a 1% agarose gel stained with ethidium bromide. Gel images were obtained using a UV transillumintor (Heqin ZF1-II, Shanghai, China).

Total RNA Isolation and 5 Rapid Amplification of cDNA Ends (RACE)
Total RNA was isolated from the liver tissues of Landes geese using RNAzol (DP-424 RNAiso Plus, Takara Co. Ltd, Dalian, China) according to the manufacturer's instructions. Using total RNA samples, we performed 5'-RACE using a SMARTer RACE Kit (Cat. No. 634860, Takara Bio USA, Inc., Mountain View, CA, USA), according to the manufacturer's instructions. Briefly, 5.5 µL of Buffer Mix I was prepared by mixing 4.0 µL of 5× First-Strand Buffer, 0.5 µL of 100 mM dithiothreitol, and 1 µL of 20 mM dNTPs. Buffer Mix I was placed on ice until use. Simultaneously, a 12 µL volume of Buffer Mix II was prepared by mixing 5 µL of 1 µg/µL RNA, 1 µL of 12 µM 5'-CDS Primer A, and 5 µL H 2 O, followed by incubating the mixture at 72 • C for 3 min, cooling at 42 • C for 2 min on a Veriti™ 96-Well Thermal Cycler (Applied Biosystems instruments, Singapore), spinning at 14,000 × g for 10 s, and adding 1 µL of 24 µM SMARTer II A Oligonucleotide. Subsequently, 8 µL Master Mix was prepared by mixing 5.5 µL of Buffer Mix I, 0.5 µL of 40 U/µL RNase inhibitor, and 2 µl of (100 U) SMART-Scribe reverse transcriptase. For the first strand cDNA synthesis, a 20 µL reaction solution was prepared by mixing 8 µL Master Mix and 12 µL Buffer Mix II, which was followed by incubation at 42 • C for 90 min, heating at 70 • C for 10 min, and storage at −20 • C until use.
After the first-strand cDNA synthesis, cDNA was amplified with the Universal Primer Mix (UPM) and a gene-specific primer (named as GSP1) on a Veriti™ 96-Well Thermal Cycler (Applied Biosystems instruments, Singapore). The reaction volume was 50 µL, which contained 15. After 5'-RACE, nested PCR was performed with a total volume of 50 µL reaction solution, which contained 17 µL H 2 O, 25 µL 2 × SeqAmp™ buffer, 1 µL SeqAmp DNA polymerase, 1 µL of the Universal Primer Short (UPS), 1 µL of the nested primers (named GSP2/3/4, separately), and 5 µL of the RACE product. The PCR conditions were 25 cycles at 94 • C for 30 s, 68 • C for 30 s, and 72 • C for 3 min.

Cloning of PCR Product
For cloning, the PCR products were first separated by electrophoresis on a 1% agarose gel stained with ethidium bromide, followed by excision and purification of the separated products using the NucleoSpin Gel and PCR Clean-Up Kit (Cat. No. 740986.20, Takara Bio USA, Inc., Mountain View, CA, USA) according to the manufacturer's instructions. The purified products were cloned into the linearized pRACE vector (PUC-19) using the In-Fusion HD Cloning kit (Cat. No. 634860, Takara Bio USA). The vectors with the inserts were transformed into Stellar™ competent cells, which were immediately placed on plates filled with Luria-Bertani (LB) culture medium containing 100 µg/mL ampicillin. The transformed cells were subsequently incubated at 37 • C overnight. Twelve colonies were selected from the plates and individually incubated overnight in LB broth in tubes at 37 • C. After harvesting the cultured cells from the colonies, we extracted DNA from the cells using the TIANprep Rapid Mini Plasmid Kit (DP105-03, TIANGEN Biotech, Co, Ltd., Beijing, China) according to the manufacturer's instructions and sent it to General Biosystems (Anhui) Co., LTD. (Chuzhou, China) for sequencing.

Bioinformatics Analysis of Goose IGFBP5 mRNA Sequence
After acquiring the full-length mRNA sequence of goose IGFBP5 gene, we applied the "ORF finder" program of National Center for Biotechnology Information (NCBI) (http://www. ncbi.nlm.nih.gov/gorf/orfig.cgi, accessed on 15 February 2021) to predict the open reading frame (ORF). The conserved domain architecture retrieval tool was obtained from the National Center for Biotechnology Information server (http://www.ncbi.nlm.nih.gov/ Structure/cdd/docs/cdd_search.html, accessed on 21 February 2021). Sequence alignment analysis was performed using BLAST software (https://blast.ncbi.nlm.nih.gov/Blast.cgi, accessed on 21 February 2021).

Acquiring the Missing Part of the First Intron Sequence of Goose IGFBP5 Gene
After the 5 -UTR of the goose IGFBP5 gene was obtained from 5 -RACE analysis, the primers were designed on the basis of the 5'-UTR sequence and the known sequence of the first intron of the goose IGFBP5 gene in GenBank. Using the primers and genomic DNA isolated from goose liver samples, we performed general PCR as described above. PCR products were sent to General Biosystems (Anhui) Co., Ltd. (Chuzhou, China) for sequencing. The complete genomic sequence of the goose IGFBP5 gene was assembled using the sequence acquired (submitted to the GenBank database and assigned with accession No. MW362807) and the known genomic sequence (NW_013185861.1) in GenBank.

Quantitative PCR Analysis
Total RNA was isolated from the tissue and cell samples using RNAzol (DP-424, RNAiso Plus, Takara Co. Ltd, Dalian, China), and the concentration and quality of the isolated RNA were determined using a NanoDrop spectrophotometer (NanoDrop Technologies, Inc. Wilmington, DE, USA). Using the RNA samples, we synthesized cDNA using HiScript Q RT SuperMix (R123-01, Vazyme Biotech Co., Ltd., Nanjing, China) according to the manufacturer's instructions. Subsequently, the relative mRNA abundance of the target genes was determined by quantitative PCR, as previously described [18]. Briefly, 10 µL of AceQ qPCR SYBR Green, 0.4 µL 10 µM forward primer, 0.4 µL 10 µM reverse primer, 0.4 µL of 50× ROX reference Dye 1, 2 µL (1 ng) cDNA template, and 6.8 µL water were mixed, and PCR was performed using a QuantStudioTM 3 Real-Time PCR instrument (A28132, Applied Biosystems instruments, Singapore) under the following conditions: 95 • C for 5 min, 40 cycles of 95 • C for 10 s, and 60 • C annealing for 30 s. The glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene was used as an internal control for normalization. The cycle threshold (Ct) was determined the using supplied software. Relative mRNA abundance was calculated using the 2 −∆∆Ct method [19] and presented as fold-change over the control.

Isolation and Treatment of Goose Primary Hepatocytes
Goose primary hepatocytes were isolated from Landes goose embryos on the 23rd day of hatching, as previously described [8]. Isolated primary hepatocytes were cultured as previously described [20]. In brief, the isolated cells were seeded at a density of 1 × 10 6 cells per well and cultured in complete culture medium (containing high glucose DMEM culture medium, 10% fetal bovine serum, 1% penicillin-streptomycin solution (100 IU/mL), and 10 µL EGF (20 ng/mL)) overnight. For transfection assays, the cells at 80% confluence were transfected using Lipofectamine™ 2000 (Cat. No. 11668-027, Invitrogen™, San Diego, CA, USA) according to the manufacturer's instructions. The overexpression vector was constructed with the full coding sequence of the goose IGFBP5 gene and pcDNA3.1 (+) vector (Cat. No. V79020, Thermo Fisher Scientific Co., Shanghai, China).

Transcriptome Analysis
Total RNA was isolated from goose primary hepatocytes transfected with goose IGFBP5 gene overexpression vector and an empty vector was used for cDNA synthesis by reverse transcription, followed by establishment of a cDNA library and sequencing. The original reads obtained by sequencing were filtered using the Fast QC software, and the reference goose genome of was used for annotation, followed by calculating the significance of the difference in the expression of each gene between the cells transfected with the empty vector and the overexpression vector using DESeq software (http://bioconductor.org/packages/ release/bioc/html/DESeq.html). Genes with|log2 (fold change)|>1 and p-value < 0.05 were the criteria for DEG selection. With the identified DEGs, GO and KEGG enrichment analyses were performed using the ClusterProfiler software package. The standard for significant enrichment was set at a p-value < 0.05.

Statistical Analysis
All values are expressed as mean ± standard error. Statistical significance was determined using Student's t-test for two-group comparisons using SPSS 16.0 (SPSS China, Shanghai, China). Statistical significance was set at p < 0.05.

Acquiring the Full Length of Goose IGFBP5 mRNA
To investigate the function of the goose IGFBP5 gene, we required a full-length mRNA sequence. However, the mRNA sequence of goose IGFBP5 in GenBank was incomplete because it lacked the 5 -untranslated region (5 -UTR). To address this, we synthesized complementary DNA (cDNA) and amplified it by 5 -rapid amplification of cDNA ends (5 -RACE) with UPM and GSP1 (Table 1, Figure 1A) using total RNA samples isolated from the liver of Landes geese. As the RACE product showed a high level of background with no specific bands in the electrophoresis analysis ( Figure 1B), nested PCR was performed using UPS and GSP2/3/4 primers (Table 1, Figure 1A) with the RACE product as the template, and a strong band of approximately 400 bp was detected by electrophoresis ( Figure 1B). This band was excised, purified, and ligated with a linearized pRACE vector (PUC19), followed by the transformation of Stellar™ competent cells with the ligation product. A total of 12 bacterial colonies were isolated and sequenced, and a 420 bp sequence fragment was identified, which was bordered with UPS and GSP3 primer sequences. Using the online BLAST program, this sequence fragment was aligned against the incomplete mRNA (XM_013196711.1) sequence of the goose IGFBP5. The result indicated that the fragment sequence was identical to the sequence of XM_013196711.1 in the overlapping region of 86 nt. On the basis of this alignment analysis, we assembled the full-length goose IGFBP5 mRNA sequence (990 nt), including the 5'-UTR (375 nt), coding sequence (438 nt), and 3 -UTR (177 nt) ( Figure 2). The newly assembled sequence was assigned the GenBank accession number MW351707. The possibility of contamination of goose genomic DNA in the 5'-RACE analysis was excluded by performing general PCR with the template cDNA and primers for the genomic sequence of the Actinin alpha 1 (ACTN1) gene ( Figure S1).

Bioinformatics Analysis of Goose IGFBP5 mRNA Sequence
Multiple sequence alignment of goose IGFBP5 mRNA with its counterparts from other animals (including duck, chicken, quail, human, mouse, and cattle) revealed a high similarity (77.23-96.81%) among all the orthologs, especially for the poultry species (more than 91%) ( Table 3).
On the basis of the coding sequence of the goose IGFBP5 mRNA sequence, we found the amino acid sequence of the goose IGFBP5 gene to contain 145 amino acid residues, including the start and stop codons (Figure 2). There were two conserved domains, namely, the IGFBP homolog domain (from 70 to 300 nt) and the thyroglobulin type-1 domain (from 568 nt to 783 nt) in the predicted IGFBP5 protein ( Figure 3A). The former domain may be involved in the control of proteolytic degradation of IGFBP5 protein, while the latter may be responsible for high-affinity binding to IGF. Moreover, multiple sequence alignment of goose IGFBP5 protein with its counterparts from other animals revealed 78.07-99.62% identity of amino acid sequence among the orthologs, especially for poultry species (more than 99%) (Table 3).
Notably, when IGFBP5 sequence of the goose was aligned against that of the duck, the 5 -UTR of the goose IGFBP5 had a high identity to the first exon of the duck IGFBP5 sequence, except for a short sequence of 29 bp at the beginning of the 5 -UTR of the goose IGFBP5, which did not match the duck IGFBP5 sequence ( Figure 3B). This finding suggests that the original first exon of the goose IGFBP5 may have lost its start codon in the open reading frame during evolution and became the 5 -UTR.

Acquiring the Full Length of the Goose IGFBP5 Genomic Sequence
The genomic sequence of the goose IGFBP5 gene (NW_013185861.1) in GenBank lacked an intronic sequence between the 5 -UTR and exon 1. To acquire the missing sequence, we aligned the goose IGFBP5 mRNA sequence with the genomic sequence of the duck IGFBP5 (NC_051778.1: 6664368. . . 6682347). On the basis of the sequence alignment, we designed primers (IGFBP5-P1) and used them to amplify the missing sequences by PCR using goose genomic DNA as a template (Table 1, Figure 4A). The PCR products ( Figure 4A) were sequenced, followed by sequence assembly using the known sequence of the goose IGFBP5 gene. The newly assembled sequence of goose IGFBP5 gene was aligned against the duck IGFBP5 gene, which indicated that there was 95.05% identity between the goose and duck IGFBP5 genomic sequences with 97% query coverage ( Figure 4B). The missing intronic sequence of the goose IGFBP5 was subsequently validated by PCR with primers including IGFBP5-P1, IGFBP5-P2, and their combination, using genomic DNA templates isolated from goose liver samples (Table 1). As expected, electrophoresis of PCR products obtained with IGFBP5-P1 and IGFBP5-P2 primers produced specific bands of approximately 1400 and 290 bp, respectively ( Figure 4A). In contrast, use of the IGFBP5-P1 forward primer and IGFBP5-P2 reverse primer did not result in any amplicons with goose genomic DNA due to the large distance (more than 17 kb) between the 5'-UTR and exon 1 sequence, or approximately 320 bp of amplicons with cDNA samples ( Figure 4A). These results suggest that the missing sequence determined here belongs to the goose IGFBP5 gene, and that the whole sequence of goose IGFBP5 gene is 18481 bp long, containing three exons (exon 1: 187 bp, exon 2: 119 bp, exon 3: 132 bp) and three introns (intron 1: 17012 bp, intron 2: 145 bp, intron 3: 363 bp). This sequence has been submitted to GenBank and a provisional accession number has been assigned (No. MW362807).

Expression of IGFBP5 Was Suppressed in Goose Fatty Liver vs. Normal Liver
The IGFBP5 expression levels in goose fatty liver vs. normal liver were determined by quantitative PCR analysis using primers designed on the basis of the goose IGFBP5 coding DNA sequence retrieved from GenBank ( Table 2). The data indicated that the mRNA levels of IGFBP5 were significantly lower in the goose fatty liver than in the normal liver ( Figure 5), suggesting that IGFBP5 may play an important role in the development of goose fatty liver.

Identification of Genes and Pathways Affected by IGFBP5 Overexpression in Goose Primary Hepatocytes
To determine the potential role of IGFBP5 in the formation of goose fatty liver, we transfected goose primary hepatocytes with an IGFBP5 overexpression vector or empty vector as a control. The results showed that the mRNA expression levels of IGFBP5 in goose primary hepatocytes transfected with the IGFBP5 overexpression vector were significantly higher than in those transfected with the empty vector ( Figure S2), indicating that IGFBP5 was overexpressed in goose hepatocytes. Cells transfected with the IGFBP5 overexpression vector (n = 6) or empty vector (n = 6) were used for transcriptome analysis. In the duck IGFBP5 gene, there is no intronic sequence between the 5 -untranslated region (5 -UTR) and the first exon. In the goose IGFBP5 gene, the red boxes denote the sequence acquired from the 5'-RACE analysis. Although the 5 -UTR sequence in the goose IGFBP5 gene is almost identical to the first exon of the duck IGFBP5 gene, it does not contain the start codon of open reading frame. Therefore, the sequence matching to the first exon of duck IGFBP5 becomes the 5 -UTR sequence of goose IGFBP5. The red "NNN" letters denote the intronic sequence missing in GenBank reference sequence. The "29 bp" denotes that the sequence of 29 bp long does not align with any part of duck IGFBP5 sequence.
A total of 85220702 raw reads and 83519842 clean reads were generated for cells transfected with the empty vector, and 91085158 raw reads and 88900062 clean reads were generated for cells transfected with the overexpression vector ( Table 4). The values of QC30 ranged from 92.59 to 94.42%, indicating that the sequencing quality met the criteria for further analysis. To identify the DEGs affected by IGFBP5 overexpression, we calculated the relative abundance of a given gene and presented it as reads per kilobase per million reads (RPKM) with|log2 (Fold Change)|>1 and p-value < 0.05 as the criteria for DEG selection. A total of 777 DEGs were subsequently identified, including 80 downregulated and 697 upregulated, in cells transfected with IGFBP5 overexpression vector vs. empty vector. The top 10 downregulated and upregulated DEGs are listed in Tables 5 and 6. KEGG pathway analysis showed that the DEGs were mostly enriched in the focal adhesion, ECMreceptor interaction, regulation of actin cytoskeleton, MAPK signaling, and GnRH signaling pathways ( Figure 6).  Note: 1 IGFBP5-S1 and -S2 are the samples of goose primary hepatocytes transfected with goose IGFBP5 overexpression vector, while Control-S1 and -S2 are the samples of goose primary hepatocytes transfected with empty vector. 2 Q20 and 3 Q30 denote the percentages of the base pair with a Phred value greater than 20 and 30 in total base pair, respectively.   Some DEGs were randomly selected for validation by quantitative PCR using the primers designed on the basis of their mRNA sequences retrieved from GenBank (Table 2), including diazepam-binding inhibitor (DBI), cytochrome c oxidase subunit III (COX3), ribosomal protein L37 (RPL37), interferon-induced protein with tetratricopeptide repeats 5-like (LOC106044214), interferon-induced GTP-binding protein Mx transcript variant X1 (LOC106036132), phosphatase and tensin homolog (PTEN), kinase insert domain receptor (KDR), tenascin C transcript variant X5 (TNC), Wiskott-Aldrich syndrome-like transcript variant X2 (WASL), 2 -5 -oligoadenylate synthase-like protein 2 (LOC106042804), p21 protein (Cdc42/Rac)-activated kinase 2 (PAK2), and actinin alpha 1 transcript variant X5 (ACTN1). The data indicated that all of these DEGs except for ACTN1 and PAK2, that is, 10 out of 12 genes, were validated (Figure 7), suggesting that the results from the transcriptome analysis were generally reliable.

Association of Total p38 MAPK Protein Levels with IGFBP5 Expression
MAPK signaling pathways include p38, ERK, and JNK MAPK cascades. In this study, the MAPK pathway was significantly enriched in DEGs. As the antibodies used in this study did not work well with goose ERK and JNK MAPK proteins (no expected blots based on their molecular weights), only the protein expression of p38 MAPK was determined in goose primary hepatocytes transfected with IGFBP5 overexpression vector vs. empty vector by immunoblotting. Total p38 protein levels were higher in goose primary hepatocytes transfected with IGFBP5 overexpression vector than in those transfected with empty vector, whereas they were lower in goose fatty liver than in normal liver ( Figure 8). These findings suggest that total p38 protein level is associated with IGFBP5 expression in goose liver cells.

Discussion
Complete mRNA sequences are important for functional gene studies. On the basis of the mRNA sequence, one is able to predict the start point of transcription, the regions of upstream regulatory elements, the amino acid sequence encoded by the gene, chemical modifications, and other features of the protein. In this study, the missing part at the 5 -end of goose IGFBP5 mRNA was determined using 5 -RACE analysis. Full-length goose IGFBP5 mRNA was assembled, and the amino acid sequence of goose IGFBP5 protein was predicted. The protein contains the IGFBP homolog domain and thyroglobulin type-1 domain based on alignment analysis with human IGFBP5 protein. The thyroglobulin type-1 domain may be responsible for the degradation of this protein [21], whereas the IGFBP homolog domain may be responsible for its binding to IGF [9]. Moreover, the missing part of the first intron in the goose IGFBP5 gene was determined using the 5 -UTR sequence acquired from 5'-RACE. Taken together, the complete mRNA and genomic sequences of the goose IGFBP5 gene provide a solid foundation for functional gene studies.
In this study, the comparisons of the mRNA and protein sequences among different animals, namely, goose, duck, chicken, quail, human, mouse, and cattle, indicate that IGFBP5 is relative conserved among species, especially among avian species. This is consistent with previous reports showing that IGFBP5 is evolutionally the most conserved IGFBP family member [12,22]. Interestingly, comparison of the goose and duck IGFBP5 sequences revealed that the goose IGFBP5 gene has three exons and three introns, which is different from the duck IGFBP5 gene with four exons and three introns. This difference is probably due to the mutative loss of the start codon originally existing in the first exon of the goose IGFBP5 gene during evolution. Moreover, there was no sequence alignment between the duck IGFBP5 sequence and the beginning of the 5'-UTR of the goose IGFBP5 gene (29 bp), which also suggests that mutations occurred at the beginning of the 5'-UTR of the goose IGFBP5 gene during evolution. However, it is unknown as to whether these mutations contribute to functional changes in goose IGFBP5.
This study revealed that IGFBP5 mRNA expression was inhibited in goose fatty liver compared to normal liver, suggesting that IGFBP5 may play a role in the development of goose fatty liver. Several previous studies have suggested that IGFBP5 is involved in the development of fatty liver disease. First, serum IGFBP5 levels in individuals with fatty liver are higher than those in healthy individuals [14]. Second, IGFBP5 overexpression in mouse liver cells can reduce hepatocyte proliferation, reduce liver enzyme levels in serum, and slow down the progression of fibrosis [23]. Third, HepG2 cells treated with free fatty acids (an in vitro model of NAFLD) had lower IGFBP5 expression levels than those in untreated cells [24] and mice treated with a high-fat diet (in vivo model of NAFLD), also showing lower IGFBP5 expression levels than in mice fed a normal diet [24]. Finally, overexpression of IGFBP5 reduced the levels of lipogenesis-related proteins (FASN, SREBP-1c, and ACC1), elevated the expression of fatty acid β-oxidation-related genes (PPARα, CPT1A, and ACOX1), and decreased intracellular lipid droplets in liver cells [24]. Consistently, the low expression of IGFBP5 in goose fatty liver versus normal liver suggests that IGFBP5 may play an important role in the development of goose fatty liver (e.g., the lower expression of IGFBP5 may promote lipid accumulation in the goose liver).
It is noteworthy that the lower expression of IGFBP5 in goose fatty liver vs. normal liver seems to be inconsistent with the higher level of serum IGFBP5 in humans with fatty liver vs. healthy individuals [14]. This inconsistency may reflect the uniqueness owned by geese. NAFLD contains a sequence of pathological symptoms ranging from simple steatosis, steatohepatitis to fibrosis, or even cirrhosis. In humans, severe hepatic steatosis is generally accompanied with pathological symptoms. In contrast, goose fatty liver, which contains a large amount of fat, is generally physiological, suggesting that goose may have a unique protective mechanism. This uniqueness may result from a long-term evolution of its ancestor as a migratory bird. For example, some new regulatory mechanisms (e.g., new transcription factors, cis-elements, or epigenetic modifications) may emerge and play a role in the development of goose fatty liver during evolution. In terms of this point, goose fatty liver may serve as a NAFLD model for discovering the protective mechanism, which might provide new ideas for treating or preventing NAFLD in humans and other animals.
Inhibition of IGFBP5 expression in goose fatty liver vs. normal liver could be partially attributed to the effects of the factors associated with NAFLD. For example, hyperinsulinemia is usually associated with NAFLD, while it is known that insulin can suppress the expression of IGFBP5 in the liver. On the other hand, IGFBP5 overexpression can reduce hepatic steatosis, and thus lower expression of IGFBP5 may promote goose fatty liver. Additionally, this study showed that IGFBP5 overexpression could activate the MAPK pathway, leading to incidence of inflammation, and thus lower expression of IGFBP5 might contribute to the anti-inflammation mechanism of goose fatty liver. In summary, there may be a cycle from IGFBP5 expression to NAFLD. However, this speculation needs to be validated by further investigation.
To address the mechanism by which IGFBP5 participates in goose fatty liver, we performed transcriptome analysis of goose liver cells transfected with an IGFBP5 overexpression vector vs. an empty vector. The data suggested that IGFBP5 functions mainly through enriched KEGG pathways, which were identified in this study. The MAPK signaling pathway was investigated, and its involvement in the development of goose fatty liver was demonstrated by p38 MAPK immunoblotting analysis, in which IGFBP5 overexpression in goose primary hepatocytes induced the protein expression of total p38 MAPK. In contrast, the suppression of IGFBP5 was accompanied by the inhibition of p38 MAPK expression in goose fatty liver. The other MAPK signaling pathways including the ERK and JNK MAPK signaling pathways might also be regulated by IGFBP5 as transcriptome analysis of goose primary hepatocytes transfected with IGFBP5 overexpression vector and empty vector identified some DEGs in these pathways. For example, KDR in the ERKpathway [25,26] and TNC in the JNK/c-Jun pathway [27]. Unfortunately, the antibodies to goose ERK and JNK MAPK proteins did not work well as there were no expected immunoblots according to their molecular weights. The MAPK signaling pathway regulates cell growth, differentiation, and apoptotic cell death [28]. Therefore, the suppression of IGFBP5 in goose fatty liver may inhibit apoptosis [29] through the MAPK pathway. Further, we identified the focal adhesion, ECM-receptor interaction, and regulation of the actin cytoskeleton pathways, which are associated with inflammation and fibrosis in the liver [30][31][32]. Thus, downregulation of IGFBP5 expression probably prevents the occurrence of hepatocellular carcinoma [29], inflammation, and fibrosis in goose fatty livers. In addition, the AKT and AMP-activated protein kinase (AMPK) pathways and intracellular signaling pathways were identified, which can be activated by IGFBP5 [13]. In HepG2 cells treated with free fatty acids, IGFBP5 overexpression promotes glucose uptake and glycogenesis and suppresses lipid accumulation by activating the IRS1/AKT and AMPK pathways [24]. Moreover, the GnRH signaling pathway was also identified, which has been reported to be involved in liver fibrosis and cirrhosis modulated by melatonin [30,33,34]. Melatonin can inhibit biliary hyperplasia in cholestatic livers by interacting with GnRH receptors [30,33].
The regulation of the enriched pathways identified in this study may be mediated by the IGF/IGFR system. Some studies have shown that IGFBP family members can bind to IGF1/2 and participate in a variety of cellular processes, including cell growth and proliferation, migration, cell survival, apoptosis inhibition, and protein synthesis by modulating the biological availability of IGF1/2 [3]. However, IGF/IGFR-independent mechanisms should not be neglected. Several lines of evidence indicate that IGFBP5 plays a role in IGF-dependent or IGF-independent manners [13]. Whether the influence of IGFBP5 on these pathways in the development of goose fatty liver is IGF-dependent or -independent warrants further investigation.
It is worth noting that previous studies have shown some contradictory findings. For example, some studies have shown that IGFBP5 can inhibit cell growth in breast cancer, but others have indicated that IGFBP5 can promote cell growth [35]. Similarly, overexpression of IGFBP5 is associated with reduced fibrosis in the mouse liver [24]; however, the present study showed that downregulation of IGFBP5 was not associated with increased fibrosis in goose fatty liver. These contradictions may be due to the differences in animal species, cell types, or the conditions where IGFBP5 functions, which require further investigation. In addition, although there are some disagreements in the literature regarding the role of IGFBP5, from the results of this study, it can be concluded that IGFBP5 participates in biological processes through the pathways identified in this study.

Conclusions
This study obtained full-length mRNA and genomic sequences of the goose IGFBP5 gene, demonstrated that suppressed expression of IGFBP5 was accompanied by a decrease in total p38 MAPK protein expression in goose fatty liver, and revealed several pathways that could mediate the biological role of IGFBP5 in goose liver cells. These findings provide a foundation to address the role of IGFBP5 and related mechanisms in the development of goose fatty liver.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10.3390/ agriculture12030347/s1, Figure S1: The gel image showing exclusion of goose genomic DNA contamination from 5 -RACE assay. Figure

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.