Liver Expression of IGF2 and Related Proteins in ZBED6 Gene-Edited Pig by RNA-Seq

Simple Summary Zinc finger BED-type containing 6 (ZBED6), as a regulatory factor, has different regulatory mechanisms in animal development. The intron of insulin-like growth factor 2 (IGF2) regulates the development of animal muscle and adipose by combining with the binding site of ZBED6. As a member of the insulin-like growth factor family, IGF2 plays an important role in embryonic growth and development, cell proliferation, muscle growth and genome imprinting. In order to further study the regulatory mechanism of ZBED6 on IGF2, we detected the expression of IGF2 and related genes in ZBED6 single allele knockout (ZBED6-SKO) pig tissues and analyzed differently expressed genes of the transcriptome of ZBED6-SKO pig liver. The results showed that the partial knockout of ZBED6 could affect the secretion of IGF2 in pig liver but had no significant difference at the protein level. This research provides a new idea for the interaction between IGF2 and ZBED6. Abstract Zinc finger BED-type containing 6 (ZBED6), a highly conservative transcription factor of placental mammals, has conservative interaction of insulin-like growth factor 2 (IGF2) based on the 16 bp binding sites of ZBED6 on the IGF2 sequence. IGF2 is related to embryo growth and cell proliferation. At the same time, its functions in muscle and adipose in mammals have been widely mentioned in recent studies. To further investigate the mechanism of ZBED6 on IGF2, we detected the expression of IGF2 and related genes in ZBED6 single allele knockout (ZBED6-SKO) pig tissues and analyzed the transcriptome of ZBED6-SKO pig liver. Through RNA-seq, we captured nine up-regulated genes and eight down-regulated genes which related to lipid metabolism. The results showed that the mRNA of IGF2 had an upward trend after the partial knockout of ZBED6 in liver and had no significant difference in protein expression of IGF2. In summary, ZBED6-SKO could affect the secretion of IGF2 in pig liver and its own lipid metabolism. Our research has provided basic information for revealing the regulatory mechanism of the interaction between ZBED6 and IGF2 in mammals.


Introduction
Zinc finger BED-type containing 6 (ZBED6) gene is a transcription regulatory factor newly discovered in recent years [1], which is closely related to the development of muscle and adipose and could regulate the expression of various important function genes. It is unique to placental mammals and the sequence is highly conserved [2]. Insulin-like growth factor 2 (IGF2) is a type of cell-like regulatory factor with insulin-like structure and widely expresses in various tissues during the embryonic stage. In adulthood, the expression of IGF2 is significantly suppressed; it mainly synthesizes in liver and transports to various parts of body. IGF2 has many conservative functions in mammals [3][4][5][6][7]. Its functions are mainly through the combination with insulin-like growth factor 1 receptor (IGF1R) and insulin-like growth factor 2 receptor (IGF2R), or it can compete with insulin-like growth factor binding protein (IGFBP). IGFBPs regulate the circulating levels, half-life and activity of IGF2 in target tissues by competitively binding [1,8]. The stable expression of IGF2 is of great significance for animal growth and health. IGF2 disorders result in various diseases such as growth retardation, diabetes, neurodegenerative diseases, atherosclerosis, osteoporosis and cancer [1]. In animal production, the abnormal expression of IGF2 is usually associated with animal production performance [9,10]. The transcriptional regulation of IGF2 depends on the recognition site specific to the transcription factor, so the genetic variation of IGF2 can regulate the growth and development of animals by increasing or decreasing its corresponding transcription factor binding site. ZBED6 could combine the GCTCG located in the third intron segment of the IGF2, thereby affecting the development of muscle and adipose [7,11]. The mechanism was not only verified in mice, but the modification of the ZBED6 binding site in the third intron of pigs' IGF2 by CRISPER/Cas9 also produced consistent results [12].
Currently, most studies focus on the downstream regulation of IGF2, while there is little research on the mechanism of the upstream regulation. In this study, we selected four ZBED6 single allele knockout (ZBED6-SKO) pigs and three wild type (WT) pigs and constructed the expression patterns of IGF2 and related genes in the liver, spleen, heart, semitendinosus and semimembranosus. The RNA-seq was used to analyze the expression of IGF2 and related genes in the liver. Quantitative real-time PCR (qRT-PCR) and western blotting were used to quantified IGF2 and related binding proteins in the tissues of ZBED6-SKO pigs and WT pigs, and hematological and blood biochemistry indices were detected. We investigated the effect of ZBED6 knockout on liver IGF2 metabolism, expected to provide novel insight for understanding the function of ZBED6 and IGF2 interaction and provide the theoretical basis for lean pig molecular breeding.

