Transcriptome analysis alk-SMase knockout mice reveals the effect of alkaline sphingomyelinase on liver.

Alkaline sphingomyelinase (alk-SMase) is phospholipase that creates ceramides and inactivates platelet activating factors during metabolism, and is linked to digestion and cancer prevention. There have been few studies that completely investigate the linked function and identify the genes related to alk-SMase. In this work, RNA sequencing was performed to investigate the function of alk-SMase. Using RNA-seq data, we discovered 95 differentially expressed genes in the liver of wild type (WT) mice and alk-SMase (gene NPP7) knockout (KO) mice. Differentially expressed genes were functionally associated with the inflammatory response, steroid metabolic process and TNF signaling pathway related terms, regulation of carbohydrate biosynthetic process, and Cytochrome P450-arranged by substrate type using Metascape, according to the results of gene ontology functional enrichment analysis. In addition, an integrated PPI and KEGG network was used to investigate the relationship of differentially expressed genes, It was discovered that one module, which contained 21 nodes and 23 interactions, was significantly related to Cytochrome P450 family, which as mediator of phospholipase. To corroborate the RNAseq results, we identified 12 important genes with high expression levels and significant differences in RNAseq for quantitative real-time polymerase chain reaction (qPCR) verification; 7 genes showed difference significantly (P < 0.05). Our research is the first to conduct a comprehensive genome-wide analysis of the alk-SMase knockout model. Alk-SMase is involved in sphingomyelin hydrolysis. In liver tissues lacking alk-SMase expression, the expression levels of seven genes, including Cyp4a14 and Cyp4a10, were altered.

alkaline sphingomyelinase

Biological general repository for interaction data sets

Molecular Complex Detection

Gene Ontology

Kyoto encyclopedia of genes and genomes

knockout

wild type

Introduction

Intestinal alk-SMase has phospholipase C activity, which can hydrolyze SM to ceramide, inactivate platelet activating factor (PAF) and limit the generation of lysophosphatidic acid (LPA), all of which are linked to colon cancer inhibition [33,35]. Alk-SMase was discovered in 1969 [22], which is specifically expressed in human liver and intestinal mucosa [23]. Alk-SMase in the intestine have been extensively studied, but few have been studied in liver. Alk-SMase as a novel member of the phosphodiesterase (NPP) family and called NPP7 [10]. Alk-SMase is similar to other NPP members and is anchored to the cell membrane surface with a short hydrophobic domain. The rest of the enzyme that includes the catalytic domain is exposed outside the cell [15]. Bile salts and trypsin can release and produce this enzyme.

Recently, some important biological effects of alk-SMase was reported, they have purified the protein [6,11] and cloned the gene [10,32], and found essential roles of the enzyme in SM digestion [23], colonic inflammation [3], cell proliferation [16] and cholesterol absorption [13]. The liver is known to be an essential organ of lipid metabolism, such as fatty acid catabolism, cholesterol metabolism, and phospholipid metabolism. Studies on the metabolism of phospholipids in the liver have previously focused on phosphatidylcholine rather than sphingomyelin, but the studies of sphingomyelin has increased in recent decades. The liver is with relatively high levels of sphingomyelin.

Therefore, human bile alkaline phosphatase research appears to be lagging behind. The main challenges are a lack of animal models and cell lines that express alk-SMase in the liver, which has been shown to reduce carcinogenesis. We discovered phenotypic differences between alk-SMase KO mice and WT mice in these investigations attributable to enzyme deletion. However, the role of associated gene connections and the transcriptome mechanism remain unknown.

Transcriptome analysis provides a comprehensive picture of gene expression data, allowing researchers to better understand how normal cell phenotypes transform and progress under various environments [14]. However, the molecular networks and mechanisms of these interacted genes in alk-SMase KO mice remain mainly unknown. Here, we have performed the first RNA sequencing on alk-SMase knockout and wild-type mice in liver tissue, which characterized the transcriptional landscape of alk-SMase knockout mice with unprecedented resolution to understand the biological function of alk-SMase. Differentially expressed genes between alk-SMase knockout and wild-type groups were identified by Cuffdiff [31], the functional enrichment analysis of differentially expressed genes and transcripts was conducted to understand the biological function, and further explored the related genes of alk-SMase from the transcriptome level by constructing gene to gene regulatory network, which merging an integrated KEGG and PPI networks. We looked into the biological relationship to learn more about alk-SMase properties and molecular control mechanism. These findings may provide a theoretical foundation for further research into the phenotypic and other abnormalities seen in alk-SMase knockout mice.

