Differential transcription profiles of long non-coding RNAs in primary human brain microvascular endothelial cells in response to meningitic Escherichia coli

Abstract

Accumulating studies have indicated the influence of long non-coding RNAs (lncRNAs) on various biological processes as well as disease development and progression. However, the lncRNAs involved in bacterial meningitis and their regulatory effects are largely unknown. By RNA-sequencing, the transcriptional profiles of host lncRNAs in primary human brain microvascular endothelial cells (hBMECs) in response to meningitic Escherichia coli were demonstrated. Here, 25,257 lncRNAs were identified, including 24,645 annotated lncRNAs and 612 newly found ones. A total of 895 lncRNAs exhibited significant differences upon infection, among which 382 were upregulated and 513 were downregulated (≥2-fold, p < 0.05). Via bioinformatic analysis, the features of these lncRNAs, their possible functions, and the potential regulatory relationships between lncRNAs and mRNAs were predicted. Moreover, we compared the transcriptional specificity of these differential lncRNAs among hBMECs, human astrocyte cell U251, and human umbilical vein endothelial cells, and demonstrated the novel regulatory effects of proinflammatory cytokines on these differential lncRNAs. To our knowledge, this is the first time the transcriptional profiles of host lncRNAs involved in E. coli-induced meningitis have been reported, which shall provide novel insight into the regulatory mechanisms behind bacterial meningitis involving lncRNAs, and contribute to better prevention and therapy of CNS infection.

Introduction

The central nervous system (CNS) is a shielded environment, which is protected by a variety of physiological barriers. Bacterial meningitis is a serious disease involving CNS infection, with high mortality and morbidity, and is often associated with a poor prognosis and lifelong sequelae in the survivors¹. The blood–brain barrier (BBB) is composed of brain microvascular endothelial cells (BMECs), astrocytes, pericyte, and microglia, and maintains CNS homeostasis by regulating the transport of nutrients, molecules, as well as certain cells from the blood to the brain^2,3,4. The formation of tight junctions is a major characteristic of BMECs, which largely blocks the entry of circulating bacterial pathogens, and is widely reported to be critical for the structure as well as the function of the BBB⁵. However, some bacteria, such as Streptococcus pneumoniae, Neisseria meningitidis, Haemophilus influenzae type b, Group B streptococcus and Escherichia coli K1, can nonetheless break through the BBB via different strategies, thereby invading the CNS and causing bacterial meningitis^{6,7,8,9,10,11}.

Accumulating evidence suggests that CNS-infecting bacteria enhance BBB permeability by damaging the tight junctions of BMECs¹², and several molecules have been reported to mediate this enhancement of BBB permeability, such as the matrix metallopeptidase (MMP) family, transforming growth factor (TGF) β1, vascular endothelial growth factor (VEGF) A, interleukin (IL)-1β, IL-6 and tumour necrosis factor (TNF)-α^13,14,15. In the final stage of this process, a large number of inflammatory cells gain access to the cerebral parenchyma, triggering central inflammatory storm and CNS dysfunction¹⁶. Therefore, there is a need to reveal the mechanism underlying the penetration of the BBB by CNS-infecting bacteria, especially regarding how they interact with the host BBB and regulate host targets, thereby contributing to infection. Numerous studies have reported CNS-infecting bacterial invasion of the BBB, the receptors on BMECs with which these bacteria interact, the regulation of host molecules that mediate BBB disruption, as well as characterisation of microRNAs that regulate tight junction expression¹⁷. However, to our knowledge, the potential long non-coding RNAs (lncRNAs) and their possible regulatory mechanisms in CNS-infecting bacterial invasion of the host BBB have hardly been identified, but the discovery of such lncRNAs is beneficial for understanding their regulatory mechanisms in central infection.

Since the beginning of the twenty-first century, an increasing number of studies have indicated that non-coding RNAs have a wide variety of biological functions. It is generally believed that more than 40% of the genome is transcribed into mRNAs, of which less than 2% can encode proteins¹⁸. Non-coding RNAs are classified according to their length, with those less than 200 nucleotides (nt) being classified as microRNAs, while the remainder are lncRNAs¹⁹. lncRNAs were initially believed to have no biological function, but recent studies have suggested that they are widely involved in X-chromosome silencing, genomic imprinting, chromatin modification, transcriptional activation, interference, nuclear transport and other important regulatory processes^20,21,22,23. By RNA sequencing (RNA-seq) technology and bioinformatic analysis, thousands of lncRNAs have been discovered in organisms from fruit flies to humans and their functions have been predicted. However, the functions of only a small number of lncRNAs have been verified in the laboratory, largely within the field of oncology²⁴.

Research has also shown that the transcription of lncRNAs has tissue specificity and plays important roles in brain development²⁵. Associations of the expression, sequence and configuration of lncRNAs, as well as their abnormal binding with proteins, have also been revealed to be associated with CNS diseases²⁶. In the case of bacterial meningitis, it has been reported that microRNAs can affect the integrity of the BBB by regulating the expression of tight junctions²⁷. However, little is known about the involvement of lncRNAs in BBB damage and brain pathology. In addition, the existence of lncRNAs involved in bacterial infection of the CNS, as well as their potential molecular regulatory mechanisms, is completely unclear.

