Transcriptome sequencing combined with bioinformatics predicts potential genes and pathways associated with bupivacaine-induced apoptosis.

This study aimed to explore the potential genes and pathways associated with bupivacaine-induced apoptosis. Human neuroblastoma cell line SH-SY5Y was used in this study. The effect of bupivacaine on cell viability of SH-SY5Y was detected by Cell Counting Kit-8. Transcriptome sequencing was performed for SH-SY5Y cells that were treated and untreated with bupivacaine based on the HiSeq 4000 sequencing platform. The sequencing results were analyzed using bioinformatics methods, including differentially expressed genes (DEGs) identification, functional enrichment analysis, protein-protein interaction (PPI) network analysis and module analysis. The cell viability of SH-SY5Y cells decreased significantly after bupivacaine treatment (p < .01). Based on the HiSeq 4000 sequencing platform, we obtained a global overview of the transcriptome of SH-SY5Y treated with/without bupivacaine. Bioinformatics analysis identified 335 up-regulated and 294 down-regulated DEGs in bupivacaine group. They were significantly enriched in cell cycle-associated functions and pathways and cAMP signaling pathway. In the PPI network, proliferating cell nuclear antigen (PCNA), v-Akt murine thymoma viral oncogene homolog 3 (AKT3), cyclin-dependent kinase inhibitor 1A (CDKN1A) and cell division cycle 6 (CDC6) had high topology scores. Module analysis obtained two sub-network modules (cluster 1 and cluster 2). PCNA, CDC6, CDKN1A and AKT3 may play important roles in bupivacaine-induced apoptosis. Additionally, bupivacaine may also induce apoptosis via pathways of cell cycle and cAMP signaling pathway.

Introduction

Local anesthetics are useful compounds, exerting effect by blocking nerve impulse propagation, but increasing studies suggest that local anesthetics are potentially neurotoxic (Brown et al. 1995; Kasaba et al. 2003). More and more studies indicated that neurotoxicity of local anesthetics was associated with apoptosis (Boselli et al. 2003; Johnson et al. 2004; Perez-Castro et al. 2009). For instance, dibucaine, tetracaine and procaine have been reported to cause neuroblastoma apoptosis (Tan et al. 2002). Arita et al. (2000) found that dibucaine could induce apoptosis in promyelocytic leukemia cells. Additionally, Gold et al. (1998) have suggested that lidocaine can induce death in dorsal root ganglion neuron of rats. Bupivacaine is a local anesthetic as well, which has been showed to induce apoptosis in human histiocytic lymphoma U937 cells and in Schwann cells (Arai et al. 2002). However, the underlying molecular mechanisms have not been fully understood.

It has been reported that bupivacaine is a potent uncoupler of mitochondrial oxidative phosphorylation, which could decrease the ATP production and collapse the mitochondrial membrane potential and then lead to apoptosis (Terada et al. 1990). Park et al. (2005) reported that bupivacaine-induced apoptosis was mediated by the generation of reactive oxygen species. Moreover, a recent study suggested that p38 MAPK pathway is associated with bupicavaine-induced apoptosis (Lu et al. 2011). Although some progresses have been made about the mechanisms of bupivacaine-induced apoptosis, it is still important to further explore the potential biomarkers related to the neurotoxicity of bupivacaine. Transcriptome sequencing is a rapidly developing approach for unbiased transcriptome analysis, which can generate an unprecedented global view of the transcriptome (Mäder et al. 2011; Zhang et al. 2013). Bioinformatics analyses play an integrative role in genomics (Bellgard et al. 2003). In this study, after investigating the effects of bupivacaine on the apoptosis of human neuroblastoma cells SH-SY5Y, we performed transcriptome sequencing for SH-SY5Y cells that were treated and untreated with bupivacaine. Then the results of transcriptome sequencing were analyzed by bioinformatics. We aimed to explore the potential genes and pathways associated with bupivacaine-induced apoptosis through transcriptome sequencing combined with bioinformatics.

Methods and materials

Cell line and cell culture

Human neuroblastoma cell line SH-SY5Y was obtained from cell bank of Chinese Academy of Sciences (Shang Hai, China). Cells were cultured in Dulbecco Modified Eagle Medium medium containing 10% FBS and 1% penicillin and streptomycin (Gibco, USA) in an atmosphere of 5% CO2 at 37°C.

Cell viability assay

