Introduction

Human polymorphonuclear leukocytes (PMNs or neutrophils) are the most abundant white cells in the blood and serve as the first line of host defense. They are terminally differentiated cells with an extremely short lifespan (8–20 hours in circulation and 1–4 days in tissue) and a daily turnover of 1011 cells1. Spontaneous apoptosis of neutrophil plays an important regulatory role in this dramatic turnover2. It is generally agreed that the death of PMN under physiological condition is spontaneous apoptosis, also known as constitutive apoptosis3, which can be mimicked in vitro by culturing the cell in the absence of sufficient amount of survival cytokines4. It was believed that mature PMNs were nearly transcriptionally inert and only a few proteins were selectively expressed5. However, a systems biology approach revealed that the lifespan of PMN could be regulated at the transcriptional level6. Multiple lines of evidence suggest that PMN death is regulated by both intracellular death/survival signaling pathways and a variety of extracellular stimuli3. The general consensus is that many factors can delay constitutive neutrophil apoptosis, but nothing can completely prevent it7.

A growing body of evidence suggests that rather than being “transcriptional noise” and by functioning as primary regulators to control the expression of many target genes, numerous non-coding RNAs (ncRNAs) are involved in almost every aspect of cellular processes including differentiation, transcription, metabolism, and apoptosis8,9. As machine learning technology entered the bioinformatics field, computational models for non-coding RNAs research such as LRSSLMDA, IMCMDA, LRLSLDA and MDHGI10,11,12,13 have advanced greatly in recent years. Recent studies suggest that microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) can regulate the lifespan of the short-lived myeloid cell14,15. For example, by regulating the transcription of pro-apoptotic gene Bcl2l11(Bim), the lncRNA Morrbid serves as a regulator of the lifespan of PMN16. It has also been reported that some lncRNA-related features could improve the diagnosis of some diseases such as non-small cell lung cancer, multiple myeloma, and bladder cancer17,18,19. It is therefore important to understand the regulatory role of lncRNAs in PMN apoptosis.

The timing of gene expression is important for apoptosis and therefore precisely controlled20. We first used different bioinformatics methods including short time-series expression miner (STEM) analysis screened statistically significant time-dependent gene expression profiles21. Then, using a time-series PMN apoptosis microarray datasets, we performed analyses of STEM and the lncRNA-mRNA co-expression network to identify key genes involved in PMN apoptosis. And finally, we experimentally confirmed the expression profiles of two mRNA (NFκB1 and BIRC3) and two lncRNAs (THAP9-AS1 andAL021707.6) during PMN spontaneous apoptosis in a time-dependent manner.

Results

Human spontaneous neutrophil apoptosis pattern

To obtain the spontaneous neutrophil apoptosis pattern, we performed flow cytometry to quantify apoptotic cells by measuring PS externalization on the PMN surface. Over time, PMN apoptosis increased gradually, having a relatively slow raise between 3 and 6 h, and a relatively quick increase after 6 h (Fig. 1).

Figure 1
figure 1

Human spontaneous neutrophil apoptosis. Apoptosis was assessed at the indicated time points (0, 3, 6, 12, 24 h) using Annexin V-FITC staining and flow cytometry. Pooled data are the mean ± SD (n = 4).

Full size image

The baseline of microarray data and lncRNA expression profiles

The public dataset GSE37416 was obtained from GEO and used in this study22. Twenty samples were selected from GSE37416 (GSM918524-GSM918527, GSM918532-GSM918535, GSM918540-GSM918543, GSM918548-GSM918551, and GSM918556-GSM918559), including samples at 0, 3, 6, 12, and 24 h (n = 4 each). The data analysis pipeline is shown in Fig. 2. A total of 4655 lncRNAs probe sets (matched with 4156 lncRNAs) were selected and re-annotated using the Affymetrix HG-U133 Plus 2.0 array from the BioMart database (Supplementary Table 1).

Figure 2
figure 2

Data analysis flow diagram.

Full size image

Transcriptome changes of differentially expressed lncRNA and mRNA