Materials and methods

Animals and samples

We constructed Alk-SMase KO mice using the Cre-LoxP system [28]. Exon 2 of the NPP7 gene on chromosome 11 is flanked by two loxP sites. The Cre recombinase reaction causes a frameshift, which results in the production of a new stop codon and the termination of translation. To obtain an inbred strain, KO mice were kept on a C57BL/6 background and backcrossed eight times. Mollegaard provided the matched C57BL/6 WT mice of the same age and sex. The Animal Ethical and Welfare Committee of Harbin Medical University in Daqing Campus authorized the study's ethics authorization, which is HMUDQ-2015-068. Three matched pairs of liver tissue were taken from three WT mice and three KO animals, respectively.

RNA extraction, sequencing and data processing

Total RNA was extracted with TRizol method according to manufacturer's instructions. RNA samples were deeply sequenced using Illumina Hiseq2500 platform after further cleaning-up using Qiagen RNeasy mini columns, and 50bp single end reads were generated [2,20]. Quality control was performed by FastQC, trim-galore was used to move the adapter and get the clean reads. Bowtie2 [19] and Tophat 2.0.6 [30] was chose to mapping the reads to the mouse reference genome (mm10) from UCSC database，resulting in alignment rate of 96.83%–97.63%. The accepted_hits.bam files were generated by the tophat mapping, which were used to assembled into transcripts using Cufflinks (Version 2.2.1) [18]. After the alignment, we used ANNOVAR to annotate the position on the reads alignment. The annotation priority is the exonic region (UTR5/UTR3, gene upstream/gene downstream) > splicing region > intronic region > gene intergenic region. Then we identified the distribution of reads on the genome. We used Cuffmerge to combine all transcripts of six samples to generate a consensus transcriptome, and then normalized transcript expression levels by Cuffnorm.

Identification of the differentially expressed genes

Based on | log2 (fold change) | > 1and adjusted p value < 0.05, Cuffdiff 2.2.1 was used to identify differentially expressed genes and transcripts in liver tissue across 3 WT mice and 3 KO mice samples. Three wild type samples and three KO samples were applied to identify the genes with significant up regulated gene expression and up regulated gene expression. The heatmaps of the most significantly differentially expression genes and transcripts were generated with R command heatmap2.

Gene ontology and KEGG pathway enrichment analysis

In this work, functional enrichment analysis of the differentially expressed genes between KO mice and WT mice were implemented using by metascape [17], We first identified all statistically enriched terms, including GO and KEGG terms (p-value < 0.05 and FDR <0.01). Furthermore, the important terms of enrichment were hierarchically clustered to form a network map based on the Kappa statistical similarity between genes. Then 0.3 kappa score was used as the threshold to put the tree into term clusters.

Construction of gene regulator network and identification of the significant genes

The protein-protein interaction data was obtained from two PPI database, Biological general repository for interaction data sets (BioGRID) [5]. To get the more comprehensive interaction information, we chose the union data of the two databases as background. And the KEGG network was established based on the relationship between genes in the pathway using our written Perl programs. Eventually, gene regulator network was constructed by merging PPI and KEGG network. The 95 differentially expressed genes were set as seed genes and then mapped into the integrated background network to extract the PPI and KEGG sub-network. The sub-network was composed of these seed genes and the nearby genes within one step distance from these seed genes in integrated background network. The network was constructed using Cytoscape_3.3.0 [26], applied to construct the biological networks. Network module analysis and hub gene identification were subjected from the interaction network by using the plugin MCODE (Molecular Complex Detection) [4] method in Cytoscape_3.3.0.

RT-PCR verification

