Transcriptomics reveals the effects of NTRK1 on endoplasmic reticulum stress response-associated genes in human neuronal cell lines

Cloning and plasmid construction

Overexpressed NTRK1 lentiviruses were purchased from GenePharma Co., Ltd. (Suzhou, China). Transcript information included NM_002529.4, and the lentiviral vector was LV5 (EF-1a/GFP/Puro/Amp).

Cell culture and transfection

The SH-SY5Y cell lines (Procell Life Science & Technology Co., Ltd., Shanghai, China) were cultured at 37 °C with 5% CO₂ in a minimum essential medium (MEM)/F12 (Procell Life Science & Technology Co., Ltd., Wuhan, China) supplemented with 10% fetal bovine serum (FBS; Gibco, New York, NY, USA), 100 μg/mL streptomycin, and 100 U/mL penicillin. SH-SY5Y cells were infected with lentiviruses (multiplicity of infection (MOI) = 35). Infected cells harvested for 48 h were subjected to quantitative reverse transcription-polymerase chain reaction (RT-qPCR), Western blot, and Cell Counting Kit-8 (CCK-8) assays.

Assessment of gene overexpression

NTRK1 overexpression was evaluated using Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as a control gene. In the present study, cDNA synthesis was performed by standard procedures, and RT-qPCR was carried out on the Bio-Rad S1000 with Hieff® qPCR SYBR® Green Master Mix kit (Low Rox Plus; YEASEN, Shanghai, China). As shown in Table S1, the primer sequences were listed. Using the 2^−ΔΔCT method, the transcript concentration was normalized to the GAPDH mRNA level (Livak & Schmittgen, 2001). By using GraphPad Prism, we conducted comparisons via the paired Student’s t-test (GraphPad Software Inc., San Diego, CA, USA).

Western blotting

In ice-cold lysis buffer (1× phosphate-buffered saline (PBS), 0.1% sodium dodecyl-sulfate (SDS), 0.5% NP-40, and 0.5% sodium deoxycholate) supplemented with a protease inhibitor cocktail (Roche, Basel, Switzerland), SH-SY5Y cells were lysed and incubated on ice for 30 min. Afterward, we boiled samples with 1× SDS sample buffer for 10 min, separated them by 10% SDS-polyacrylamide gel electrophoresis (SDS-PAGE), and transferred them to membranes. With TBST buffer (20 mM Tris-buffered saline and 0.1% Tween-20) containing 5% non-fat milk powder, membranes were incubated for 1 h at room temperature with FLAG (1:1,000; Cat. No. 2368S; Cell Signaling Technology Inc., Danvers, MA, USA), GAPDH (1:5,000; Cat. No. 60004-1-Ig; Proteintech, Rosemont, IL, USA), GRP78 (1:3,000; Cat. No. 11587-1-AP; Proteintech, Rosemont, IL, USA), p-IRE1 (1:500; Cat. No. AF7150; Affinity, West Bridgford, UK), XBP1s (1:2,000; Cat. No. 24868-1-AP; Proteintech, Rosemont, IL, USA), and ATF6 (1:1,000; Cat. No. DF6009; Affinity, West Bridgford, UK) primary antibodies, followed by incubation with horseradish peroxidase (HRP)-conjugated secondary antibody. Enhanced chemiluminescence reagent (ECL; Cat. No. 170506; Bio-Rad Laboratories Inc., Hercules, CA, USA) was used to visualize the immunoblots.

RNA extraction and sequencing

In SH-SY5Y cells, total RNA was extracted with TRIzol reagent (Cat. No. 15596026; Invitrogen, Carlsbad, CA, USA), as described by Chomczynski & Sacchi (1987). After RNA extraction, DNaseI was used to digest DNA. With the Nanodrop™ One Cspectrophotometer system (Thermo Fisher Scientific Inc., Waltham, MA, USA), RNA quality was assessed using A260/A280. Electrophoresis of 1.5% agarose gel was used to confirm RNA integrity. Final quantification of qualified RNAs was performed with the Qubit™ RNA Broad Range Assay kit (Cat. No. Q10210; Life Technologies Corp., Carlsbad, CA, USA) using the Qubit 3.0 Fluorometer.

In addition, KCTM Stranded mRNA Library Prep kit for Illumina (Cat. No. DR08402; Wuhan Seqhealth Co., Ltd., Wuhan, China) was used to prepare stranded RNA sequencing libraries. On the Novaseq 6000 sequencer (Illumina; Chicago, IL, USA) with PE150 model, 200–500 bps PCR products were enriched, quantified, and finally sequenced.

RNA-seq data processing

First, raw reads with more than 2-N bases were discarded. Then, FASTX-Toolkit (Ver. 0.0.13) was used to trim adaptors and low-quality bases from raw sequencing read. The short reads <16-nt were also dropped. After that, An alignment of clean reads to the GRCh38 genome was performed by HISAT2 (Kim, Langmead & Salzberg, 2015) with four mismatches allowed. Reads that were uniquely mapped were counted and fragments per kilobase of transcript per million fragments mapped (FPKM) were calculated (Trapnell et al., 2010).

