Renal Transcriptomics Reveals the Carcinogenic Mechanism of Ethyl Carbamate in Musalais

Introduction

Musalais is a traditional fermented wine produced by the Uygur people in southern Xinjiang of China. As a traditional wine with regional characteristics, musalais is protected under the scope of the national intangible cultural heritage.¹ Ethyl carbamate (EC) is naturally produced during the fermentation process of food and alcoholic beverages, especially in wine made from stone fruits.² It has been demonstrated to induce carcinogenesis and was classified as a group 2A carcinogen by the International Agency for Research on Cancer of the World Health Organization,^3–5 resulting in concerns on the safety of EC in fermented wine. Previous studies have demonstrated the carcinogenic potency of EC on lung tumorigenesis.⁶ In addition. EC-induced consequent mutation and DNA adducts were observed in organs such as mammary glands, liver, ovary, heart, and forestomach.⁷ EC is produced during the fermentation process of musalais, and its content increases with the duration of fermentation. In recent decades, substantial effort has focused on fermentation characteristics and production engineering,⁸ but safety control of EC in musalais has attracted minimal attention, and the precise effect of EC in musalais remains to be elucidated.

Transcriptomics, especially gene-array technology, enables the detection of the mRNA levels of multiple genes in parallel. It provides genome-wide information on biological processes and related molecular pathways and is an effective tool for screening compounds in the assessment of toxicological potential. This information can be used to identify biomarkers for compounds and evaluate tissue-specific toxicity, allowing earlier and better prediction of toxicity.⁹ In this work, rats were treated with musalais containing EC at varying contents (0.1, 1, and 10 mg/kg). To investigate the toxicity of EC in musalais, the liver and kidney of the rats were subjected to transcriptomic sequencing (RNA-seq). Differentially expressed genes (DEGs) were identified and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed on these DEGs to examine the important biological functions affected by EC in musalais. Furthermore, the key disease-related biomarkers of EC in musalais were predicted.

Materials and Methods

Ethical Statement

Experiments involving rats were approved by the Institutional Review Board of Wuhan Myhalic Biotechnology Co., Ltd. based on the ethical Guidelines for Animal Care and Use of the Model Animal Research Institute (approve number: HLK-20181209-01).

Animals and Experimental Design

Sixteen healthy Sprague Dawley rats aged 8–10 weeks were purchased from the Hubei Provincial Center for Disease Control and Prevention. The rats were randomly divided into four groups (4 rats per group): control (Con), High, Middle (Mid), and Low. Rats in the Con group were fed normally without any treatment. Rats in the High, Mid, and Low groups were subjected to intragastric administration of musalais containing 10 mg/kg, 1 mg/kg, and 0.1 mg/kg EC, respectively, once a day (2 mL each time) for a week. Musalais was spiked with varying EC concentrations to achieve total content of 0.1, 1, and 10 mg/kg in the Low, Mid, and High groups, respectively. Thereafter, the rats were anesthetized and sacrificed to obtain the liver (L) and kidney (K). The collected tissues were maintained at −80 °C for subsequent sequencing.

Preparation of Musalais Containing EC

EC (Dr. Ehrenstorfer, Germany; purity ≥99.4%) was diluted in musalais and its the concentration was determined by gas chromatography-mass spectrometry using a GCMS-QP2010 Plus apparatus (Shimadzu, Japan). Briefly, 1 μL of sample was injected into a DB-WAX chromatographic column (30 mm × 0.25 mm, 0.25 μm, Agilent, USA). The initial column temperature was 100 °C and increased to 120 °C at a rate of 50 °C/min, then held at 120 °C for 5 min and increased to 200 °C at a rate of 50 °C/min and maintained for 1.5 min. High-purity helium gas (>99.999%) with a flow rate of 1.0 mL/min was used as the carrier gas. Mass spectrometry was performed in selected-ion monitoring mode with a transfer line at 280 °C, electron voltage of 70 eV, and source temperature of 230 °C.