Materials and Methods
All experimental procedures were performed in accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals approved by the State Council of the People's Republic of China. The study was approved by the Institutional Animal Care and Use Committee of Northwest A&F University (permit number: NWAFAC1019).

RNA Isolation and Library Preparation
Total RNA was isolated via the RNAiso plus protocol (Takara Biomedical Technology Co., Ltd., Beijing, China). RNA quality were evaluated using Qubit2.0 RNA (Life, USA) and a cDNA library prepared using a Hieff NGS™ MaxUp Dual-mode mRNA Library Prep Kit for Illumina ® (YIASEN Biotech Co., Ltd., Shanghai, China). Libraries were quantified (Qubit2.0 DNA Assay Kit, YIASEN Biotech Co., Ltd., Shanghai, China). Among them, only the liver was submitted for sequencing (Illumina Xten, San Diego, CA, USA).

Primers Design and qRT-PCR
qRT-PCR was performed to verify the RNA-seq expression pattern of IGF1, IGF2, IGF1R, IGF2R, IGFBP1, IGFBP2, IGFBP3, IGFBP4, IGFBP5, IGFBP6 and IGFBP7. ACTB was used as the housekeeping gene (Table 1). 1000 ng total RNA was performed to reverse transcription. SYBR green was utilized in qRT-PCR (Y480 Real-Time PCR Detection System, Roche) detection (Takara, Japan). Primer concentrations were 10 nM. Amplification protocol was 95 • C for 30 s, 50 cycles at 95 • C to denature and 60 • C for 30 s to anneal. Melt curve analysis was from 55 to 95 • C, with an increment of 0.5 • C every 5 s. Samples were run in triplicate. All data were normalized to ACTB and calculated with the 2 −∆∆Ct method [13].

ZBED6 SKO Efficiency in RNA Level
ZBED6 gene-edited pig were disrupted by CRISPR/Cas9. The gene editing pig was SKO type, and it could express two types of transcripts, one was normal type and another was KO transcript, which was a frame shift mutation caused by 1 bp deletion (chro9: 64521346del "T"). The WT pig only has normal type. Two pairs of primers were the same for the amplification of the gene-editor locus in the WT pig. However, for the ZBED6-SKO pig, the primer ZBED6-1F/R was used to detect normal transcripts and ZBED6-KO transcripts, and the primer ZBED6-2F/R was used to detect normal transcripts. Two pairs of primers were referenced to the wild type respectively, and the different value between the red primer and green primer for the qRT-PCR of the ZBED6-SKO pig could be seen as KO efficiency. To ensure the accuracy of the amplification, the green primer was set as a primer amplification refractory mutation system (ARMS), and high specific AceTaq ® DNA Polymerase (Vazyme Biotech Co., Ltd., Nanjing, China) was used to enhance the reliability of the ARMS (Figure 1).