PCA mapping was used to mathematically reduce the dimensionality of the spectrum of gene expression values and obtain information regarding the total amount of dataset variation23. The PCA data (Fig. 3) shows that the samples from 0, 3, 6, 12, and 24 h are clearly segregated from each other as apoptosis progressed. Microarray hybridization results showed that about 3050 mRNA and 220 lncRNA were differentially expressed when they were compared with 0 h as control (P < 0.01, logFC > 1). The differentially expressed mRNAs and Gene Ontology (GO) analyses at each time point are shown in Table 1 and Fig. 4, and the differentially expressed lncRNAs and transcript types and distribution were shown in Fig. 5. The majority of differentially expressed lncRNAs were antisense and long intergenic non-coding RNAs. These differentially expressed mRNAs were enriched mainly in the processes related to cellular response to DNA damage stimulus, apoptosis process, the inflammatory response, and catabolic processes.

Figure 3
figure 3

Summary of microarray data during apoptosis. PCA graph depicting variation and clusters of all genes in the data set at five time points in 2 dimensions.

Full size image
Table 1 Overview of the differential expression analysis of time-series lncRNAs and mRNAs at the gene level.
Full size table
Figure 4
figure 4

Overview of time-series mRNA differential expression (DE) analysis during neutrophil apoptosis. (A) Heatmap of differentially expressed genes for mRNAs. (B) Top GO terms significantly enriched in genes that are differentially expressed compared to 0 h as control.

Full size image
Figure 5
figure 5

Overview of time-series lncRNA differential expression (DE) analysis during neutrophil apoptosis. (A) Heatmaps of differentially expressed genes for lncRNAs. (B) The transcript types and distribution of 220 differentially expressed lncRNAs.

Full size image

Temporal cluster analysis of significantly differentially expressed genes

We normalized the expression data to 0 h (control) and identified temporal gene expression profiles using STEM. Within the 50 model profiles, we identified 12 significant clusters containing a total of 2342 genes (Fig. 6A–L). The profile boxes depict the gene expression patterns over five time points. The profile number on the top left corner of each profile box was assigned by STEM, and the number on the bottom left represents the adjusted P-value. We found genes with continuous downregulation pattern are in profile 9. Their expression increased in the early stage and dropped later, and they were distributed across profiles 45, 37, and 44. Meanwhile, profiles 1, 3, 4, 5, 7, 13, 33, and 8 showed biphasic expression patterns. GO function analysis for profiles shows that there are two interesting clusters comprised of genes related to apoptosis, namely profiles 4 and 45. The lncRNAs in the two profiles are depicted in Fig. 6M,N. Genes in profiles 4 and 45 are summarized in Supplementary Table 2. Additionally, the expression patterns of profiles 4 and 45 appear to be opposite to each other. GO analysis was performed for the genes in these two clusters based on the DAVID database. As shown in Table 2, the high-enrichment GO terms included: inflammatory response, NIK/NF-κB signaling, negative regulation of apoptotic signaling pathway, microtubule polymerization or depolymerization, execution phase of apoptosis, apoptotic mitochondrial changes, etc.

Figure 6
figure 6

STEM clustering of the differentially expressed genes. (A–L) All significant profiles based on the P values of the numbers of genes assigned vs. expected. Profile 45 (0, 2, 1, 0, 0): 249 genes assigned, 76 genes expected, P-value = 1.6 E-57 (significant); Profile 4 (0, −2, −1, 0, −2): 178 genes assigned, 34 genes expected, P-value = 4.5E-69 (significant); (M) expression pattern of 16 lncRNAs assigned in Profile 45; (N) expression pattern of 6 lncRNAs assigned in Profile 4.

Full size image
Table 2 GO results for profiles 4 and 45.
Full size table

PPI network and lncRNA-mRNA co-expression network analyses for genes in profile 45 and 4

The PPI network of differentially expressed genes was constructed by STRING and visualized by Cytoscape (Fig. 7A), and the most significant module was identified using Cytoscape plug-in MCODE (Fig. 7B). A lncRNA-mRNA co-expression network (Fig. 7C) was built to identify interactions between mRNAs and lncRNAs in profiles 4 and 45. The results showed that gene in these profiles highly interacted to each other, and hub genes in this module were mainly belonged to the NFκB signaling pathway such as NFKB1, RELA, BIRC2, and BIRC3. Among the co-expressed lncRNAs, we found two top-degree lncRNAs, ENST00000504520 (THAP9-AS1) and ENST00000609212 (AL021707.6), that had the highest Pearson’s correlation coefficients (>0.95) vs. NFκB1.