To fully elucidate the interaction between CNS-infecting bacteria and the host BBB, as well as the regulatory mechanisms involved in this, in the current study, we applied RNA-seq and bioinformatic approaches to identify potential host lncRNAs active in primary hBMECs in response to the infection of meningitic E. coli strain PCN033, a brain isolate that has been demonstrated to disrupt the BBB as well as inducing CNS inflammatory responses²⁸. The characteristics, potential functions, and transcriptional specificity of these lncRNAs were analyzed, and the possible correlation between lncRNAs and mRNAs was predicted. Moreover, we demonstrated significant regulatory effects of the proinflammatory cytokines on meningitic E. coli induction of lncRNAs. These observations together suggest the novel concept that lncRNAs are involved in bacterial meningitis, which should provide more regulatory host targets and contribute to further study on the pathogenic mechanisms involved in this disease, as well as better prevention and therapy for it.

Results

Identification of differentially expressed lncRNAs and mRNAs in primary hBMECs upon meningitic E. coli infection

To analyze the differences in transcriptomic profiles in primary hBMECs in response to infection, total RNAs from infected or uninfected primary hBMECs were prepared for RNA-seq following the workflow shown in Fig. 1. The RNA-seq data were generated from three uninfected and three infected groups, in which 91.7–99.3 million raw reads and 88.2–89.8 million clear reads per sample were obtained. The base percentage composition and quality distribution of each sample were analyzed (Supplemental Fig. 1). A total of 43,262 RNAs were obtained, among which 25,257 were lncRNAs and 18,005 were coding mRNAs (Table 1). Among these 25,257 lncRNAs, 24,645 lncRNAs had previously been annotated and 612 were potential novel lncRNAs that contained 757 transcripts (Supplemental Table 1 and Supplemental Fig. 2). A total of 1261 RNAs were found to have significant changes in their expression levels (increase by ≥2-fold or decrease to ≤0.5-fold, at p ≤ 0.05), which included 895 lncRNAs and 366 mRNAs (Table 1). The heat maps showed the trends of change of the lncRNAs as well as mRNAs in primary hBMECs upon infection (Fig. 2), among which the levels of 382 lncRNAs increased by ≥2-fold and 513 decreased to ≤0.5-fold. The most significantly upregulated lncRNAs (by >5-fold) and the most significantly downregulated lncRNAs (to <0.125-fold) are listed in Tables 2 and 3. Among the novel unannotated lncRNAs, XLOC_026572 (3.09-fold) and XLOC_026285 (2.75-fold) were the most highly upregulated, while XLOC_175351 (0.37-fold) and XLOC_070887 (0.37-fold) were the most downregulated (Supplemental Table 2).

Table 1 Overview of RNA-seq data in this study.

Table 2 The most significantly up-regulated lncRNAs in primary hBMEC upon infection (>5-fold).

Table 3 The most significantly down-regulated lncRNAs in primary hBMEC upon infection (<0.125-fold).

Figure 1

Overview of the RNA-sequencing processes and bioinformatic analysis.

Figure 2

Heat maps showing unsupervised clustering of meningitic E. coli-infected primary hBMECs and uninfected primary hBMECs based on differentially expressed lncRNAs (A) and coding genes (B). The expression profiles are displayed with three samples in each group. Red represents high and green represents low relative expression. One-way analysis of variance was used for the statistical analysis. LncRNAs or mRNAs with an increase in expression of ≥2-fold or a decrease to ≤0.5-fold and with a false discovery rate-adjusted p value ≤ 0.05 were considered statistically significant.

We also demonstrated a total of 366 coding mRNAs for which the expression increased by ≥2-fold or decreased to ≤0.5-fold, among which 198 mRNAs were significantly upregulated and 168 downregulated (Table 1). Detailed information on the differentially expressed mRNAs is listed in Supplemental Table 3. The GO, WEGO term, and KEGG pathway analyzes highlighted some substantially enriched pathways in this profile, including toll-like receptor, p53, NF-kappa B, MAPK, NOD-like receptor, HTLV-I infection, osteoclast differentiation and other important signalling pathways (Supplemental Fig. 3), many of which are closely associated with cell growth, differentiation, environmental adaptation, inflammatory responses and other important cellular physiological/pathological processes.

LncRNA Characterization