RNA-seq Analysis and Statistical
FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used for quality control and filtering of raw data. Trimmomatic (http://www.usadellab.org/cms/page=trimmomatic) was used to remove reads containing adapter, reads containing ploy-N and low-quality reads to get clean reads. The clean reads were mapped to the reference genome (Sscrofa11.1) using HISAT2 tools (https://ccb.jhu.edu/software/hisat2/index.shtml). Statistics of the comparison results were produced using RSeQC (http://rseqc.sourceforge.net/). To perform differential expression analysis, differentially expressed genes (DEGs) were analyzed using the R software package DESeq according to the parameters fold change >2 and false discovery rates (FDR) <0.05. Finally, the results of expression difference analysis were visualized. Gene ontology (GO) enrichment analysis of up-regulated and down-regulated genes was performed using the topGO package (http://www.geneontology.org), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (https://www.kegg.jp/) and euKaryotic Ortholog Groups (KOG) classification enrichment analysis (https://www.ncbi.nlm.nih.gov/COG/) were performed using the clusterProfiler package to determine the biochemical metabolic pathways and signal transduction pathways in which DEGs are mainly involved. In general, when the p value < 0.05, it is considered that the function is significantly enriched. DEGs with significant levels of FDR were used for GO and KEGG enrichment analysis. Moreover, the pig KEGG database was used as reference data sets for gene set enrichment analysis (GSEA) to get differently expressed gene sets. GSEA was performed using the software GSEA v4.0.3 (www.broadinstitute.org/gsea), and 333 pathways of pigs were selected in the KEGG database as the reference gene sets (https://www.kegg.jp/). The three filtering criteria for differential gene sets selection were absolute value of normalized enrichment score (NES) > 1, p value < 0.05 and FDR < 0.25. The chi-square test and analysis of variance (ANOVA) were used with SPSS software 18.0 (IBM, Armonk, NY, USA). Statistical testing was carried on the records (p < 0.05 and p < 0.01).
The high-throughput sequencing of the transcriptome sequence data has been saved in the NCBI Sequence Reading Archive (SRA), where the accession number of the pig liver samples is PRJNA657644 (SRR12480633-SRR12480639). Animals 2020, 10, 2184 7 of 18

Detection of ZBED6 Gene SKO Efficiency in Pig Tissues and Blood Routine Examinations
At DNA level, the genotypes of seven pig samples were identified through sanger sequencing, No.65, No.68 and No.103 were WT, and No.58, No.64 and No.65 and No.90 were ZBED6-SKO; the genotypes of all samples were the same as the record. At RNA level, a two-step qRT-PCR method was used to identify the partial knockout efficiency in five tissues, including liver (57.82%), spleen (44.06%), heart (59.71%), semitendinosus (76.15%) and semimembranosus (63.38%). Compared with the WT group, the mRNA of the ZBED6 gene was significantly down-regulated in all tissues ( Figure 1). In the blood routine examination, the KO group was higher than the WT group for the level of MONO# and MONO%, (p < 0.05) (Table S1). No difference was found in the other indices between the three groups (p > 0.05) (Table S2).

Expression of IGF2 and Related Genes in ZBED6-SKO Pig Tissues
qRT-PCR was performed to detect IGF2 and related genes in the studied tissues. There was no significant difference in IGF1, IGF2, IGF1R, IGF2R, IGFBP1, IGFBP2, IGFBP3, IGFBP4, IGFBP5, IGFBP6 and IGFBP7. It is worth noting that the expression of the IGF2 gene in the liver of ZBED6-SKO pigs had a certain up-regulation (p = 0.06) (Figure 1). The expression levels of IGF2 and related genes were analyzed via liver transcriptome sequencing, and no significant difference was found in related genes. Moreover, IGFBP6 is not expressed in the liver, which can be mutually confirmed with the results of qRT-PCR ( Figure 2). The total expression of the IGFBP family was analyzed and found that the results were also not significant between the ZBED-SKO group and the WT group. In the analysis of the expression of insulin receptors (INSR), it was found that the expression of INSR increased nearly twice in the ZBED6-SKO group (Figure 3). At protein level, the results showed that there is no significant difference in ZBED6-SKO pig tissues. Moreover, there was a large inter-individual difference of IGF2 ( Figure 1).

Quality Evaluation of Sequencing Data
For the RNA-seq, the total reads count, total bases count and average read length of the samples are shown in Table 2. Calculation of the Q scores revealed that the quality of the sequencing data met the standards for further analysis. In addition, the GC content distribution was homogeneous. Table  2 indicates that the sequencing data were well filtered for further research.

Quality Evaluation of Sequencing Data
For the RNA-seq, the total reads count, total bases count and average read length of the samples are shown in Table 2. Calculation of the Q scores revealed that the quality of the sequencing data met the standards for further analysis. In addition, the GC content distribution was homogeneous. Table 2 indicates that the sequencing data were well filtered for further research.

DEGs and Gene Enrichment in ZBED6-SKO Pig Livers
To get further insight into the effect of ZBED6-SKO pig livers, differential expression gene analysis and functional enrichment analysis were performed between ZBED6-SKO and WT pig livers. Using the p value (p < 0.05) and FDR value (FDR < 0.05) as the different standards for analysis, we obtained 139 up-regulated genes and 97 down-regulated genes as well as 9 up-regulated genes and 8 down-regulated genes respectively (Table 3 and Figure 4). A total of 26 significant differential expression pathways were obtained, including 10 up-regulated pathways: lysosome, phospholipase D signaling pathway, glycerophospholipid metabolism, phosphatidylinositol signaling system, systemic lupus erythematosus, type I diabetes mellitus, ether lipid metabolism, primary immunodeficiency, natural killer cell mediated cytotoxicity and Th17 cell differentiation. There were 16 down-regulated pathways: nicotine addiction, P53 signaling pathway, ferroptosis, endometrial cancer, oxidative phosphorylation, cysteine and methionine metabolism, alpha-linolenic acid metabolism, amino sugar and nucleotide sugar metabolism, pancreatic cancer, chronic myeloid leukemia, carbon metabolism, bacterial invasion of epithelial cells, metabolism of xenobiotics by cytochrome P450, renal cell carcinoma, small cell lung cancer and rig-I-like receptor signaling pathway ( Figure 5). Among them, phospholipase D signaling pathway, glycerophospholipid metabolism, phosphatidylinositol signaling system and ether lipid metabolism were all related to lipid metabolism. It is worth mentioning that when the standard of DEGs in GO and KEGG analysis is reduced from FDR to P value, the DEGs are significantly enriched to the gene sets and signaling pathways including fatty acid biosynthesis, TNF signaling pathway, ATP binding, transcription export complex, T cell differentiation in thymus, response to cytokine, I-kappaB kinase/NF-kappaB signaling, skeletal muscle cell differentiation and negative regulation of transcription from RNA polymerase II promoter, some of which are related to fat development (Table S3). This is consistent with the results of GSEA analysis.

Discussion
High sequence conservation of the ZBED6 takes place during the development of placental mammals [1]. The ZBED6-IGF2 interaction is also conservative and based on the 16 bp binding sites of ZBED6 located in the IGF2 sequence [1,10,12]. A mountain of research has shown that the ZBED6-IGF2 axis has an important role in the development of muscle and other tissues [12,14], and the level of IGF2 in the blood was the main focus of attention in their studies. Although liver is the main organ secreting IGF2, there have been few reports about the expression of IGF2 and IGFBPs in the liver. Therefore, the purpose of this study was to explore the expression pattern of IGF2 and their binding proteins in ZBED6 gene-edited pig, to provide a reasonable explanation of the functions of the ZBED6-IGF2 axis in liver, muscle and other tissues.
At RNA level, we analyzed the partial knockout efficiency of the ZBED6 gene and found that the partial knockout efficiency of ZBED6 in different tissues ranged from 44.06% to 76.15%. Therefore, we speculate that this is due to the existence of some compensation mechanisms between different tissues after ZBED6 knockout [15]. As a member of the insulin-like growth factor family, IGF2 is widely involved in many physiological and metabolic processes in the body and plays an important role in cancer development, neuromodulation, diseases of glucose metabolism, osteoporosis, muscle development and fat deposition [16,17]. Due to the partial knockout of ZBED6, IGF2 mRNA was up-regulated in the liver, but this was not statistically significant (p = 0.06). It may be caused by the gradual decline in the expression of IGF2 in adulthood and the increase in degree of methylation. In the other tissues, the expression of IGF2 was not significantly different between the ZBED6-SKO group and the WT group, which was caused by the lower expression level of IGF2 than liver. However, the expression of IGF2 in the studied tissues was not changed significantly at the protein level, which may be due to the decrease of the total level of IGFBP expression in the liver. This phenomenon could have a certain relationship with the maintenance of IGF2 steady state. Similar to our study, Younis found that when IGF2 expression was up to 100 fold, no significant changes were found in some phenotypes of the male samples [7]. This could be related to some unknown mechanisms or regulatory factors on IGF2 regulation. In addition, even if we chose half-sibling or full-sibling individuals, there are still large differences in the expression of IGF2 at the protein level of different individuals, which may be affected by individual differences and limited by current statistical methods. IGF2 enters the various tissues with blood circulation to perform corresponding functions. IGF2 mainly regulates its circulating level, half-life and activity in target tissues by combining with IGF1R and IGF2R or binding to IGFBPs [18][19][20]. IGF2 has an important relationship with insulin; IGF pathway signals could be mediated through IGF receptors or insulin receptors [21]. In this study, RNA-seq analysis of IGF2 receptors (IGF2R) and INSR found that IGF2, IGF2R and INSR were all up-regulated in the liver, but the IGFBPs were not up-regulated. Therefore, we speculate that the liver can use nearby free IGF2 directly and act through IGF2 receptors instead of relying on IGFBPs via autocrine patterns.
There were 17 differently expressed genes in knockout pig liver based on RNA-seq data (FDR < 0.05). Co-expression of SUMF1 and sulfatase will result in a significant synergistic increase in enzyme activity, which can be used to treat polysulfatase deficiency [22]. NIPSNAP1 activates and regulates central sensitization through extracellular signal-regulated kinase (ERK), thereby promoting the occurrence of inflammatory pain. It shows that NIPSNAP1 can be used as a new therapeutic target for pathological inflammatory pain [23]. SEC61A2 mediates the translocation and insertion of membrane proteins (including potassium channels) in the endoplasmic reticulum [24]. In the central nervous system, VCAM1 regulates myelination of oligodendrocytes, and VCAM1 knockout mice have reduced myelin thickness [25]. ABCA7 can transport phospholipids, and the order of transporting phospholipids is phosphatidylcholine (PC)≥lysoPC>sphingomyelin (SM) = phosphatidylethanolamine (PE) [26]. ANKS4B plays a key role in intestinal epithelial cell microvilli assembly during intestinal cell differentiation [27]. AQP1 is a water-selective transport protein affecting the permeability of cell membranes [28]. CatSpere is a subunit of CatSper. The Catsper channel is unique to sperm and is essential for the hyperactivity of sperm flagella [29]. EMR4 is a novel epidermal growth factor (EGF)-TM7 molecule [30]. ZSWIM3 is the zinc finger chelation domain of SWIM, which could play a role in DNA binding and protein binding interactions [31]. Among them, VCAM1 and ABCA7 are both involved in lipid production and transportation [25,26]. ZBED6 can regulate cell proliferation, apoptosis, cell cycle, cell aggregation and neuronal differentiation [2,11,14], which could be related to the effect of VCAM on cell proliferation and thus delaying the differentiation of oligodendrocytes in vivo. In addition, ZBED6 also affects cell-to-cell N-cadherin connection, which participates in the exchange of materials and signal transduction with the external environment [32]. This might be related to the ability of DEGs CRACR2A and USP6NL to control endocytosis and signal transduction [33,34].
Significantly, expression genes (FDR <0.05) were used to evaluate gene functions in ZBED6 SKO pig liver, and no GO, KEGG or KOG signal pathways being enriched. GSEA can reveal many common biological pathways when single gene analysis does not find that the experimental group is associated with the control group. GSEA considers all of the expressed genes in an experiment, instead of selected genes with significant differential expression, which can preserve gene-gene correlations. This method uses the consistency of gene expression to obtain biological information in the data [35]. Among all the up-regulated pathways, there are the phospholipase D signaling pathway, glycerophospholipid metabolism, phosphatidylinositol signaling system and ether lipid metabolism related to lipid development and metabolic regulation, which shows a certain relationship with the function of IGF2. Phospholipase D (PLD) is a member of the superfamily of phospholipases, and phospholipases are essential for intracellular and extracellular signal transduction. Phosphatidic acid is the main metabolite of PLD in cells, which is a precursor of many lipids in the de novo pathway in cells. Moreover, phosphatidic acid can also be converted into two other biologically active lipid molecules, diacylglycerol (DAG) and lysophosphatidic acid (LPA) [36]. Glycerophospholipid is a kind of phospholipid and its metabolism has been found to play a role in the treatment of non-alcoholic fatty liver by drugs [37]. The phosphatidylinositol signaling system is a complex cell regulatory system composed of enzymes, phospholipid messengers and their binding proteins. In this pathway, extracellular signal molecules bind to G-protein-coupled receptors on the cell surface, and phospholipase C hydrolyzes 4,5-bisphosphatidylinositol (PIP2) on the plasma membrane into inositol-1,4,5-triphosphosate (IP3) and diacylglycerol (DG); the extracellular signal is converted into an intracellular signal, which reflects the importance of phosphatidylinositol metabolism in cell regulation [38]. Ether lipids are a unique class of glycerophospholipids, and their tissue distribution will be different. The highest levels of blood lipids are found in the brain, heart, spleen and white blood cells, while the intracellular lipid content in the liver is very low [39]. Mice and humans lacking ether lipid usually show myelination defects in the central nervous system and peripheral nervous system [40]. The myelin sheath is a layer of fatty tissue wrapped around the axons of certain neurons. It has an insulating effect and improves the conduction velocity of nerve impulses, as well as protecting the axons [41]. Since we use FDR < 0.05 as the standard to obtain fewer DEGs, the standard for DEGs is reduced from the FDR value to the P value for GO and KEGG analysis. The results showed that DEGs were enriched in fatty acid biosynthesis and skeletal muscle cell differentiation. Moreover, in the detection of blood biochemical indices, the KO group was lower than the WT group for the level of triglyceride (p < 0.05), which can be mutually confirmed with the results of GSEA.
In this study, we found that the partial knockout efficiency of ZBED6 in different tissues of ZBED6-SKO pigs was from 44.06% to 76.15%. The partial knockout of ZBED6 increased the expression of IGF2 in ZBED6-SKO pigs, but it did not reach a significant level (p = 0.06). This could be due to the gradual decline in IGF2 expression in adulthood and the increase in methylation. At the protein level, there was a large inter-individual variation in the expression of IGF2, and the results between the ZBED6-SKO group and WT group were not significant, which could be related to the maintenance of IGF2 metabolism. Based on RNA-seq data, we performed gene enrichment analysis and GSEA in the liver of ZBED6-SKO pigs, and the results showed that the enriched pathways and gene sets were related to lipid development and metabolic regulation, which indicates a certain relationship with the function of IGF2 ( Figure 6).

Conclusions
The DEGs in ZBED6-SKO pig liver is related to lipid metabolism, production and transportation. The expression of ZBED6 is down-regulated in ZBED6-SKO pigs in five tissues and the expression of IGF2 is up-regulated in the liver, indicating that the partial knockout of ZBED6 has a certain relationship with the secretion of IGF2 in the liver and its own lipid metabolism. At the level of protein, the expression of IGF2 is not significantly different in ZBED6-SKO pigs, which could be related to IGF2 metabolism.
Supplementary Materials: The following are available online at www.mdpi.com/2076-2615/10/11/2184/s1, Table S1: Blood physiological index of ZBED6 KO pig and wild type, Table S2: Serum biochemical index of ZBED6 KO pig and wild type, Table S3: GO and KEGG pathways with significant p value; Figure S1: sequencing of ZBED6-KO locus and IGF2 3072 locus.

Conclusions
The DEGs in ZBED6-SKO pig liver is related to lipid metabolism, production and transportation. The expression of ZBED6 is down-regulated in ZBED6-SKO pigs in five tissues and the expression of IGF2 is up-regulated in the liver, indicating that the partial knockout of ZBED6 has a certain relationship with the secretion of IGF2 in the liver and its own lipid metabolism. At the level of protein, the expression of IGF2 is not significantly different in ZBED6-SKO pigs, which could be related to IGF2 metabolism.