Figure 7
figure 7

PPI network construction, module analysis, and lncRNA-mRNA co-expression network analyses. (A) The PPI network of genes in profiles 45 (light blue) and 4 (orange) was constructed using Cytoscape. (B) The hub gene module was obtained from the PPI network having 13 nodes and 71 edges with gene groups attached (diamonds). (C) The co-expression network was constructed with 13 hub mRNAs (light blue dots) and 20 co-expression lncRNAs (green arrowheads) in profiles 45 and 4.

Full size image

Validation of the key spontaneous neutrophil apoptosis-associated mRNAs and lncRNAs

Finally, we wanted to experimentally validate the changes of bioinformatics-identified key mRNAs and lncRNAs during neutrophil spontaneous apoptosis. The qRT-PCR results (Fig. 8A,B) showed the levels of NFκB1 and BIRC3 increased dramatically during the first 3 hours of incubation and gradually decreased thereafter. Of note, these results are consistent with that from the microarray. Similarly, we estimated the changes of the two lncRNAs, THAP9-AS1 and AL021707.6, and found that these lncRNAs exhibited the exactly same pattern as NFκB1 and BIRC3 mRNAs (Fig. 8C,D).

Figure 8
figure 8

Validation of the expression of key genes. Neutrophils were purified and cultured for 0 h, 3 h, 6 h, 12 h and 24 h. qRT-PCR was performed to measure the relative expression of NFκB1 (A) and BIRC3 (B), and lncRNAs including AL021707.6 (C) and THAP9-AS1(D). β-ACTIN was used as an internal control. Experiments were performed three times. The t-test was used for analysis of each data set and comparison between different groups. P < 0.05 was considered as the level of statistical significance. *P < 0.05; **P < 0.01.

Full size image

Discussion

Spontaneous neutrophil apoptosis plays critical roles in neutrophil homeostasis and resolving inflammation. Recently, many efforts have been made to decipher the molecular mechanisms of spontaneous PMN apoptosis. Most of the earlier studies have been only focusing on the roles of protein factors. More recent researches indicate that due to their extended lengths, lncRNAs could regulate protein expression by binding and sequestering their mRNAs. That is why some lncRNAs are also known as miRNA sponges24. It has been suggested that both miRNAs and lncRNAs are involved in the regulation of myeloid cell lifespan14,15. More recent findings suggest that lncRNA Morrbid is capable of affecting the lifespan of PMN by regulating the transcription of its neighboring pro-apoptotic gene, Bcl2l11 (Bim)16. Since gene expression dynamics are characterized by a time-dependent pattern, the expression profiles of genes at a single time point are insufficient to fully characterize the role of lncRNAs in apoptosis. We aimed to identify molecular events governing apoptosis using STEM to assess lncRNA and mRNA expression profiles.

Comparison of the mRNA transcriptional data of PMNs at 0 h and later time points, along with GO analysis, revealed differential expression of apoptosis-related genes. We also used the STEM platform to investigate how gene expression profiles change continuously during apoptosis process. We selected 50 predetermined temporal model profiles and quantified the genes assigned to each profile. Some distinct gene expression patterns were noticed during apoptosis. For example, genes modulating cell survival such as NFκB1 and BIRC3 in profile #45 significantly increased from hour 0 to hour 3 and reduced afterward. Meanwhile, genes involved in apoptosis in profile #4 significantly decreased between hour 0 and hour 3 and increased afterward. These findings are in consistent with that reported by Ward et al. that NFκB plays a central role in promoting PMN survival25. BIRC3, a member of the inhibitor of apoptosis (IAP) family, is at the downstream of NFκB signaling pathway and changed parallel to that of NFκB. BIRC3 is of the major anti-apoptosis genes and are partly responsible for sustained neutrophilia26. Mechanistically, IAP members can repress apoptosis by inhibiting caspase-3, -7, and -927,28.