We observed that the lncRNA transcripts from primary hBMEC are shorter than the mRNAs (Fig. 3A), and the lncRNAs also tend to contain fewer exons (Fig. 3B). Both lncRNAs and mRNAs have various numbers of transcripts. The majority of lncRNAs identified in our study contain one to two transcripts, while most mRNAs have one to four (Fig. 3C). The open reading frame lengths of mRNAs were greater than those of lncRNAs (Fig. 3D). Fifty-three of the newly identified lncRNAs were differentially transcribed with an increase by ≥2-fold or a decrease to ≤0.5-fold, accounting for 8.7% of the total (Supplemental Table 2), but only 3.4% of the annotated lncRNAs exhibited significant changes (Supplemental Table 4). The overall distribution and number of these differentially expressed lncRNAs and mRNAs are displayed in a volcano plot (Fig. 3E). Box plots were also created to show the characteristics of lncRNA and mRNA expression in each sample (Fig. 3F). Using PhastCon software, we found that lncRNAs were less conserved than mRNA transcripts and the sequences of the newly identified lncRNAs are conserved at similar levels to those of the annotated lncRNAs (Fig. 3G). The level of conservation of lncRNAs was also analyzed among different species, including human, mouse, rat and fruit fly; we observed that lncRNAs, as well as mRNAs of human and mouse, were more similar to each other than those of rat and fruit fly (Fig. 3H).

Figure 3

LncRNA genomic features in primary hBMECs.

(A) Transcript lengths of 25,257 lncRNAs (red) and 18,005 mRNAs (green). (B) Exon numbers of lncRNAs and mRNAs. (C) Transcript numbers of lncRNAs and mRNAs. (D) Open reading frame length distribution of lncRNAs and mRNAs. (E) Volcano plots of the overall distribution of lncRNAs and mRNAs. (F) Box plots showing the expression density distribution of lncRNAs and mRNAs in each sample. (G) The cumulative distributions of median PhastCons score for various kinds of transcript. (H) Conservation score comparison for lncRNAs and mRNAs in human, mouse, rat and fruit fly.

Functional annotation and correlation analysis of the LncRNAs

Since lncRNAs largely function in gene regulation, we next annotated the lncRNAs identified in this study based on the function of the genes that they regulate, by trans target gene analysis, cis target gene analysis and competing endogenous RNA (ceRNA) analysis, etc.²⁹. It has been reported that lncRNAs can regulate overlapping and neighbouring genes³⁰. Considering the correlation between lncRNA and genes, we selected genes that overlap or are 100 kb upstream or downstream of the lncRNAs as candidate target genes regulated by lncRNAs. To investigate the relationship between lncRNAs and their neighbouring coding genes, we analyzed gene pairs formed by lncRNAs and their neighbouring genes, and identified 15,957 coding gene/coding gene pairs (2,467 in divergent) and 40,237 lncRNA/coding gene pairs (5997 in divergent) (Fig. 4A). Divergent lncRNAs were a special kind of lncRNAs that are transcribed in the opposite direction to their nearby coding genes, and could regulate the expression of their chromosomal neighbour genes. For example, recent study has indicated that divergent lncRNAs could regulate their adjacent mRNA in pluripotent cells³⁰. Among these coding genes near lncRNAs, there was a relative higher number of genes involved in the HIF-1 signalling pathway, MAPK signalling pathway, regulation of actin cytoskeleton, as well as tight junction (Fig. 4B). The novel lncRNAs and their associated mRNAs identified by cis target gene analysis are shown in Supplemental Table 5, while the annotated lncRNAs and their associated mRNAs are shown in Supplemental Table 6. The lncRNAs and their associated mRNAs, both of which exhibited significant changes in their expression levels, are listed in Table 4; the correlation coefficient between lnc-ANKRD37-1 and ANKRD37 was the highest.

Table 4 The correlation of both significantly changed lncRNAs and mRNAs (correlation > 0.95).

Figure 4

Analysis of lncRNA functional annotation.

(A) Number of gene pairs formed by lncRNAs and mRNAs located close together in the genome. (B) KEGG annotation of the functions of genes in the vicinity of lncRNAs. Red shows higher expression and blue shows lower expression. The number of differentially expressed genes in the signalling pathway is shown in parentheses.

Apart for the above-mentioned genes near lncRNAs, lncRNAs can also regulate gene expression by a trans effect over a long distance³¹. The results for such lncRNAs are shown in Supplemental Table 7. Five pairs of genes, such as lnc-NEUROD2-1 and ZNF546, were identified in the trans analysis.

Small RNAs are a large class of regulatory molecules that are found in almost all organisms. A previous study suggested that the function of small RNAs is controlled by nearby lncRNAs³². Here, small RNAs potentially associated with the lncRNAs and within 100 kb of them in the genome were also analyzed, the results of which are shown in Supplemental Table 8. Recent studies have shown that miRNAs can regulate not only gene expression, but also the transcription of lncRNAs²⁵. By using the miRNA target gene prediction software miRanda, we observed that large numbers of miRNAs can bind to different lncRNAs (Supplemental Table 9).

Recent genomic studies have suggested that some lncRNAs are precursors of miRNAs, for which the transformation into mature miRNAs occurs via cutting by the Drosha and Dicer enzymes³³. We thus compared our lncRNAs to miRBase to predict the possible miRNA precursors of lncRNAs, and found that seven lncRNAs have the potential to be precursors of miRNA (Supplemental Table 10).