Total RNA was extracted with Trlzol reagent (Invitrogen) and reverse transcribed with the one-step RT-PCR kit (TransGen Biotech Co,Ltd.,China) according to the manufacturer's instructions. All primers were obtained from invitrogen and designed using reference sequences published by the National Center for Biotechnology Information. Briefly, PCR was performed in a reaction volume of 25 μL containing 200 ng cDNA, 12.5 μL SYBR Green RT-PCR Master Mix, and 1.25 μL of each of two primer solutions (10 μM). Quantification of gene expression was performed using the ABI PRISM 7500 Sequence Detection System (Applied Biosystems, Foster City, CA) with SYBR Green (TransGen Biotech Co, Ltd, China). The mRNA levels were determined by qRT-PRC in triplicate for each of the independently prepared RNAs and were normalized to the levels of β-actin expression. Six differential expressed genes were chose for the experiment. T-test was used for significance testing between wild type and KO groups, p-value≤0.05.

Results

Analysis of transcriptome characteristics

The whole transcriptome sequencing was conducted on three matched samples of alk-SMase knockout mice and wild type mice. In total, close to 120 million single-end reads were obtained by HiSeq2500 ultra-high-throughput sequencing systems. In addition, Table 1 listed the read mapping rate in all six samples, which are more than 95%. The quality control of fastq data produced by high throughput sequencers is performed by FastQC software. RNA-seq reads mapped to a reference genome, while Cufflinks was used to assemble these mapped reads into probable transcripts and then generate transcriptional loci. Cuffmerge was chose to combine the transcripts gtf file generated by each Cufflinks into the comprehensive transcripts annotation file. The expression of all genes was normalized by Cuffnorm were in agreement, with high correlation (Fig. 1 and Supplementary Fig. S3). We annotated the position on the reads alignment and then calculated the distribution of reads on the genome (Supplementary Fig. S1). Although gene expression patterns between the two group of samples were very similar, numerous differences were observed.

Table 1

The overall read mapping rate.

Sample	Mapping rate	Input reads	Mapped reads
K82G	97.39%	22,444,506	21,859,552
K84G	97.27%	21,035,226	20,461,432
K85G	97.05%	20,771,090	20,159,341
W81G	96.83%	16,845,418	16,311,244
W83G	97.63%	17,074,312	16,669,266
W90G	97.09%	17,271,567	16,768,238

Fig. 1

The genes global statistics of gene expression in KO samples and WT samples. (A)Count vs dispersion plot by condition in KO samples and WT samples. (B) Gene expression distributions density graph of different samples. (C) Boxplot of gene expression in KO samples and WT samples. (D) The correlation graph of two types of different samples.

Identification of differentially expressed genes and transcripts in KO mice

We performed transcriptome analyses to identify the differentially expressed and functional pathways to characterize differentially expressed genes and transcripts by bioinformatics methods. Cuffdiff and fold change were used to investigate the differentially expressed genes and transcripts. A total of 95 genes were determined differentially expressed in liver tissue of KO mice compared with WT mice, including 58 up-regulated and 39 down-regulated mRNAs. Hierarchical cluster analysis (Fig. 2) showed that liver samples from KO mice were clustered into one cluster, and WT mice were clustered into another cluster, indicating that the expression levels of differentially expressed mRNA in KO mice. As shown on the volcano plot (Supplementary Fig. S2) [7] a number of genes were affected by gene knockout, it also displayed the results of differentially expressed mRNAs, black dots indicates no statistically significant differences in gene expression and the red dots in the figure represent differentially expressed genes with statistical significance.

Fig. 2

Hierarchical clustering analysis of differentially expressed genes and transcripts in KO samples and WT samples. (A) The heat-map of 95 differentially expressed genes, red and green colors in the heatmap represent induced and repressed genes, respectively. Scale bar denotes the log2 value of fold change. (B)The heatmap of 10 differentially expressed transcripts. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Function enrichment analysis of differentially expressed genes