CCK-8 assay

Cell proliferation was measured using a CCK-8 assay kit (MCE Co., Ltd., Beijing, China). SH-SY5Y cells were seeded into 96-well plates in triplicate at a density of 2 × 10⁴ cells/well and cultured for 72 h. At 37 °C for 2 h, the cells were then incubated in 10% CCK-8 solution diluted in the MEM/F12. An enzyme-labeling instrument (Bio-Rad Laboratories Inc., Hercules, CA, USA) was used to measure the absorbance (optical density (OD) value) at 450 nm wavelength.

Analysis of DEGs

DEGs were screened with DESeq2 from R Bioconductor (Love, Huber & Anders, 2014). To identify DEGs, we used a P value of 0.05 and a fold-change (FC) of 1.5 or 0.67 as cut-off criteria.

Functional enrichment analysis

With KOBAS 2.0, GO and KEGG pathway enrichment analyses were conducted to sort DEGs into functional categories (Xie et al., 2011). To determine the enrichment of each term, the hypergeometric test and Benjamini-Hochberg false discovery rate (FDR) controlling procedure were applied.

PPI network construction, modules, and identification of top 10 genes

Search tool for the retrieval of interacting genes (STRING) database (http://stringdb.org/) (Szklarczyk et al., 2019) was employed to construct PPI networks by mapping the DEGs into PPI data. The symbols of the DEGs were imported into the database, and high-resolution bitmaps were generated. The bitmap only included interactors with a combined confidence score of 0.4 or higher. The hotspot module was obtained in large PPI networks by the Molecular Complex Detection (MCODE) which is a Cytoscape plugin. In this study, MCODE parameters were as follows: the degree of cut-off = 2; cluster finding, haircut; node score cut-off = 0.2; k-core = 2; and the maximum depth = 100 (Bader & Hogue, 2003). The modules with established score >9 and the number of nodes >9 were screened out. The GO enrichment analysis was used to analyze genes in the modules by DAVID. CytoHubba, a plugin for Cytoscape 3.9.1, was used to calculate the degree of connectivity and identify the top 10 genes in the PPI networks (Shannon et al., 2003). The PPI networks, modules, as well as the top 10 genes were visualized based on their node degree using Cytoscape software.

RT-qPCR for validation of the hub genes

An RT-qPCR analysis was performed and normalized with GAPDH to clarify the validity of the top 10 hub genes. RT-qPCR was performed using the same cell lines as RNA-seq. For the amplification program, denaturing was carried out for 30 s at 95 °C, followed by 40 cycles of denaturing at 95 °C for 10 s and 60 °C for 30 s, then annealing and extension at 95 °C for 15 s, 60 °C for 60 s, and 95 °C for 15 s. For each sample, the experiment was conducted in triplicate. Table S1 illustrates the list of the primer sequences used for RT-qPCR analysis.

Functional annotation of the hub genes

The ClueGO (ver. 2.5.8) and CluePedia (ver. 1.5.8) Cytoscape plugins were performed to decipher the enrichment analysis of the hub genes. GO-biological process (GO-BP) and pathways were considered statistically significant at P < 0.05.

Statistical analysis

According to RNASeqPower, this experimental design has a statistical power of 0.96. The cell biology and RT-qPCR data between the two groups were compared using an unpaired two-tailed t-test. All values were presented as the mean ± standard deviation (SD). Except for the RNA-seq experiment, each experiment was performed at least three times. There was a statistically significant difference at P < 0.05.

NTRK1 overexpression promoted the proliferation of SH-SY5Y cells

A lentiviral vector expressing the NTRK1 gene (NTRK1-OE) or an empty lentiviral vector (not expressing the NTRK1 gene) (NC) was used to transfect SH-SY5Y cells. The results of the RNAseq, RT-qPCR, and Western blotting showed that NTRK1 was significantly overexpressed in SH-SY5Y cells (Figs. 1A and 1B). Moreover, the overexpression of NTRK1 significantly promoted SH-SY5Y cell proliferation (Fig. 1C). The results validated that NTRK1-OE was successfully constructed and NTRK1 favored the survival of the neurons.

NTRK1 overexpression altered the gene expression profiles in SH-SY5Y cells

RNA-seq was performed to examine the NTRK1-mediated transcriptional regulation in SH-SY5Y cells. In total, six RNA-seq libraries were prepared from NTRK1-OE and NC SH-SY5Y cells, with three biological replicates per group (NTRK1_1, NTRK1_2, and NTRK1_3; NC_1, NC_2, and NC_3). After removing low-quality reads and sequence adaptors, a clean pair-end read was obtained for 69.1 million samples on average. Approximately 62.8 million read pairs were uniquely mapped per sample to the human genome (Table S2).

Furthermore, gene expression was calculated using these uniquely mapped reads. A proprietary pipeline was used to calculate FPKM, which represents gene expression levels. There were 22,154 genes expressed (FPKM > 0), and 10,364 genes expressed at FPKM > 1 in at least one sample (Table S3). A correlation matrix using FPKM values for all 22,154 genes was calculated using Pearson’s correlation value (more than 0.99) between NTRK1-OE and the control, which indicated the similar expression of most genes. Moreover, NTRK1-OE and NC samples were clearly differentiated by principal component analysis (PCA) based on FPKM values for all detected genes (Fig. 1D), and revealed that SH-SY5Y cells overexpressing NTRK1 showed altered gene expression profiles.

For further analysis of NTKR1’s transcriptional regulation, DESeq2 was performed to identify the genes that are differentially expressed between NTRK1-OE and NC cells, with a cut-off FC ≥ 1.5 or ≤0.67 and an FDR of 5%. The gene expression level in SH-SY5Y cells was extensively regulated by NTRK1 overexpression (Fig. 1E), as shown by 193 upregulated and 226 downregulated genes (Table S4). Additionally, the NTRK1-OE and NC samples were clearly distinct in the hierarchical clustering analysis of normalized FPKM values of DEGs, and the three replicate datasets also showed a high level of consistency (Fig. 1F). These results indicated that NTRK1 overexpression significantly changed the transcript expression level of a series of genes in neurons.

NTRK1 overexpression affected the biological process in SH-SY5Y cells

The GO and KEGG pathway enrichment analyses of all 419 DEGs were carried out to determine their potential roles. A significant increase in the expression levels of 193 genes was observed. According to the GO enrichment analysis, these genes were mostly enriched in response to endoplasmic reticulum (ER) stress, protein folding in ER, IRE1-mediated unfolded protein response (UPR), collagen fibril organization, ATF6-mediated UPR, cellular protein metabolic process, response to unfolded protein, ER unfolded protein response, folic acid metabolic process, and extracellular matrix (ECM) organization (Fig. 2A). The KEGG pathway analysis further confirmed that the upregulated genes were highly enriched in protein processing in ER (Fig. 2C). Moreover, upregulated genes were also highly enriched in a range of biochemical processes and p53 signaling pathway. A total of 226 genes downregulated by NTRK1 were mainly enriched in a series of cellular parts and cellular processes, and pathways associated with cell proliferation and migration, such as Dopaminergic synapses, ECM-receptor interactions, transforming growth factor-β (TGF-β) signaling pathway, and Wnt signaling pathway (Figs. 2B and 2D). Overall, these data suggest that a range of biological processes are affected by NTRK1 overexpression in human neuronal cell lines, in which response to ER stress occupies great importance.

PPI network analysis and hot modules

A PPI network based on the DEGs was constructed on STRING analysis for functional association and sequentially visualized on Cytoscape, and the results are shown in Fig. 3A. Furthermore, there are 294 nodes and 1,062 edges in the network.

Using the Cytoscape plugin “MCODE”, two functional modules (modules 1 and 2) were detected with scores >9 and nodes >9. Module 1 covered 143 edges and 23 nodes, with a score of 13, containing 21 upregulated and two downregulated genes (Fig. 3B). Module 1 functional enrichment analysis revealed that ER stress-related functions dominated the most significant enrichment results in GO-BP (Table 1). Having 42 edges and 10 nodes, module 2 scored 9.33, with 10 genes downregulated (Fig. 3C), and these genes were mainly enriched in mitotic sister chromatid segregation (Table 1). The results of hot modules in the PPI network suggest that ER stress-related functions exert vital roles in neurons with NTRK1 overexpression.

Identification of the top 10 genes and validation of hub genes

Using the cytoHubba plugin, the top 10 target genes were identified based on scores, and interactors of the top 10 genes were reconstructed (Fig. 3D). These top 10 genes with detailed information are presented in Table 2.

The expression levels of the top 10 genes were quantified using RT-qPCR to verify the effects of NTRK1-OE on them. Overall, seven hub genes were verified, including five upregulated and two downregulated DEGs (Fig. 4A). The upregulated DEGs included COL1A1 (|log2FC| = 0.63), P4HB (|log2FC| = 0.65), HSPA5 (|log2FC| =1.27), THBS1 (|log2FC| = 0.73), and XBP1 (|log2FC| = 0.90); the downregulated DEGs included CCND1 (|log2FC| = 0.73) and COL3A1 (|log2FC| = 1.31). Hub genes tended to be enriched in terms of ‘response to endoplasmic reticulum stress’, ‘cellular response to transforming growth factor beta stimulus’, and ‘collagen fibril organization’ as the biological process (Fig. 4B). Furthermore, the western blot results showed NTRK1 over-expression increased expression of ER stress markers in neurons, including GRP78, p-IRE, XBP1s, and ATF6 (Fig. S2). Hence, these data further confirm that ER stress response is dominant in neurons with NTRK1 overexpression.