Finally, lncRNAs can act as ceRNAs, which have a miRNA binding site and can regulate the level of miRNA, thus affecting its target genes³⁴. Here, via a series of bioinformatic methods, we analyzed the lncRNAs as well as mRNAs that possibly target the same miRNA and found that 6941 pairs of lncRNA and mRNA may have the same miRNA binding site (Supplemental Table 11).

Real-time polymerase chain reaction (PCR) verification of lncRNAs

Thirty differentially expressed lncRNAs with the most significant changes were selected for verification of their expression levels by quantitative real-time PCR. The results demonstrated that the upregulation (Fig. 5A) or downregulation (Fig. 5B) of the lncRNAs showed the same trends as in the RNA-seq data.

Figure 5

Quantitative PCR verification of the transcription of lncRNAs in primary hBMECs after meningitic E. coli infection.

The relative expression levels of the most highly differentially expressed lncRNAs as indicated were examined using qPCR in primary hBMECs that were treated with meningitic E. coli for 3 h. Data from RNA-seq were highly consistent with the qPCR results. Data are expressed as mean ± SEM from three separate experiments.

Proinflammatory cytokines differentially regulate lncRNAs in primary hBMECs

It was previously reported that cytokines can regulate the transcription of lncRNAs; for example, the transcription of lncRNAs can be affected by IL-4 and interferon-γ in macrophages³⁵. Our early work has also demonstrated that the infection of meningitic E. coli PCN033 can induce a strong inflammatory response in vivo and in vitro²⁸. We therefore investigated the effects of stimulation by diverse cytokines on the transcription of lncRNAs in primary hBMECs. Proinflammatory cytokines, including TNF-α, IL-1β, IL-6, IL-8, GRO-α and MIP-2, were used at a concentration of 10 ng/mL. As shown in Fig. 6, lnc-ANKRD37-1, lnc-RAB11B-3, lnc-PTTG1-1 and lnc-BIRC3-1, the levels of which were significantly increased in response to meningitic E. coli infection, were differentially regulated by the stimulation of different cytokines. GRO-α and MIP-2 significantly decreased the transcription of lnc-ANKRD37-1; TNF-α and IL-1β significantly increased the transcription of lnc-PTTG1-1 and lnc-BIRC3-1, while IL-6 significantly decreased the level of lnc-PTTG1-1. Notably, lnc-RAB11B-3 was downregulated by most of these inflammatory cytokines. Likewise, lnc-PRSS16-1, lnc-OLFML3-5 and lnc-FAM21A-2, which were all downregulated by bacterial challenge in our study, also exhibited different types of regulation in response to cytokines. For example, TNF-α and IL-6 positively regulated lnc-PRSS16-1, IL-1β and GRO-α downregulated lnc-OLFML3-5, and IL-1β significantly decreased the transcription of lnc-FAM21A-2. Taken together, these findings suggest that the proinflammatory cytokines and chemokines could differentially regulate the host lncRNAs. However, the specific regulatory mechanisms of these inflammatory factors on host lncRNAs, as well as their corresponding biological effects, require further investigation.

Figure 6

IL-8, MIP-2, GRO-α, IL-1β, IL-6 and TNF-α alter the expression of lncRNAs in primary hBMECs.

Cytokines at a concentration of 10 ng/mL were used to stimulate primary hBMECs for 24 h, and the transcription of these lncRNAs were investigated via qPCR. Data are expressed as mean ± SEM from three separate experiments.

Cell-type-specificity analysis of the differentially expressed lncRNAs

We next investigated the transcriptional specificity of the differentially expressed lncRNAs among different host cells in response to infection. Since our work was conducted on hBMEC, which is a cerebral microvascular endothelial cell, we thus chose another cerebral-born cell line U251, a human astrocyte, to analyze the transcription of lncRNAs. Meanwhile, since hBMEC is the brain vascular endothelial cell, we therefore tested the transcription of lncRNAs in a different, not brain-derived vascular endothelial cell HUVEC, a human umbilical vein endothelial cell line. We found that lnc-ANKRD37-1, lnc-CDC6-3, lnc-C5-1, lnc-BIRC3-1 and lnc-RP11-582J16.5.1-2 were upregulated in all of these three cell lines (by ≥2-fold) (Figs 5A and 7A,C), but exhibited responses of various intensities to meningitic E. coli infection. We moreover observed that the levels of lnc-CXCL3-1, lnc-PERP-10, lnc-PTTG1-1, lnc-IL5-1 and lnc-RSPH9-4 were increased in both brain-derived cells (≥2-fold): hBMECs and U251 cells, but not in the peripheral endothelial cells HUVECs, suggesting that these differentially expressed lncRNAs might be specifically produced in the brain upon infection. In addition, lnc-RAB11B-3, lnc-RAB11B-2, lnc-RAB11B-1 and lnc-KRT80-4 were found to be induced in hBMECs and HUVECs (≥2-fold) but not in U251 cells (Figs 5A and 7A,C), suggesting that these differentially expressed lncRNAs might typically be transcribed in endothelial cells. More importantly, we found that lnc-SMNDC1-1 was significantly induced in hBMECs but not in U251 cells and HUVECs, suggesting that this lncRNA is involved in meningitic E. coli penetration of the BBB (Figs 5A and 7A,C).