SH-SY5Y cells were divided into four groups: without treatment (0 mM bupivacaine); treated with 0.5 mM bupivacaine; treated with 1.0 mM bupivacaine; and treated with 1.5 mM bupivacaine. Bupivacaine were obtained from Shanghai Harvest Pharmaceutical Co., LTD (Shanghai, China). The number of viable cells in four groups after treatment was evaluated using the Cell Counting Kit-8 (Dojindo, Kumamoto, Japan) following the manufacturer’s instructions. Briefly, before bupivacaine treatment, 100 μL cell suspension were inoculated into 96-well plates (1 × 104 cells per well) and cultured in a 37°C incubator overnight. Then cells were treated with different concentrations of bupivacaine for 24 h. After that, 10 μL CCK8 was added to each well and incubated for 2 h in incubator. The optical absorbance at wavelength 450 nm was measured for each well using the plate reader (BioTEK, USA).

RNA isolation, library preparation and sequencing

Total 2 mL cell suspension were inoculated into six-well plates (1.5 × 105 cells per well) in triplicate and cultured in a 37°C incubator overnight. Cells were divided into two groups (two repeat for each group): without treatment (0 mM bupivacaine) and treated with 1.5 mM bupivacaine. After 24 h treatment with bupivacaine, cells were washed twice with phosphatic buffer solution and collected into centrifuge tube. Then 1 mL Trizol was added at room temperature and the cell lysates were stored at –80°C after 15 min.

Total RNA from cell lysates was isolated using Phol/chloroform. The quality and integrity of the total RNA were evaluated using the Agilent 2100 Bioanalyzer and Agilent RNA 6000 Nano Kit (Agilent Technologies, USA), respectively. Then mRNA was purified using oligo (dT) magnetic beads, and the high-quality mRNA samples from each group were pooled for preparation of cDNA library using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina®. The mRNA was fragmented into small pieces using an additional fragmentation buffer, then the first-strand cDNA was synthesized using random hexamer-primers. The cDNA was further converted into double-stranded cDNA which was subsequently purified with AMPure XP beads. The purified double-stranded cDNA was performed end repair, dA tailing and adaptor ligation. Size selection was performed using AMPure XP beads, generating cDNA libraries. The libraries were enriched using PCR. The paired-end libraries were then sequenced on the Illumina HiSeq 4000 platform, generating paired-end reads of 150 bp.

The data are available at National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database under the accession number SRP078920.

Primary processing of raw reads and mapping of RNA-Seq reads

The obtained raw reads were quality-filtered by removing the adaptor sequences and any ambiguous or low-quality reads using Cutadapt (version 1.2.1) with default parameters (Martin 2011). After the reference genome index was built using Bowtie2 2.1.0, the obtained clean reads were aligned to the reference genome (UCSC hg19) using Tophat2 2.1.0 (Trapnell et al. 2009).

Gene expression analysis

After read mapping, the count values of reads that mapped to each gene were obtained using HTSeq 0.6.1p2 (Anders et al. 2015). The read count values were considered as the original expression values of genes. Finally, we obtained the read count matrixes of four samples (two control groups and two bupivacaine groups).

Differentially expressed genes analysis