Liver and Kidney Transcriptomic Analysis

Total RNA was extracted from kidney and liver tissues using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and purified using the Total RNA Purification Kit (TRK1001, LC Science, Houston, TX, USA) following the manufacturer’s procedure. RNA content and purity were quantified using a NanoDrop ND 1000 apparatus (NanoDrop, Wilmington, DE, USA). The RNA integrity number was assessed to be greater than 7.0 using the Agilent 2100 Bioanalyzer (Agilent, CA, USA). mRNA was purified from total RNA using poly-T oligo-attached magnetic beads with two rounds of purification. The purified mRNA was fragmented into small pieces using divalent cations at high temperature. Cleaved RNA fragments were reverse-transcribed into cDNA according to the protocol of the mRNA-Seq sample preparation kit (Illumina, San Diego, USA). The average insert size for the final cDNA library was 250–350 bp. Next-generation sequencing was performed on 150-bp paired-end reads in the library using an Illumina Hiseq 4000 sequencing platform at LC Bio (Hangzhou, China).

Data Analysis

Valid data were collected by removing reads including sequencing adaptors and sequencing primers and nucleotide with quality scores (Q) of lower than 20 from raw data obtained using the Illumina paired-end RNA-seq approach. In addition, the content of Q20 (Q ≥ 20) and Q30 (Q ≥ 30) and the GC content of valid data were calculated. All subsequent analyses were carried out using valid data. Valid data were mapped to the UCSC reference genome of rattus norvegicus (http://genome.ucsc.edu/) using the HISAT package¹⁰ and the mapped reads were assembled using StringTie software.¹¹ Perl scripts were used to reconstruct a comprehensive transcriptome, and the fragment per kilobase of exon model per million mapped reads (FPKM) was calculated as the gene expression level for each sample using StringTie. The distribution statistics of the FPKM of genes were visualized as boxplots using R software (version 3.5.2).

Differential Expression Analysis

Differential expression analysis was performed on the FPKM values of the untreated and treated groups using the edgeR package. DEGs were identified using the following criteria: genes between Con and High/Mid/Low groups displaying |log₂(fold change)| ≥ 1 and p < 0.05. The identified DEGs were visualized using volcano plots generated by R software (version 3.5.2). For hierarchical clustering, the data were Z-normalized by gene, and the Euclidean distance was chosen as the similarity to cluster expression profiles.¹² The results of cluster analysis were visualized using R software (version 3.5.2). The Gene Ontology (GO) (http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp/kegg/) databases were used to examine the enriched GO terms and pathway associated with the DEGs. GO terms and pathways with p < 0.05 were considered as statistically significant.

Results

Preliminary Analysis of Transcriptomics Sequencing Data

A total of 16 RNA-seq libraries were constructed using RNA samples from kidney and liver tissues extracted from both treated and untreated rats. The raw reads, raw bases, valid reads, valid bases, Q20, Q30, and GC percentage were collected from each library. As shown in Table 1, Q20 > 99%, Q30 > 97%, and GC content ranged from 48% to 50%. These results demonstrated the high quality of the transcriptomic sequencing data and guaranteed the reliability of subsequent analysis. In addition, the FPKM values for each library are shown in Table 2 and visualized by boxplots (Figure 1). The boxplot for each region corresponds to the maximum, upper quartile, median, lower quartile, and minimum FPKM from top to bottom. The FPKM values of liver tissue were similar among all groups, while there were differences between groups in kidney tissues. The FPKM values in both liver and kidney tissues were similar between parallel groups, indicating the excellent reproducibility of the samples.

Table 1 Sample Details from RNA-Sequencing

Table 2 Distribution Statistics of FPKM Values of Genes

Figure 1 FPKM value for each library. The boxplot for each region corresponds to maximum, upper quartile, median, lower quartile, and minimum FPKM from top to bottom. FPKM, fragments per kilobase of exon model per million mapped reads.

Identification of DEGs

To better understand the effect of EC in musalais, rats were treated with musalais containing EC at varying contents. DEGs were identified in liver and kidney tissues by comparing the FPKM values in the High/Mid/Low groups with those in the Con group, using the criteria |log₂(fold change)| ≥ 1, false discovery rate < 0.05, and statistical difference p < 0.05. DEGs with statistical differences (Supplementary 1) were illustrated using volcano plots (Figure 2, up-regulated genes are highlighted in red, down-regulated genes are highlighted in blue, and genes showing no difference are highlighted in grey) and counted (Figure 3). The kidney tissues in the High group showed the highest number of DEGs, with 435 up-regulated (Red) and 1877 down-regulated (Blue) genes. The complete list of DEGs is shown in Supplementary 2. Hierarchical clustering (unsupervised), performed to analyze the DEGs, clearly described the separation within the shortlisted genes in the respective groups and demonstrated the different patterns in liver (Figure 4) and kidney (Figure 5) tissues between treated and untreated rats.

Figure 2 Volcano plots of the DEGs. Up-regulated genes are shown in red and down-regulated genes are shown in blue. Genes that showed no difference are shown in grey. DEGs, differentially expressed genes.

Figure 3 Number of DEGs. Red represents the number of up-regulated genes and blue represents the number of down-regulated genes. DEGs, differentially expressed genes.

Figure 4 Cluster analysis of the DEGs in liver samples. Orange indicates that gene expression is up-regulated. Blue indicates that gene expression is down-regulated. The color intensity indicates the degree of up- or down-regulation. DEGs, differentially expressed genes.

Figure 5 Cluster analysis of the DEGs in kidney samples. Orange indicates that gene expression is up-regulated. Blue indicates that gene expression is down-regulated. The color intensity indicates the degree of up- or down-regulation. DEGs, differentially expressed genes.

GO Enrichment and KEGG Pathway Enrichment Analysis of DEGs

Comparing the biological function of the DEGs in the High and Mid groups with those in the Low group clearly revealed the effect of high EC content in musalais. The biological functions of a gene can be defined using three GO categories: biological process, molecular function, and cellular component. Figures 6 and 7 show the number of DEGs associated with various GO terms in liver and kidney tissues, respectively. The GO annotations indicated that the key subcategories in each group were similar in the liver and kidney, ie, cytoplasm, nucleus, integral component of membrane, and extracellular exosome.

Figure 6 GO subcategories of the DEGs in liver samples. GO, Gene Ontology; DEGs, differentially expressed genes.

Figure 7 GO subcategories of the DEGs in kidney samples. GO, Gene Ontology; DEGs, differentially expressed genes.

To better understand the biological effect of high EC content in musalais, KEGG pathway enrichment was analyzed based on the DEGs. The details of the predicted pathways are listed in Supplementary 3 and the enrichment plots of the top 20 KEGG pathways are displayed in Figure 8. For liver tissues, compared to the Low group, there was no regular pattern in the predicted KEGG pathways showing statistical differences (p < 0.05) in both High and Mid groups. However, for kidney tissues, all pathways with p < 0.05 were found in the High and Mid groups but not in the Low group. Specifically, these pathways were drug metabolism – cytochrome P450, chemical carcinogenesis, metabolism of xenobiotics by cytochrome P450, Wnt signaling pathway, and p53 signaling pathway, indicating that high EC content in musalais could affect metabolic and biological processes in rats by modulating these pathways. Further Venn diagram analysis demonstrated that in the High group, 11 genes (Atad1, Gsta5, Gsta3, Ugt2a3, Gstp1, Gsta2, Adh1, Gsta1, Gsto2, Mgst1, and RGD1562107) were simultaneously enriched in drug metabolism – cytochrome P450, chemical carcinogenesis, and metabolism of xenobiotics by cytochrome P450 pathways and one gene (Ccnd1) was enriched in both Wnt and p53 signaling pathway (Figure 9A and Table 3). In the Mid group, 10 genes (AABR07014550.1, Mgst1, Aldh1a3, Gstp1, Gsta5, Gsta1, Gstt2, Adh1, and Gsta2) were simultaneously enriched in drug metabolism – cytochrome P450, chemical carcinogenesis, and metabolism of xenobiotics by cytochrome P450 pathways and one gene (Ccnd1) was enriched in both Wnt and p53 signaling pathway (Figure 9B and Table 4). Among these genes, Mgst1, Gstp1, Gsta5, Gsta1, Adh1, Gsta2, and Ccnd1 were simultaneously observed in the High and Mid groups. These results suggested that high EC content in musalais could affect metabolic and biological processes in rats by regulating the expression of Mgst1, Gstp1, Gsta5, Gsta1, Adh1, Gsta2, and Ccnd1.

Figure 8 KEGG enrichment plots for top 20 pathways in each treatment group. KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 9 Comparison of enriched DEGs in the specific KEGG pathways in High (A) and Mid (B) groups of kidney tissues. The numbers of DEGs enriched in the corresponding pathway are indicated. DEGs, differentially expressed genes. KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table 3 Identification of Genes Enriched in Similar Pathways in High Kidney Group

Table 4 Identification of Genes Enriched in Similar Pathways in Mid Kidney Group

Discussion

Musalais is a unique wine brewed by the local people in southern Xinjiang in China. With a history of more than 2000 years, it is an indispensable part of the local culture and has been protected as a form of intangible human cultural heritage since 2007, contributing substantially to the development of local tourism and economy.^1,13,14 Musalais is produced through natural fermentation, resulting in the generation of EC,^15–17 also known as urethane. EC is an ethyl ester of carbamic acid and is mainly synthesized via the metabolism of urea and citrulline during the fermentation process.¹⁸ Accumulating evidence has demonstrated the carcinogenic effect of EC in organs such as the liver, lung, lymph, skin, and breast.^19–21 In addition, Lajovic et al demonstrated the carcinogenesis of vinyl carbamate epoxide, the ultimate carcinogen of EC metabolism,²² through both experimentally and quantum chemical simulations.²³ Thus, many countries including Canada,²⁴ the USA,²⁵ Brazil,²⁶ and South Korea²⁷ have set a maximum residue level for EC in fermented beverages. Despite that high levels of EC are found in many regularly consumed Chinese alcoholic beverages, no regulation has been implemented on EC,²⁸ and safety assessment has not been conducted on EC in alcoholic beverages.

In this work, transcriptome sequencing was performed to analyze the effects of musalais containing varying amounts of EC on gene expression in rat liver and kidney, the major organs involved in metabolism and detoxification. In addition, the biological functions of the identified DEGs were examined. Although plenty of studies demonstrated that alcohol influences gene expression in the liver and kidney,^29,30 the alcohol content in each EC treatment group in the current study was the same, resulting in similar effects on gene expression. Thus, the effect of alcohol can be excluded. GO enrichment analysis based on DEGs demonstrated that the key subcategories in the liver and kidney were similar between groups, ie, cytoplasm, nucleus, integral component of membrane, and extracellular exosome. However, KEGG pathway analysis revealed that drug metabolism – cytochrome P450, chemical carcinogenesis, metabolism of xenobiotics by cytochrome P450, Wnt signaling pathway, and p53 signaling pathway were only found in the High and Mid groups but not in the Low group in kidney tissues. Comparing the biological functions of DEGs in the High and Mid groups to those in the Low group clearly revealed the effect of high EC content in musalais on rats. Cytochrome P450, a heme-mercaptide protein, is involved in the metabolism of endogenous substances such as steroids, acids, prostaglandins, and unsaturated fats and exogenous substances including drugs, natural compounds, and carcinogens.^28,31,32 A defense system centered on the heme protein cytochrome P450 protects organisms from toxic compounds.^33,34 However, cytochrome P450-mediated biotransformation could lead to the metabolic activation of environmental chemicals into active carcinogens.³⁴ For instance, cytochrome P450 promotes the metabolism of aflatoxin B1, which is known as the most potent natural carcinogen to humans.³⁵ Studies on the carcinogenic mechanism of EC have shown that its metabolism in biological systems is mainly related to cytochrome P450. Specifically, EC can be reversibly converted by P450 to ethyl N-hydroxycarbamate, a proximate on ultimate carcinogenic metabolite of ethyl carbamate.³⁶ In addition, EC is oxidized by cytochrome P450 to vinyl-carbamate, followed by the formation of methyl epoxyethylene formate. This epoxide forms an addition polymer with DNA, RNA, and proteins in vivo, resulting in damage to the double-stranded structure of DNA and thus leading to cellular carcinogenesis.^37,38 Previous studies have found that DNA damage resulted in the activation p53. DNA damage regulates p53 partially through signaling pathways involving sensor kinases including checkpoint kinase-1/2 and the Wee1 kinase. This results in the transcriptional regulation of a variety of genes involved in apoptosis and cell cycle control,³⁹ thereby modulating cancer progression. Inversely, p53 signaling is involved in cancer progression by regulating DNA damage and biological processes such as cell cycle progression, apoptosis, metabolism, senescence, and programmed necrosis.^40,41 As a tumor suppressor, p53 is prone to inactivation by destabilizing mutations, resulting in tumor development.^41,42 Research has revealed that the Wnt signaling pathway participates in both tissue homeostasis and development through modulation of endogenous stem cells.⁴³ Abnormalities in Wnt signaling have been shown to play critical roles in the occurrence and/or development of many tumors.⁴³ Collectively, our findings suggested that high EC content in musalais could play a carcinogenic role or affect metabolic and biological processes in rats by regulating these pathways.

Venn diagram analysis of the KEGG pathways and enriched DEGs indicated that Mgst1, Gstp1, Gsta1, Gsta2, Gsta5, and Adh1 were concurrently enriched in drug metabolism – cytochrome P450, chemical carcinogenesis, and metabolism of xenobiotics by cytochrome P450 pathways in both High and Mid groups. Among these, Ccnd1 was concurrently enriched in Wnt and p53 signaling pathways. Both Mgst1⁴⁴ and Gstp1⁴⁵ have been demonstrated to exert important functions in oxidative stress and tumorigenesis. Gsta1 is a Phase II detoxification enzyme involved in the biological progression of colorectal, lung, prostate, and bladder cancer,⁴⁶ and the relationship between cancer and Gsta2,⁴⁷ Gsta5,⁴⁸ and Adh1⁴⁸ has been reported. Ccnd1 belongs to the highly conserved cyclin family, whose members are characterized by abundant expression during the cell cycle.⁴⁹ In turn, Ccnd1 has been associated with Wnt⁵⁰ and p53 signaling⁵¹ in affecting the biological progression of tumors. These genes are helpful in understanding the biological response of rats to high EC content in musalais.

Overall, our work predicted that high EC content in musalais possibly induces cancer or affects metabolic and biological processes by regulating cytochrome P450 metabolism, chemical carcinogenesis, metabolism of xenobiotics by cytochrome P450, Wnt signaling pathway, and p53 signaling pathway. This may be achieved by targeting Gstp1, Gsta5, Gsta1, Adh1, Gsta2, and Ccnd1, as revealed by transcriptomic sequencing. The findings provide a useful reference for the safety evaluation of musalais. The limitation of the current study is that only two RNA-seq libraries were constructed in each group, which might increase the random error of statistical analysis. Future research will aim to address this issue by implementing methods to minimize random statistical error. Prospective studies will also be designed to evaluate the effect of high EC content in musalais on the metabolic and biological processes of rats and investigate whether the underlying mechanism is mediated by the pathways and genes identified in this research.