Figure 7

Transcriptional levels of lncRNAs in U251 and HUVEC cells in response to meningitic E. coli challenge.

The relative expression levels of the 30 most differentially expressed lncRNAs as indicated were examined using qPCR. Lnc-CXCL3-1 and lnc-C5-1 was the most significantly upregulated lncRNAs in U251 cells. Lnc-ITGA11-1 and lnc-DIRC1-1 almost maintained their expression levels in U251 cells. In contrast, lnc-ANKRD37-1 and lnc-RAB11B-1 were the lncRNAs that increased the most in HUVECs, but lnc-CXCL3-1 and lnc-PTTG1-1 were almost unchanged in these cells. Similarly, all the downregulated lncRNAs in primary hBMECs were significantly decreased in HUVECs upon infection. The lnc-SMNDC1-1 was undetectable in both U251 and HUVEC cells, which was marked by the empty bar graphs with # symbol. Data are expressed as mean ± SEM from three separate experiments.

We also investigated the lncRNAs with significantly decreased expression in U251 cells and HUVECs in response to infection. We identified only two lncRNAs in U251 cells, lnc-ITGA11-1 and lnc-DIRC1-1, for which the levels were not affected by meningitic E. coli infection, while the remaining ones exhibited similar decreases in hBMECs, U251 cells and HUVECs (Fig. 7B and D). This suggests that the downregulation of these lncRNAs was not cell-specific in response to meningitic E. coli.

Discussion

The aim of this study is to explore the potential key host lncRNAs that are involved in bacterial meningitis. Via RNA-seq approach, we identified and characterized the differentially expressed lncRNAs in primary hBMECs in response to meningitic E. coli. A total of 25,257 lncRNAs were identified here, 24,645 of which have already been annotated and 612 of which are potentially novel. Among these, the expression levels of 895 lncRNAs were significantly changed in response to infection, among which 382 lncRNAs were significantly increased and 513 were decreased. We also investigated the potential regulatory effects of the proinflammatory cytokines on the transcription of lncRNAs, as well as lncRNA specificity in different cells. To our knowledge, this is the first demonstration of the differential induction of host lncRNAs in primary hBMECs by meningitic E. coli, implying that lncRNAs may have potential regulatory roles in the course of bacterial meningitis.

By characterizing the differentially expressed lncRNAs, we found that they possess certain characteristics in common with those reported previously. For example, the lncRNAs are shorter than mRNAs, with fewer exons, fewer transcripts and less conservation than protein-coding mRNAs. We also found that the sequences of the newly identified lncRNAs are conserved at similar levels to those of the annotated lncRNAs. Notably, via correlation analysis, we demonstrated a potential association of these lncRNAs with miRNAs. The miRNAs have been shown to play a critical role in regulating the biological function of endothelial cells, and previous research indicated that miRNAs can affect the integrity of the BBB and blood–tumour barrier by decreasing the expression of tight junctions. For example, miR-155, miR-18a and miR-34c were shown to be able to downregulate ZO-1, occludin and claudin-5 expression, thereby, enhancing barrier permeability^36,37,38. LncRNAs were also reported to be involved in brain development as well as in some CNS diseases^39,40.

Recent studies have demonstrated that lncRNAs can regulate the expression of their neighbouring genes³⁰. It is well known that tight junction expression is regulated by many factors^41,42, but whether lncRNAs are directly involved in this regulation has remained unclear. However, in our correlation results, lnc-TAF9-1 and occludin mRNA were shown to potentially be associated with tight junction regulation (Supplemental Table 6), and our previous data indicate that meningitic E. coli induction of VEGFA and Snail-1 downregulated the expression of tight junctions²⁸. However, whether VEGFA and Snail-1 can be regulated by host lncRNAs was still unclear. Fortunately, our data reveal an extremely strong correlation between lnc-RSPH9-4 and VEGFA (Table 4), implying that lnc-RSPH9-4 may involve in the induction of VEGFA by meningitic E. coli. A recent study also showed that lnc-SPRY4-IT1 regulates intestinal epithelial barrier function by modulating the expression levels of tight junctions⁴³. In embryonic zebrafish, lnc-tie-1AS directly binds mRNA tie-1 and inhibits its transcription, leading to specific defects in vascular endothelial adherens junctions and tight junctions. In addition, lnc-tie-1AS significantly inhibits VEGF-induced endothelial tube formation in HUVEC cultures by decreasing the level of Tie-1 protein and damages the endothelial cell junctions⁴⁴. Here, we propose that lncRSPH9-4 regulates bacterial-induced VEGFA, which was shown by our previous study to regulate tight junctions. However, further investigation is required to determine how this lncRNA regulates VEGFA, and even the expression of tight junctions. The adhesion molecules VCAM-1 and ICAM-1 have been reported to influence transcellular or paracellular T-cell diapedesis across the BBB⁴⁵. In a study of portal vein tumour thrombus, lnc-ICR was reported to contribute to the development of disease by regulating ICAM-1⁴⁶. Here, our correlation analysis predicted the possible association of lnc-ZGLP1-3 with ICAM-1, suggesting that the former is a regulator mediating the inflammatory response in the CNS. However, this hypothesis requires verification. In the literature, there is also support for the association of lncRNAs with small RNAs; for example, lnc-TUG1 is known to regulate blood–tumour barrier permeability by targeting miR-144. Knockdown of lnc-MALAT1 increases blood–tumour barrier permeability by upregulating miR-140^47,48. The data obtained in this study imply that there might be a regulatory relationship between the identified lncRNAs and small RNAs. Although the primary finding in this report is the differential induction of host lncRNAs in primary hBMECs in response to meningitic E. coli infection, there is a need to undertake further investigation of the specific correlation with a certain lncRNA, as well as its regulatory effect in the process of bacterial meningitis.