Through RT-PCR, we confirmed that at the early stage of apoptosis, NFκB1 is significantly upregulated along with its downstream pro-survival genes (especially BIRC3) and then gradually downregulated as apoptosis proceeds. And this regulation pattern may significantly affect the PMN apoptosis process, because the Annexin V-FITC and flow cytometry of PMN (Fig. 1) showed that PMN apoptosis rate decreased dramatically after the expression of NFκB1 and BIRC3 peaked at 3 h, and increased again when NFκB1 and BIRC3 reduced after 6 h. PPI network analysis of genes in profiles #45 and #4 also confirmed that genes in NFκB pathway are the hub of the apoptotic gene cluster. Together, these results indicate that PMN apoptosis may be strictly regulated by the levels of NFκB and IAP members. Moreover, using STEM algorithm and coexpression analysis, we identified two previously unreported lncRNAs that might be involved in the relation of apoptosis. THAP9-AS1 and AL021707.6 showed similar patterns as NFκB1 and IAP members with an extremely high correlation (PCC > 0.95). We hypothesize that these lncRNAs could be involved in regulating spontaneous neutrophil apoptosis by increasing the level of transcription factor NFκB, and its downstream IAP members to repress the caspases in the early stage of PMN spontaneous apoptosis (0–6 h) and when these lncRNAs were down-regulated at 12 h and 24 h the level of NFκB were also decreased sharply and PMN apoptosis occurs.

In summary, we presented a dynamic picture of the mRNAs and lncRNAs changes during neutrophil spontaneous apoptosis. Our study also has several limitations. For example, our findings are largely based on bioinformatics analysis of the publically available databases and therefore more functional experiments are needed to verify these results. The molecular mechanisms about how the NFkB-related genes affect neutrophil spontaneous apoptosis should also be further elucidated experimentally.

Conclusions

We systematically analyzed mRNA and lncRNA expression changes in the spontaneous neutrophil apoptosis model. The results showed that the expression of both mRNAs and lncRNAs are stage-specific. NFκB1, BIRC3, and two co-expressed lncRNAs are inversely correlate with apoptosis. These findings suggest that these factor could be essential in the regulation of spontaneous neutrophil apoptosis. Finally, our finding not only bettered our understanding in PMN apoptosis but provided potential key regulatory molecule for PMN apotosis and therapeutic targets for over-reactive inflammatory response caused by the abnormality in PMN apoptosis.

Methods

Peripheral blood neutrophil isolation

Heparinized venous blood was obtained from three healthy individuals with the protocol approved by the Institutional Review Board for Human Subjects at Southwest Medical University. Written informed consent was obtained from every subject. All experiments were performed in accordance with relevant guidelines and regulations. Neutrophils were isolated using 3.0% Dextran T-500 (Amersham Biosciences Corp, Piscataway, NJ, USA), followed by density gradient centrifuge separation29. Neutrophils were suspended in phosphate-buffered saline, counted, and diluted to 1 × 107/ml. Trypan blue, Wright-Giemsa staining, and microscopic analysis were performed to estimate neutrophil purity, and the suspensions were routinely >98% neutrophils. PMNs were diluted to 1 × 106/ml in RPMI-1640 medium (10% fetal bovine serum; Gibco, Gaithersburg, MD, USA), transferred into T-25 flask and subsequently incubated at 37 °C with 5% CO2 for 12 h or 24 h.

Detection of spontaneous neutrophil apoptosis using Annexin V in vitro

Spontaneous neutrophil apoptosis was assessed by measuring phosphatidylserine (PS) externalization on the neutrophil surface using Annexin V30. In brief, PMNs (5 × 106/ml) were left untreated in RPMI-1640 medium. Aliquots of PMNs were removed at the indicated time points. An annexin V-fluorescein isothiocyanate (FITC)/propidium iodide (PI) apoptosis detection kit (KeyGen BioTech, Jiangsu, China) was used according to the manufacturer’s instructions.

The baseline of microarray data