The read count data were first normalized using the trimmed mean of M values (TMM) of edgeR (Dai et al. 2014) package. Then the normalized data were transferred into gene expression matrix with voom (Law et al. 2014) method in limma (Ritchie et al. 2015; Smyth 2005) package. The voom could transfer read count data into log2-counts per million (logCPM) which could calculate the differential expression levels between control and bupivacaine groups using limma package. The significant p-values were calculated using moderated t-statistics and adjusted with Benjamini–Hochberg method. The adjusted p-value < .05 and |log2 (fold change)| ≥ 0.58 were considered as the cutoff values. The heatmap of differentially expressed genes (DEGs) was obtained using ComplexHeatmap (http://www.bioconductor.org/packages/devel/bioc/html/ComplexHeatmap.html) (Gu et al. 2016) package in R.

Functional enrichment analyses

The identified DEGs were performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using Database for Annotation, Visualization, and Integrated Discovery (DAVID) online tool (http://david.abcc.ncifcrf.gov/) (Da Wei Huang and Lempicki 2008). The GO mainly includes biological process (BP), molecular function (MF) and cellular component (CC) categories. The BP, MF and CC terms and KEGG pathways that enriched by at least 2 DEGs and with p-value < .05 were considered as significant enrichment results.

Protein–protein interaction network analysis

The Protein–protein interaction (PPI) pairs of proteins that encoded by DEGs were retrieved using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, http://string-db.org/) (Szklarczyk et al. 2014). The prediction methods of STRING database come from the databases of Gene Fusion, Neighborhood, Co-expression Experiments, Co-occurrence, Databases and Textmining. Based on the predicted PPI pairs, the PPI network was visualized using Cytoscape (http://www.cytoscape.org/) (Shannon et al. 2003).

Three calculation methods of network topological features were used to analyze the importance of nodes in the PPI network. The three methods included degree centrality, betweenness centrality and closeness centrality, which were used to calculate the topology scores of nodes in the PPI network using the CytoNCA (Tang et al. 2015) plugin in Cytoscape. The nodes with higher scores may play important roles in the PPI network, which were hub nodes.

PPI network module analysis

In the PPI network, the proteins with similar functions tend to cluster. Therefore, analysis of the functional clustering module in the PPI network may help us understand the unknown functions of proteins. In this study, MCODE (Cline et al. 2007) plugin in Cytoscape was used to analyze the sub-network module in the PPI network with module score >4.5. In MCODE, each module has a score. The higher the score, the closer the interaction among nodes in module. The DEGs involved in the module were then performed GO and KEGG pathway enrichment analyses.

Results

Bupivacaine inhibited cell viability

The effect of bupivacaine on SH-SY5Y cell viability was detected by Cell Counting Kit-8. As shown in Figure 1, compared with 0 mM bupivacaine, the cell viability of SH-SY5Y in 0.5 mM bupivacaine group decreased after 24 h treatment (p > .05). With the increase of the concentration of bupivacaine, the cell viability in 1.0 and 1.5 mM bupivacaine groups decreased significantly in comparison with 0 mM bupivacaine group (p < .01).

Figure 1.

Effect of bupivacaine on cell viability of human neuroblastoma cell line SH-SY5Y. 0 mM: control group (not treated with bupivacaine); 1.0 mM: treated with 1.0 mM bupivacaine; and 1.5 mM: treated with 1.5 mM bupivacaine. ** p < .01 compared with control group.

Because 1.5 mM bupivacaine had the best inhibiting effect on SH-SY5Y cell viability, we used 1.5 mM bupivacaine in the treatment groups in the subsequent experiments.

Analysis of RNA-Seq data

Based on the HiSeq 4000 sequencing platform, we obtained a global overview of the transcriptome of SH-SY5Y treated with/without bupivacaine. Across all 4 samples, between 27,443,982 and 28,274,506 raw reads with raw bases ranged from 4,116,597,300 to 4,241,175,900 were generated. After stringent filtration, 27,659,392, 26,923,806, 27,769,440 and 27,561,998 clean reads were obtained, respectively, representing 97.86%, 98.11%, 98.21% and 98.11% of the raw reads (Table 1). After mapping reads to the reference genome, we obtained the overall read mapping rates ranged from 91.20% to 92.50%, and the concordant pair alignment rates ranged from 84.3% to 86.7%.

Table 1.

The results of data filtering.

Sample	Control_1	Control_2	Bupivacaine _1	Bupivacaine _2
Raw reads number	28,263,528	27,443,982	28,274,506	28,092,912
Raw bases number	4,239,529,200	4,116,597,300	4,241,175,900	4,213,936,800
Clean reads number	27,659,392	26,923,806	27,769,440	27,561,998
Clean reads rate (%)	97.86	98.11	98.21	98.11
Clean bases number	4,148,908,800	4,038,570,900	4,165,416,000	4,134,299,700
Low-quality reads number	233,254	240,916	187,950	230,560
Low-quality reads rate (%)	0.82	0.88	0.67	0.82
Adapter polluted reads number	370,882	279,260	317,116	300,354
Adapter polluted reads rate (%)	1.31	1.02	1.12	1.07
Raw Q30 bases rate (%)	94.24	94.13	94.86	94.2
Clean Q30 bases rate (%)	94.73	94.73	95.37	94.77
Ns reads number	0	0	0	0
Ns reads rate (%)	0	0	0	0

Differential expression analysis

A total of 629 DEGs including 335 up-regulated and 294 down-regulated DEGs were identified in bupivacaine groups with adjusted p-value < .05 and |log2(fold change)| ≥ 0.58. The list of DEGs was shown in Supplementary Table 1.

Functional enrichment analysis showed that the up-regulated DEGs were significantly enriched in BP terms of GO:0045668∼negative regulation of osteoblast differentiation, GO:0035136∼forelimb morphogenesis, and GO:0048167∼regulation of synaptic plasticity, CC terms of GO:0005905∼coated pit, GO:0031594∼neuromuscular junction, and GO:0043025∼neuronal cell body, and MF terms of GO:0042803∼protein homodimerization activity, and GO:0008022∼protein C-terminus binding (Figure 2(a)). Additionally, they were significantly enriched in 7 pathways, such as hsa05215:Prostate cancer, hsa04261:Adrenergic signaling in cardiomyocytes and hsa04024:cAMP signaling pathway (Figure 2(b)).

Figure 2.

The CC, BP and MF terms, and KEGG pathways enriched by (a and b) up-regulated and (c and d) down-regulated DEGs. The x-axis represents the number of enriched gene.

The down-regulated DEGs were significantly enriched in BP terms associated with cell cycle, CC terms of GO:0005654∼nucleoplasm and GO:0042555∼MCM complex, and MF terms of GO:0004003∼ATP-dependent DNA helicase activity, GO:0051082∼unfolded protein binding (Figure 2(c)). Moreover, five pathways were significantly enriched by the down-regulated DEGs, such as hsa05130:Pathogenic Escherichia coli infection, hsa00240:Pyrimidine metabolism, and hsa03030:DNA replication (Figure 2(d)).

PPI network analysis

A PPI network with 220 nodes and 465 PPI pairs was collected (Figure 3). The top 15 nodes with high topology scores that were calculated by three centrality methods were shown in Table 2, which could be considered as hub nodes. Proliferating cell nuclear antigen (PCNA) had the highest score in three centrality methods. In addition, valosin containing protein, v-Akt murine thymoma viral oncogene homolog 3 (AKT3), cyclin-dependent kinase inhibitor 1A (CDKN1A) and cell division cycle 6 (CDC6) also had high topology scores.

Figure 3.

The PPI network constructed by the DEGs. White node stands for up-regulated gene and gray node stands for down-regulated gene.

Table 2.

The top 15 nodes with high topology scores that were calculated by three centrality methods.

Symbol	Degree	Symbol	Betweenness	Symbol	Closeness
PCNA	35	PCNA	9877.4795	PCNA	0.0408125
TUBB4B	20	VCP	6175.589	HSPA5	0.0406685
TUBB	20	AKT3	6111.1416	AKT3	0.0406082
CDKN1A	19	HSPA5	4892.837	TUBB	0.040533
HSPA5	19	PCSK1	4049.953	TUBB4B	0.040518
VCP	18	TUBB	3755.4736	VCP	0.0404956
TUBA1B	17	KCNMA1	2987.0981	TUBA1B	0.0403538
AKT3	17	CDKN1A	2955.5698	CDKN1A	0.0403389
MCM3	15	TUBB4B	2805.5461	FUS	0.0403315
CDC6	15	FUS	2766.349	TUBA1A	0.0403092
TUBA1A	15	IGF2	2758.5251	HDAC9	0.040287
MCM7	15	HDAC9	2656.2869	HSPA6	0.0402278
HDAC9	15	NPY	2652.3323	SMARCA4	0.0401909
ASF1B	13	CRHR1	2415.4155	ASF1B	0.0401761
PLK2	13	SYT1	2405.1887	TBPL2	0.0401761

Based on the threshold of module score >4.5, we obtained two sub-network modules (cluster 1 and cluster 2) (Figure 4). Cluster 1 consisted of 13 nodes and 37 interaction pairs, such as CDKN1A. Genes in this cluster were significantly enriched in GO terms such as GO:0006260∼DNA replication, GO:0000082∼G1/S transition of mitotic cell cycle, GO:0042555∼MCM complex and GO:0005524∼ATP binding, as well as pathways such as hsa04110:Cell cycle and hsa03030:DNA replication (Figure 5(a)). Cluster 2 included 15 nodes and 33 interaction pairs, such as FUS RNA binding protein (FUS). Genes in cluster 2 were significantly enriched in GO terms of GO:0000398∼mRNA splicing, via spliceosome, GO:0043209∼myelin sheath, and GO:0042288∼MHC class I protein binding, and pathways of hsa03040:Spliceosome and hsa04141:Protein processing in endoplasmic reticulum (Figure 5(b)).

Figure 4.

The sub-network module obtained from the PPI network. White node stands for up-regulated gene and gray node stands for down-regulated gene.

Figure 5.

The CC, BP and MF terms, and KEGG pathways enriched by the DEGs of (a) cluster 1 and (b) cluster 2.

Discussion

Bupivacaine, as a long-acting local anesthetic, has been used for controlling post-operative and chronic pain (Kohase et al. 2002). However, results from our present study indicated that bupivacaine significantly induce apoptosis of SH-SY5Y cells. Bioinformatics analysis of the results of transcriptome sequencing identified 629 DEGs between bupivacaine and control groups. These DEGs, such as PCNA, CDC6 and CDKN1A, were significantly enriched in cell cycle-associated pathways. Additionally, PCNA, CDC6, CDKN1A and AKT3 had higher topology scores in the PPI network.

Cell cycle is a vital process which occurs in cell, leading to the division of cell and duplication of its DNA to produce two daughter cells. It has been reported that cell cycle arrest may lead to apoptosis (Santos et al. 2011). In this study, the down-regulated DEGs, such as PCNA and CDC6, were significantly enriched in cell cycle-associated functions and pathways, suggesting the role of cell cycle in bupivacaine-induced apoptosis. PCNA is a DNA clamp, acting as a processivity factor for DNA polymerase δ in eukaryotic cells, besides it is essential for replication (Moldovan et al. 2007). Essers et al. (2005) have reported that PCNA is important for both DNA synthesis and DNA repair. Importantly, PCNA overexpression can inhibit cell apoptosis (Zhao et al. 2014). For CDC6, it plays an important role in initiating DNA replication during the S phase of the cell cycle (Dalton and Whitbread 1995). Silencing of CDC6 may inhibit the DNA replication and arrest cells in G1/S phase resulting in inability to enter S phase (Niimi et al. 2012). Therefore, in this study, we speculated that bupivacaine may promote apoptosis by down-regulating PCNA and CDC6.

In addition, DEGs in cluster 1, such as CDKN1A, were also significantly enriched in cell cycle-related functions and pathways. CDKN1A, also known as p21, is a potent cyclin-dependent kinase (CDK) inhibitor, which can bind to CDK2, CDK1, and CDK4/6 complexes and regulate cell cycle progression at G1 and S phase (Gartel and Radhakrishnan 2005). It has been reported to regulate cell cycle progression, terminal differentiation and apoptosis (Wagner et al. 2001). Additionally, CDKN1A expression is tightly controlled by the tumor suppressor protein p53 which can initiate apoptosis if DNA damage proves to be irreparable. Apoptosis induced by p53 is considered as a central mechanism of tumor suppression (Speidel 2010). Thus, CDKN1A may play a role in bupivacaine-induced apoptosis.

Interestingly, PCNA, CDC6 and CDKN1A were hubs in the PPI network. In addition to them, AKT3 had higher topology scores as well. AKT3 is a member of the AKT subfamily of serine/threonine protein kinases which are associated with a variety of BPs including cell proliferation, apoptosis, differentiation and tumorigenesis (Xu et al. 2012). It was found to be enriched in cAMP signaling pathway in this study. The cAMP signaling pathway is an important intra cellular messenger delivery system, which is implicated with the abnormal proliferation and apoptosis of malignant lymphoid cells. In malignant lymphoid cells, manipulating this pathway could induce cell cycle arrest and apoptosis (Wang 2014). We speculated that bupivacaine may induce apoptosis by targeting AKT3 via cAMP signaling pathway.

Although several genes and pathways were identified in the present study, no experiment was performed to validate the expression level of these DEGs, including PCNA, CDC6, CDKN1A and AKT3 in bupivacaine-treated cells, which was really a limitation of our study. In the future, we will conduct experimental validations with larger sample size to validate our results.

In conclusion, PCNA, CDC6, CDKN1A and AKT3 were differentially expressed in bupivacaine-treated group, which may play important roles in bupivacaine-induced apoptosis. Additionally, bupivacaine may also induce apoptosis via pathways of cell cycle and cAMP signaling pathway. Findings of this study may provide a novel insight into the molecular mechanisms underlying the neurotoxicity of bupivacaine and may provide important implications for future research.

Geolocation information

No.29, Longtan Road, Taishan District, Tai’an, Shandong Province, China