It has been reported that cytokines can regulate the transcription of lncRNAs³². For example, in Kawasaki disease, TNF-α induces vascular endothelial cell apoptosis by causing the overexpression of pregnancy-induced non-coding RNA⁴⁹. Here, E. coli-induced cytokines, such as TNF-α, IL-1β, IL-6, IL-8, GRO-α and MIP-2, were selected to stimulate the primary hBMECs. Upon stimulation, we observed various regulatory effects of these cytokines and chemokines on the lncRNAs. Significantly, both TNF-α and IL-1β induced marked upregulation of lnc-PTTG1-1 and lnc-BIRC3-1. In terms of the functions of the genes associated with these lncRNAs, PTTG1 is an important oncogenic transcriptional factor implicated in various malignancies involved in epithelial–mesenchymal transition⁵⁰. BIRC3 is closely related to the apoptosis inhibitory factor family and plays pivotal roles in the regulation of apoptosis and NF-κB signalling⁵¹. Given this background, the induction of lnc-PTTG1-1 and lnc-BIRC3-1 by cytokines, the specific regulatory mechanisms involved, as well as their associated phenotypes in hBMECs are under our investigations. Meanwhile, our correlation analysis indicated that there is a strong correlation between lncRNAs and cytokines. For example, lnc-CXCL3-1 and GRO-3 exhibited synchronous increases of more than 20-fold, with a correlation coefficient of up to 0.989, providing a potential novel direction for future research.

Previous studies have suggested that lncRNAs exhibit tissue-specific and species-specific expression patterns in multiple diseases or development models^52,53. In this study, we demonstrated that lncRNAs exhibit various transcription patterns in primary hBMECs upon meningitic E. coli challenge; those lncRNAs with the greatest changes in expression were also confirmed by qPCR. We also investigated these significantly differentially expressed lncRNAs in another two cell lines, U251 and HUVECs, to analyze their transcriptional specificity. Lnc-CXCL3-1, lnc-PERP-10, lnc-PTTG1-1, lnc-IL5-1 and lnc-RSPH9-4 were excessively transcribed in cerebral cell lines, primary hBMECs and U251, but almost unchanged in HUVECs. In terms of the functions of the genes associated with these lncRNAs, CXCL3 is a member of the chemokine family and involved in many pathophysiological processes, such as inducing leukocyte migration to specific inflammatory sites and promoting carcinogenesis⁵⁴. PERP was first identified as a p53 effector, and has since been shown to play roles in development, caspase activation and cancer⁵⁵. IL-5, one of the most important interleukins, is mainly secreted by macrophages and T cells. It has been confirmed to have multiple biological functions, including promoting the differentiation of B-1 cells⁵⁶. In contrast, lnc-RAB11B-1, lnc-RAB11B-2 and lnc-RAB11B-3, which are all around the mRNA Rab11, were shown to be upregulated in two kinds of endothelial cells, primary hBMECs and HUVECs, but not in U251 cells. Rab11 is reported to be associated with recycling endosomes and has been confirmed to regulate vesicular trafficking through the recycling of endosomal compartments, the plasma membrane and early endosomes to the trans-Golgi network⁵⁷. In the literature, it is indicated that meningitic E. coli transmigrates through the BBB enclosed in a vacuole without multiplying intracellularly, and it must maintain the ability to traffic in BMECs as live bacteria, during which it must prevent lysosomal fusion, thereby avoiding lysosomal enzyme degradation¹⁶. However, it is currently poorly understood whether Rab11 is required in this process and which regulatory mechanisms involving Rab11 may be associated with this. Notably, we observed the specific induction of lnc-SMNDC1-1 in hBMECs, but not in the other two cell lines, in response to meningitic E. coli. The level of lnc-SMNDC1-1 increased more than 8-fold upon infection, and we were able to correlate this lncRNA with the protein-coding gene MXI1. However, the specific role of lnc-SMNDC1-1 and its possible regulatory effects, as well as the mechanisms involved, remain to be uncovered.