The Gene Expression Omnibus (GEO) database was searched to identify available datasets. After careful validation, the dataset GSE37416 containing neutrophil cell culture sampled at 0, 3, 6, 12, and 24 h were selected for this study22. The downloaded raw data were normalized using the log scale robust multi-array analysis with default settings31. All sample datasets were hybridized with the HG-U133 plus 2.0 Array (Affymetrix, Santa Clara, CA, USA), including 54,675 probe sets that are widely used in biological research (http://www.affymetrix.com/analysis/index.affx). Principal component analysis (PCA) was performed to visualize data variance32.

lncRNA annotation pipeline

To obtain lncRNA expression data from the HG-U133 Plus 2.0 Array, the annotation was downloaded from the BioMart database (http://asia.ensembl.org/biomart/martview/). Only the probes annotated as lncRNAs were selected; transcript ID, chromosome location, strand, biologic types, and other annotation information were also downloaded33.

Identification of differentially expressed lncRNA and mRNAs

The expression data of raw CEL files were normalized, log2 transformed and the background was adjusted utilizing a Bioconductor package Robust MultiArray Average (RMA)34, R 3.5.1 software. Next, a set of probe ID-centric gene expression values was retrieved for downstream analyses. The normalized data were then analyzed with linear models for microarray data (LIMMA), a modified t-test incorporating the Benjamini–Hochberg multiple hypotheses correction technique35. Probe sets without corresponding gene symbols or genes with more than one probe set were removed or collapsed to max probe through GSEA 3.036. LogFC (fold change) >1 and adjusted P-values < 0.01 were considered statistically significant. A Gene Ontology (GO) analysis was performed using the Database for Annotation Visualization and Integrated Discovery (DAVID)37.

STEM and GO analyses

We used the STEM21 clustering algorithm to identify temporal gene expression profiles during spontaneous neutrophil apoptosis with a maximum number of model profiles set as 50 and a maximum unit change in model profiles between time points set at 2. Gene expression values were transformed to log ratios relative to the expression value at 0 h. Then, each gene was assigned to the filtering criteria of the model profiles, and the correlation coefficient was determined. Standard hypothesis was performed using the true ordering of time points, and determined the p-value using the number of genes assigned to the model profile and the expected number of assigned genes (adjusted p-value, 0.05 by Bonferroni correction). The boxes in the figures were colored if the profiles were statistically significant. We also used DAVID v6.8 to identify the GO biological processes involving genes with significant profiles.

PPI network construction, module analysis, and lncRNA-mRNA co-expression network analyses

The PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) (version 10.5) online database38. Interactions between proteins may offer insights into pathogenic mechanisms. In this study, the PPI network of DEGs from profiles 45 and 4 was constructed using the STRING database, and interaction with a combined score >0.4 was considered statistically significant. Cytoscape (version 3.6.1) is an open source bioinformatics software platform for visualizing molecular interaction networks39. The Cytoscape plug-in Molecular Complex Detection (MCODE) (version 1.5.1) is an algorithm for clustering a given network based on the topology to identify densely connected regions40. The PPI networks were constructed using Cytoscape, and the most significant module in the PPI networks was identified using MCODE. The selection criteria were as follows: degree cut-off = 2, MCODE scores >7, Max depth = 100, node score cut-off = 0.2, and k-score = 2. It is believed that genes with the same biological function or that regulate the same pathway may have similar expression patterns. The lncRNA-mRNA co-expression network was built based on the Pearson correlation coefficient (PCC) analysis of the lncRNA and mRNA expression levels41. The PCC was calculated for each lncRNA-mRNA, mRNA-mRNA, and lncRNA-lncRNA pair using package diffcoexp in R software42, and significant lncRNA-mRNA pairs with the absolute PCC > 0.9 and q-value < 0.01 were selected to construct the co-expression network using the Cytoscape 3.6.1 program.

Quantitative real-time PCR

Total RNA was extracted using TRIzol reagent followed by reverse transcription and quantitative real-time polymerase chain reaction (qRT-PCR) using the Power SYBR Green PCR Mastermix (Applied Biosystems, Foster City, CA, USA) as described previously43. The mean expression values were calculated based on the mean expression of the housekeeping gene β-actin. The primers used in this research are listed in Table 3. Statistical analyses were done using one-way analysis of variance, followed by Student–Newman–Keuls posthoc tests using SPSS 16.0 software (SPSS Inc., Chicago, IL, USA).

Table 3 The qRT-PCR primers used in this research.
Full size table