To explore the potential biological functions of the differentially expressed genes, functional enrichment analysis was performed for the genes using the Metascape functional annotation tool. With a p-value≤0.05 as the threshold, the results indicated that the up regulated genes were mainly involved in inflammatory response associated pathways and biological functions, such as, inflammatory response, metal sequestration by antimicrobial proteins, neutrophil dereanulation inflammatory response, cell chemotaxis, response to peptide, TNF signaling pathway, cellular response to inerleukin-1, regulation of peptidase activity, metal ion homeostasis, response to peptide, steroid metabolic process and so on (Fig. 3A). These down regulated genes were enriched in six functional terms, including monocarboxylic acid metabolic process, Cytochrome P450 - arranged by substrate type, circadian regulation of gene expression, regulation of carbohydrate biosynthetic process, monocarboxylic acid biosynthetic process and carboxylic acid transport (Supplementary Fig. S4). Enrichment Map graph revealed to the association of the significantly enriched pathways and go terms as networks.

Fig. 3

Function enrichment analysis for differentially expressed genes. Function enrichment analyses were performed using Metascape with p < 0.05. (A) Function enrichment analysis for up-regulated expressed genes. (B) Enrichment map revealed the connection of different functional terms, which terms shared the same genes. The color of nodes correspond to different functions. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Identification of significant genes by EMODE in integrated PPI and KEGG network

As a result of the complexity of gene regulation, we expected to find the interaction and regulation between gene and gene in NPP7 knockout mice. In recent years, the establishment of biological networks has become important to the study of systems biology, in this study, we selected an integrated protein-protein interaction and KEGG pathway-based biological network as a background network, which is more accurate and straightforward in demonstrating gene relationships. The mouse protein-protein interaction data sets BioGRID, and pathways that differentially expressed genes were enriched. We mapped 95 differentially expressed seed genes into the integrated background network to extract the PPI and KEGG sub-network. The sub-network was composed of these seed genes and the nearby genes within one step distance from these seed genes in integrated background network, including 516 node sand 857 edges (Fig. 4A), in which some nodes were presented with high degrees (hub nodes) while the majority of the nodes were presented with low degrees. The hub genes in this subnetwork included Cyp2c55, Pklr, Cyp4a14, Cyp4a10, Cyp26a1, Cyp7a1 Cyp2a5, Ppp1r3b, Fasn and Gstt3 (Table 2), they are all down-regulated expressed.

Fig. 4

The subnetwork were constructed by mapping DEGs to PPI and KEGG integrated network. (A) The sub-network was extracted from the PPI-KEGG integrated by one-step neighbor nodes. The size of the point represents the size of the node's degree. (B) Function module identified from the PPI-KEGG integrated network.

Table 2

Hub genes in the PPI and KEGG network.

Gene	Regulation	Degree
Cyp2c55	down	141
Pklr	down	123
Cyp4a14	down	112
Cyp4a10	down	112
Cyp26a1	down	78
Cyp7a1	down	50
Cyp2a5	down	42
Ppp1r3b	down	17
Fasn	down	17
Gstt3	down	16

In light of the genes connected in the network and genes in the same pathway thought to be with similar function and biologically related, dividing these highly connected genes into groups by network analysis may be useful for potential functional processed in a manner complementary to standard differential expression analysis. A plugin MCODE in Cytoscape_3.3.0 was used to further identify the functional modules from the sub-network. The parameters were set degree cutoff = 2 and max.depth = 100 to find functional modules and finally obtained the module with gene number larger than four [1,27] (Fig. 4B). The module contains 21 nodes and 23 edges.

Validation of differentially expressed genes by real-time PCR

To conform the result of RNA-seq data, we used quantitative real-time PCR (qPCR) to examine the expression of 12 important genes with high expression levels and significant differences in RNAseq, including highly expressed transcripts in the sub-network. As shown in Fig. 5, the findings of qPCR experiment suggested that all transcripts were consistent with those from RNA-seq, and 7 genes showed difference significantly by t-test (P < 0.05), such as Cyp4a14 and Cyp4a10 (Fig. 5). However, since differential expression analysis was applied to three biological replicates, respectively, partial transcripts show opposite expression patterns between different repeats, showed with individual differences. Furthermore, despite the fact that the mean of the expressions in the figures displays considerable disparities, the difference was not statistically significant.

Fig. 5

Results of rtPCR performed on RNA extracted from WT mice and NPP7 knockout mice.