In conclusion, these findings indicate that lncRNAs are endogenous regulators that are responsible for certain pathophysiological phenotypes in bacterial meningitis. However, the precise regulatory mechanisms involving individual lncRNAs require more experimental evidence. In previous studies, the differential lncRNA transcriptomic profiles in mouse BMECs after cerebral ischemia were reported⁵⁸. In addition, in a stroke model, the regulation of lncRNAs in the BBB was studied⁴⁰. However, to the best of our knowledge, the host lncRNA transcriptomic profiles upon meningitic bacterial penetration of the BBB have not been reported. This is the first attempt via high-throughput RNA sequencing technology to characterize the differentially induced lncRNAs in primary hBMECs in response to meningitic E. coli challenge. Proinflammatory cytokines and chemokines have been shown to be the potential stimuli that mediate the transcription of lncRNAs in primary hBMECs. Moreover, we compared the transcription of these differentially expressed lncRNAs in human primary hBMECs, U251 cells and HUVECs in vitro, and found that some of the lncRNAs exhibited different expression profiles among these cell lines. Taking these findings together, the discovery of host lncRNAs, as well as novel lncRNA transcripts, provides new insights into the process of bacterial meningitis infection, and suggests potential new targets to improve the prevention of this disease in the future.

Methods

Bacterial strain and cell culture

The meningitic E. coli strain used in this study, PCN033, was isolated from swine cerebrospinal fluid in Hunan Province, China, in 2006⁵⁹. The strain was routinely grown aerobically at 37 °C in Luria–Bertani medium overnight.

Primary hBMECs which we used in this study were kept within ten passages were kindly provided by Prof. Kwang Sik Kim from Johns Hopkins University School of Medicine and routinely cultured in RPMI1640 supplemented with 10% heat-inactivated foetal bovine serum (FBS), 10% Nu-Serum, 2 mM L-glutamine, 1 mM sodium pyruvate, nonessential amino acids, vitamins, and penicillin and streptomycin (100 U/mL)⁶⁰. The human astrocyte cell line U251 was cultured in Dulbecco’s modified Eagle’s medium supplemented with 10% heat-inactivated FBS, and HUVECs were maintained in EGM-2 medium (Lonza, Walkersville, MD, USA) with 10% heat-inactivated FBS. All of the cells were incubated in a 37 °C incubator under 5% CO₂ until monolayer confluence. Confluent cells were washed three times with Hanks’ balanced salt solution (Corning Cellgro, Manassas, VA, USA) and starved in serum-free medium (1:1 mixture of Ham’s F-12 and M-199) for 16–18 h before further treatment.

In vitro infection and cytokines stimulation assays

Meningitic E. coli strain infection of primary hBMECs, as well as U251 cells and HUVECs, was performed in accordance with previously described methods^61,62. Briefly, E. coli overnight culture was resuspended and diluted in serum-free medium and added to the starved confluent primary hBMEC monolayer grown in 10-cm dishes at a multiplicity of infection of 10 (approximately 10⁸ colony-forming units per dish) to allow invasion at 37 °C for 3 h. For cytokine stimulation, recombinant human IL-8, MIP-2, GRO-α, IL-1β, IL-6 and TNF-α were purchased from Novoprotein Scientific (Shanghai, China) and used at a final concentration of 10 ng/mL to stimulate the primary hBMECs for 24 h. Finally, cells were washed three times with chilled PBS and subjected to RNA extraction by using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), in accordance with the manufacturer’s instructions.

RNA sequencing and bioinformatic analysis

LncRNA library preparation

For library construction, ribosomal RNA was removed using the Ribo-zero rRNA Removal Kit to obtain approximately 2 μg of total RNA. The library was constructed with the NEB Next^® Ultra™ Directional RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, USA), following the manufacturer’s recommendations. Briefly, fragmentation was carried out using divalent cations at an elevated temperature in NEBNext First-Strand Synthesis Reaction Buffer (5×). First-Strand cDNA was synthesised using random hexamer primer and M-MuLV reverse transcriptase. Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In the reaction buffer, dNTPs with dTTP were replaced by dUTP. Remaining overhangs were converted into blunt ends via the activities of exonuclease/polymerase, followed by adenylation of the 3′ ends of DNA fragments. Finally, NEB Next Adaptor with a hairpin loop structure was ligated to prepare for hybridisation. Fragments of approximately 200 base pairs (bp) were selected and extracted using the AMPure XP system (Beckman Coulter, Beverly, MA, USA). Library quality was assessed on an Agilent Bioanalyzer 2100 system. The libraries were sequenced at the Total Genomics Solution Institute (Shenzhen, China) on an Illumina Hiseq 2500 platform and 125-bp paired-end reads were generated.

Raw reads quality control

Raw reads in fastq format were subjected to a process for removing adaptor and low-quality reads based on the following criteria: (1) reads that aligned to adaptors or primers with no more than two mismatches; (2) reads with >10% unknown bases (N bases); and (3) reads with >50% low-quality bases (quality value ≤10) in one read. Finally, the filtered reads were used for further analysis.

Novel LncRNAs prediction

High-quality reads were first aligned to the human genome (UCSC, hg19) with Hisat (v0.1.6) and then reconstructed to transcripts with StringTie (v1.0.4). To achieve “union” of the transcripts, cuffcompare (v2.1.1) was used to remove the replicated transcripts in different samples or conditions. Subsequently, several filtering steps were applied to isolate the most robust transcripts, in accordance with the method proposed by Prensner⁶³. First, transcripts with a total length of less than 200 nt were discarded. Second, transcripts for which the peak expression across all samples was less than 2.0 and those present in only one sample considered as a ‘background’ transcript were discarded. Third, transcripts overlapping with known mRNAs and lncRNAs on the same strand were removed. Finally, the remaining transcripts without protein-coding potential predicted by lncRNA prediction software (CPC v0.9-r2, CNCI v2.1 and PFAM) were classified as novel lncRNAs.

Analysis of the differential expression of lncRNAs and mRNAs

HTSeq v0.6.1 was used to count the number of reads mapped to each gene for quantification of the gene expression level and normalised using the FPKM (fragment per kilobase per million fragments) method. The R package DESeq2 (v1.4.5) was used to determine the differential expression. The resulting p-values were adjusted using Benjamini and Hochberg’s approach for controlling the false discovery rate⁶⁴, and genes or lncRNAs with an adjusted p-value of < 0.05 and an increase in their expression of ≥2-fold or a decrease to ≤0.5-fold were considered differentially expressed.

Analysis of cis and trans target gene

For cis role of lncRNA, the gene flanked within 100 kb of lncRNA was chosen. And then the pearson correlation coefficient was calculated between gene and lncRNA using the expression level (FPKM). Finally, the gene-lncRNA pairs with absolute value of the correlation coefficient more than 0.6 were considered as candidate lncRNAs that regulated corresponding cis target genes. For trans role of lncRNA, the pearson correlation coefficient was calculated between gene and lncRNA using the expression level (FPKM). The gene-lncRNA pairs with absolute value of the correlation coefficient more than 0.6 were chosen. And then filtered with sequence homology (blastn, E-value < 1.0E-10 and identity >99 and matched length ≥20 bp) between gene-lncRNA pair.

Competing endogenous RNA analysis

Based on the competing endogenous RNA hypothesis (“ceRNA hypothesis”), lncRNAs can function as mircoRNA “sponge” to interact directly with mircoRNAs and prevent them from binding to mRNAs. The miRNAs binding sites in lncRNAs were predicted with miRanda⁶⁵. And the miRNAs binding sites of mRNAs were download from miRTarBase database which collectes the experimentally validated microRNA-target interactions.

GO and KEGG enrichment analysis

The R package goseq (v1.16.2) was used to perform the GO enrichment analysis, and GO terms with a corrected p-value of above 0.05 were excluded. For KEGG analysis, the genes were mapped directly to the KEGG database. Then, the enriched pathways were obtained using a q-value cutoff of 0.05 with the R hypergeometric function and R q-value package.

Conservation analysis

PhastCons scores of human (hg19), mouse (mm10), rat (rn5) and fruit fly (dm6) were downloaded from the UCSC database. To assign a conservation score to a transcript, the mean PhastCons score for the exons and introns of each transcript model was calculated. The conservation score was compared among the protein-coding sequence, lncRNA and corresponding introns⁶⁶.

RNA isolation and quantitative real-time PCR analysis

Total RNA was extracted from primary hBMECs, U251 cells and HUVECs using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Contaminating genomic DNA was removed by DNase I treatment (NEB, Ipswich, MA, USA). Aliquots (1 μg) of the total RNA in each sample were subjected to cDNA synthesis using PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara Bio Inc., Tokyo, Japan). Real-time PCR was performed with a ViiA™ 7 Real-Time PCR System (Applied BioSystems, Foster City, CA, USA) using Power SYBR Green PCR Master Mix (Applied BioSystems), in accordance with the manufacturer’s instructions. Primers for the quantitative real-time PCR are listed in Supplemental Table 12. The amplification conditions were as follows: 50 °C for 2 min and 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. The products were then applied to a melt curve stage with denaturation at 95 °C for 15 s, annealing at 60 °C for 1 min and slow dissociation by ramping from 60 °C to 95 °C at 0.05 °C/s to ensure the specificity of the PCR products. The expression levels of the target genes were normalised to the expression of GAPDH. Each assay was performed independently in triplicate.

Statistical analysis

Data are expressed as mean ± standard error of the mean (SEM), and the significance of differences between groups was evaluated by two-way analysis of variance. P-values < 0.05 were considered significant. Graphs were plotted and analyzed using GraphPad Prism version 6.0 (GraphPad Software, La Jolla, CA, USA).