Discussion

In this study, we performed an comprehensive analysis to investigate the interaction and connection between alk-SMase (NPP7) and other genes and the functional of alk-SMase (NPP7) in the liver by sequencing the global transcripts expression of liver tissue in alk-SMase knockout mice and wild-type mice.The goal was to figure out how alk-SMase works at a molecular level in various conditions.

They discovered that alk-SMase plays critical roles in SM digestion, colonic inflammation, cell proliferation, and cholesterol absorption, among other biological processes. Colon cancer and colitis were reported to have lower alk-SMase activity. Although the role of alk-SMase in the development of colon cancer has been investigated, the enzyme's biological importance and interactions in the liver have not been thoroughly investigated. Gene expression analysis is a critical tool for determining and understanding the mechanisms underlying illness phenotypes in a variety of conditions [36]. The regular analysis in transcriptome profile is screening differentially expressed genes. 95 differentially expressed genes were identified between the wild type mice and alk-SMase knockout mice, and these genes were related to inflammatory response and TNF signaling pathway, TNF is a multifunctional pro-inflammatory cytokine that has been found to be closely related to lipid metabolism, coagulation and endothelial function. TNF also is an anti-cancer agent [21,25], which is related to cell differentiation [8,9], apoptosis [24] and inflammatory reaction [29,34]. Gstt3, one of the down-regulated expressed genes concentrated in the TNF signal pathway, was discovered to be a gene for glutathione metabolism and glutathione transferase, as well as genes associated to amino acid, protein, lipid, and glucose metabolism. As a result, we deduced that variable Gstt3 expression affects the phenotypic alterations of our alk-SMase knockout mice. The absence of NPP7 caused alterations in gene expression in the intestine and hepatic transcriptomes in gene knockout mice. Changes in intestinal function are caused by the loss of expression of the intestinal alkaline sphingomyelin, which affects bile acid metabolism, bile acid enterohepatic circulation, and changes in the expression of other associated genes in the liver. We analyzed genes differentially expressed in the liver of knockout mice, and the functions were mainly enriched in pathways such as metabolism.

We built a network based on protein-protein interaction and pathway information to validate the putative biological links of DEGs. The sub-network revealed one functional module with a gene number greater than four. For example, in the sub-network (Fig. 4B), several major Cytochrome P450 genes (Cyp2c55PklrCyp4a14Cyp4a10Cyp26a1Cyp7a1Cyp2a5) in model of gene knockout mice and intestine-liver interaction were decreased expression, which corresponded to decreasing levels of macrophages and dendritic cells in the lamina propria [12]. Furthermore, in later tests, we can investigate the function of each module in greater depth. We discovered several genes that have a major role in the function related to liver metabolism through a literature review of the key genes identified by the network.

Conclusions

In conclusion, we studied the genome-wide transcriptome of alk-SMase knockout mice with high resolution. The biological significance of alk-SMase was studied. The findings gave fresh information on the functional modules of alk-SMase that are involved in liver disorders.

Ethics approval and consent to participate

The animal studies were approved by Animal Ethical and Welfare Committee of Harbin Medical University in Daqing Campus (Approval number HMUDQ-2015-068) and this study was carried out in accordance with the guidelines and regulations established by this committee.

Availability of data and material

The complete RNA sequencing data generated and analyzed in this study have been submitted to the Genome Expression Omnibus, the accession number is GSE128349 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE128349). To review GEO accession GSE128349:Go to https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE128349.Enter token ehuzwogerxanjsz into the box.

Funding

This work was supported by grants from the Fundamental Research Funds for the Provincial Universities [grant number DQYWH201602].

Natural Science Foundation of Heilongjiang Province [grant number LH2020C058] and Health Council Research Program of Heilongjiang Province [grant number 2019–075].

Author contributions

P Zhang conceived and designed the experiments. J Zhu and L Wang collected the data. X Li and YQ Cheng performed the analysis. J Zhu and DX Lan participated in the discussion of the algorithm. P Zhang and YC Li edited the manuscript and revised the manuscript. All authors have finalized the draft of the paper.

Declaration of competing interest

The authors declare no competing financial interests.