Skip to main content

RNA-seq co-expression network analysis reveals anxiolytic behavior of mice with Efnb2 knockout in parvalbumin+ neurons


Anxiety disorders are the most common psychiatric disorders, and the change in the activity of the prefrontal cortex (PFC) is considered as the underlying pathological mechanism. Parvalbumin-expressing (PV+) inhibition contributes to the overall activity of the PFC. However, the molecular mechanism underlying the excitation-inhibition imbalance of PV+ neurons in the PFC is unknown. Efnb2 is a membrane-bound molecule that plays an important role in the nervous system through binding the Eph receptor. To investigate whether the loss of Efnb2 in PV+ affects anxiety, we examined the behavior of wild type and Efnb2 in PV+ neurons knockout (KO) mice. We monitored the defensive responses to aversive stimuli of elevated plus maze (EPM) and found that KO mice exhibited obvious fearless and anxiolytic behaviors. To further investigate the underlying regulatory mechanism, we performed RNA sequencing, analyzed the differentially expressed genes (DEGs), and constructed the weighted gene co-expression network analysis (WGCNA). The WGCNA identified 12 characteristic modules. Among them, the MEgreen module showed the most significant correlation with KO mice of EPM stimuli. The Gene Ontology enrichment and the Kyoto Encyclopedia of Genes and Genomes enrichment analysis revealed that this was related to the distal axon, Ras signaling pathway and insulin signaling pathway. Furthermore, the whole-cell voltage clamp recordings also proved that Efnb2 gene knock-out could affect synaptic function. Together with the transcriptomic analysis of mice with Efnb2 knockout on PV+ neurons, our findings suggest that Efnb2 gene in the PV+ neuron of PFC may be a crucial factor for fear and anxiety, which provide an insight into anxiety pathophysiology.


Anxiety is one of the most common mental disorders that can produce adverse cognitive effects and contribute to the clinical manifestation of anxiety disorders. Anxiety disorders are of considerable morbidity [1]. They are commonly comorbid with other mental disorders like major depression, dysthymia, and so on. These disorders will undoubtedly increase the rates of disability and the risk of mortality. The balance of excitation and inhibition in prefrontal cortex (PFC) plays a fundamental role in the emotion regulation of anxiety [2]. The excitation and inhibition balance is regulated by the interaction of glutamatergic and GABAergic neurons. Specifically, the prefrontal GABAergic interneurons, especially the Parvalbumin-expressing (PV+) neurons, are the primary regulators of the projective excitatory neurons’ spiking activity and emotional behavior [3]. The PV+ neurons account for more than 50% of the total population of GABAergic neurons in the frontal cortex and control the region near or on the cell bodies of excitatory pyramidal cells by providing strong, fast spiking inhibition [4, 5]. Therefore, understanding the underlying mechanism remains a priority to promote the development of efficacious and reliable pharmacological anxiolytic therapies.

Ephrins are membrane-bound molecules that play an important role in the nervous system through the binding of the Eph receptor [6, 7]. Ephrin-Eph signaling is involved in mediating the communication among neurons and between neurons and glial cells [8]. For example, the Ephrin B2 (Efnb2) extracellular portion associates with NMDA neurotransmitter receptors and promotes their clustering at synapses following ephrin-B stimulation [9]. Efnb2, a more powerful player compared with ephrin b1 and b3, was found to be expressed in a large proportion of forebrain inhibitory neurons [10]. It has recently been implicated in GABAergic circuit development and considered as a potent regulator of neuronal development and synaptic function [11,12,13,14]. Important advances in the studies of Efnb2 in the developing nervous system included axon and dendrite growth [11, 12], dendritic spine morphogenesis, and synapse formation [13, 14]. Reverse signaling by postsynaptic Efnb2 plays an essential role in synaptic plasticity in mutant mice [8, 15]. However, the underlying mechanisms of the pathogenic variation of the Efnb2 in the anxiolytic effect is currently unclear. Combining the existing researches of PV+ neurons and Efnb2 [2, 9], we speculated that the Efnb2-mediated dysfunctions of specific PV+ neurons lead to the imbalanced neural networks.

Efnb2 is an important factor in the neurological system, but the association between Efnb2 gene and anxious stress remains controversial. In this paper, we knocked out the Efnb2 gene of PV+ interneurons in the mouse PFC brain region. To investigate whether genetic Efnb2 variants affect fear and anxiety behaviors, we tested the behavior of mice using the experiments of marble burying, cold stress defecation, and elevated plus maze (EPM) between wild-type (WT) and knock-out (KO) mice. The EPM is an animal model of fear and anxiety and is related to physiology, pharmacology, and behavior [16, 17]. Moreover, EPM has been used as a stress model to investigate the regulatory mechanism of anxiety [18]. Placing an animal on an elevated platform, especially the open arms of the maze, could induce stress in animals and hence increase fear and anxiety.

In this study, we deleted Efnb2 specifically in PV+ neurons and observed fearless and anxiolytic behaviors in KO mice. Furthermore, by performing RNA sequencing, weighted gene co-expression network analysis (WGCNA) and whole-cell patch-clamp recording, we aimed to determine gene co-expression networks underlying fearless and anxiolytic behaviors and identify key modules associating with it, which may provide a treatment strategy and therapeutic target to regulate fear and anxiety symptoms.


Mice with Efnb2 knockout on PV+ neurons show obvious anxiolytic behavior

We generated mice lacking Efnb2 only in PV+ neurons by crossing mice with Cre driven by the PV and mice carrying Efnb2 flanked LoxP sites (Fig. 1a). Cre-recombination and mouse genotypes were confirmed by RT-PCR analysis. The Efnb2LoxP primers resulted in gene products of 500 bp. The PVcre primers resulted in gene products of 400 bp. The Mutant-Efnb2 primers resulted in gene products of 350 bp. The mice with loxP-flanked alleles in Efnb2 were WT group, the mice with Efnb2 gene knockout were KO group (mice 1–3, 5, 7, 8; Fig. 1b).

Fig. 1

Behavioral results of wild-type (WT) and Efnb2 on PV+ neurons knock-out (KO) mice. a Schematic drawing of the breeding strategy for a generation of KO mice. b PCR-based analysis of mice genotypes. c Mouse movement track diagram. Black: open arms; grey: closed arms; red: tracks. d Results of entries into open arms (male, n = 30, 27; Mann Whitney test, #P ≤ 0.0001). e Results of time in open arms (male, n = 30, 27; Mann Whitney test, #P ≤ 0.0001). f marbles buried diagram. g marbles buried number (male, n = 12, 12; Mann Whitney test, #P ≤ 0.0001). h Fecal number (male, n = 12, 12; Mann Whitney test, **P ≤ 0.01)

To investigate whether ablation of Efnb2 gene on PV+ neurons was required for the anxiolytic effect elicited on PFC, we performed a behavioral examination through a series of testing trials with marble burying, cold stress defecation, and EPM between WT and KO mice. We chose these behavioral paradigms, because we were mostly interested in determining the innate anxiolytic role of Efnb2. Mice are free to avoid the aversive stimulus or approach the potential threat. These behavioral tests all showed that KO mice exhibited obvious fearless and anxiolytic behaviors. In the EPM experiments, both entries into open arms (%) and time spent in open arms (%) of KO mice were significantly higher than that of WT mice (Mann Whitney test, #P ≤ 0.0001; Fig. 1c, d, and e). In the marble burying and cold stress defecation experiments, KO mice had lower number of marbles buried (Mann Whitney test, #P ≤ 0.0001; Fig. 1f and g) and fecal number (Mann Whitney test, **P ≤ 0.01; Fig. 1h) than the WT group.

Differential expression genes among WT, KO, WT_EPM and KO_EPM

To further explore the molecular determinants of the neural network and behavioral abnormalities in Efnb2 knockout mice, RNA sequencing was performed on the mPFC of the WT (male; n = 4), KO (male; n = 4), WT mice after 1 h of the EPM test (WT_EPM; male; n = 3), and KO mice after 1 h of the EPM test (KO_EPM; male; n = 3). Among the samples, 32.35–53.19 million raw reads were generated. After removing low-quality reads, 31.60–50.92 million clean reads were obtained. Approximately 97.0% of clean reads per sample were successfully mapped to the mouse reference genome (Additional file 1). Finally, 55385 transcripts were obtained and the Fragments Per Kilobase of transcript per Million mapped reads (FPKM) were calculated. The principal component analysis (PCA) found a clear separation between the KO_EPM and WT_EPM, and an overlap was found between the KO and WT (Fig. 2a). We also analyzed the differentially expressed genes (DEGs; Additional file 2) between each group with Deseq2 package in R. Based on the statistical analysis (|Fold change|> 2, Padj < 0.05), 235 (KO vs. WT) genes were dysregulated, of which 116 were up-regulated in the KO group of mice. In other comparisons, 435 (KO_EPM vs. WT_EPM), 935 (KO_EPM vs. KO) and 1051 (WT_EPM vs. WT) genes were dysregulated. We found that the number of dysregulated genes induced by EPM was higher than of dysregulated genes induced by genotype (Fig. 2b).

Fig. 2

Analysis of differentially expressed genes (DEGs). a Principal component analysis (PCA) of WT, KO, WT_EPM, and KO_EPM. b The numbers of up-regulated (blue) and down-regulated (yellow) DEGs between KO vs. WT, KO_EPM vs. WT_EPM, KO_EPM vs. KO and WT_EPM vs. WT. c Venn diagram of DEGs (KO vs. WT) and DEGs (KO_EPM vs. WT_EPM). d The top 5 GO terms of 249 genes (DEG1), which were unshared DEGs in KO_EPM vs. WT_EPM. e Venn diagram of DEGs (KO_EPM vs. KO) and DEGs (WT_EPM vs. WT). f The top 5 GO terms of 507 genes (DEG3), which were unshared DEGs in KO_EPM vs. KO. GO terms are sorted by the P value. Gene number represents the number of gene types. BP biological process, CC cellular component, MF molecular function

The Venn diagram of DEGs (KO vs. WT) and DEGs (KO_EPM vs. WT_EPM) showed that there were 381 genes, which were unshared DEGs in KO_EPM vs. WT_EPM (DEG1, Fig. 2c). GO enrichment analysis of the 381 genes showed that they were enriched in the ribonuclease P RNA binding, Toll-like receptor binding, MHC protein binding, peptide binding, and others (Fig. 2d). There were 507 genes, which were unshared DEGs in KO_EPM vs. KO (DEG3, Fig. 2e). The GO enrichment results of the 507 genes revealed the actin binding, structural constituent of myelin sheath, PDZ domain binding, and others (Fig. 2f; Additional File 3).

WGCNA identify the co-expression module

To identify gene expression profiles associated with fearless and anxiolytic behaviors, we performed a WGCNA to the RNA-seq data [19, 20]. In WGCNA, 11,527 genes (Additional file 4) of the transcriptome dataset are clustered according to co-expression. Average linkage hierarchical clustering was applied to the topological overlap matrix. The branches of highly correlating genes were formed, which were cut and assigned a color. Finally, 12 modules (MEbrown, MEtan, MEblack, MEgreen, MEpurple, MEturquoise, MEred, MEgreenyellow, MEmagenta, MEyellow, MEblue, and MEpink) were identified (Fig. 3a, Additional file 5). The MEgrey module represents unassigned genes. The Module Eigengene (ME) was calculated for each module; it functioned as a representative of the module. The correlations among modules were calculated. Their heatmaps were shown in Fig. 3b. The biweight midcorrelation between each ME and sample trait (binary variables: WT, KO, EPM_WT, and EPM_KO) was calculated based on the bicorAndPvalue() function of WGCNA package in R. Among the 12 modules, the MEgreen module had the most significant correlation with the KO_EPM (Fig. 3c).

Fig. 3

Weighted gene co-expression network analysis (WGCNA). a Cluster dendrogram of co-expression networks. b Heatmap of WGCNA module-to-module correlation. c Heatmap of the correlation between WGCNA module eigengene and sample traits. Within each table cell, upper values indicate the biweight midcorrelation between ME and the traits. Lower values in brackets are the P-values for the association test

The analysis of MEgreen module

GO enrichment analysis of MEgreen module consisted of three parts, namely, biological process (Fig. 4a), cellular component (Fig. 4b), and molecular function (Fig. 4c). We selected the top 20 GO terms to sort the P value. The potential biological processes revealed by the MEgreen module included cellular component disassembly, regulation of cell growth, protein-containing complex disassembly, learning or memory, and exocytosis. The cellular components included distal axon, synaptic vesicle membrane, exocytic vesicle membrane, growth cone, and exocytic vesicle. The molecular function included small GTPase binding, nuclear import signal receptor activity, Ras GTPase binding, Ras guanyl-nucleotide exchange factor activity, guanyl-nucleotide exchange factor activity, and others. GO terms were related to neural regulation include distal axon, synaptic vesicle membrane, exocytic vesicle membrane, nuclear import signal receptor activity, growth cone, and others (Fig. 4a–c, Additional File 6). KEGG enrichment analysis of MEgreen module revealed the higher-level systemic functions that were potentially involved, including autophagy, insulin signaling pathway, Ras signaling pathway, and others (Fig. 4d, Additional File 7).

Fig. 4

The analysis of MEgreen module genes. a The top 20 GO terms relevant to the biological process of MEgreen module. b The top 20 GO terms relevant to the cellular component of MEgreen module. c The top 20 GO terms relevant to the molecular function of MEgreen module. d KEGG enrichment analysis of MEgreen module. The top 20 KEGG terms sorted by the P value. e Venn diagram of MEgreen module genes, DEG1 (unshared DEGs in KO_EPM vs. WT_EPM) and DEG3 (unshared DEGs in KO_EPM vs. KO). f The top 5 GO terms of the shared genes among MEgreen module, KO_EPM1 and KO_EPM2. GO terms are sorted by the P value. Gene number represents the number of gene types. BP biological process, CC cellular component, MF molecular function

Venn diagram showed that there were 17 genes which were shared among MEgreen module, DEG1 and DEG3 (Fig. 4e). GO enrichment analysis of the 17 genes showed that they were enriched in the voltage-gated anion channel activity, DNA polymerase activity, transmembrane transporter complex, regulation of AMPA receptor activity and others (Fig. 4f). The Sult1a1, Cnih3, Clcnka and Chtf8 were the main genes of the enrichment.

Efnb2 knockout on PV+ neurons effects on excitatory post-synaptic potential

To investigate whether Efnb2 knockout on PV+ neurons of PFC affects the neuronal transmission of pyramidal neurons in medial prefrontal cortex (mPFC), we recorded the spontaneous excitatory post-synaptic currents (sEPSCs) and spontaneous inhibitory post-synaptic currents (sIPSCs) of neurons in the PFC of WT and KO mice. In interneurons from WT mice, the average sEPSCs frequency was 0.5494 ± 0.0846 Hz (n = 11). The average sEPSCs frequency from WT mice was significantly reduced, i.e., 0.2572 ± 0.0566 Hz (n = 13, **P < 0.01; Fig. 5a and b). We investigated the overall incoming inhibitory events by recording sIPSCs of PFC neurons in WT and KO mice, and no difference was found (Fig. 5c and d).

Fig. 5

Frequency of spontaneous excitatory post-synaptic currents (sEPSCs) are reduced after Efnb2 knockout. a Example traces showing sEPSCs recorded in interneurons of PFC for WT (left) and KO (right). b Cumulative probability plot (left) for sEPSCs event interval. Average sEPSCs frequency for the two conditions (right). Data represent mean ± SEM. c Example traces showing spontaneous inhibitory post-synaptic currents (sIPSCs) recorded in interneurons of PFC for WT (left) and KO (right). d Cumulative probability plot (left) for the event interval of IPSP. Average sIPSCs frequency for the two conditions (right). Data represent mean ± SEM. (male, Unpaired t test, **P ≤ 0.01)


Although a previous study has shown that Efnb2 is expressed in PV+ neurons of PFC [10], its function and relationship with disease progression are unknown. The PV+ neurons of PFC are the primary regulators of emotional behavior [5]. Whether the Efnb2 gene on PV+ neurons can regulate the mental disorder has not been studied. We showed that Efnb2 gene in PV+ neurons of PFC could affect fear and anxiety. The Efnb2 knockout mice exhibited obvious fearless and anxiolytic behaviors.

The PCA results showed a significant difference between KO_EPM and other three groups. In order to further explore the molecular mechanism of fearless behavior, we analyzed the DEGs and identified the genes which were affected by the genotype and the EPM. The DEG1 which were mainly induced by the genotype in KO_EPM mice revealed the molecular function of Toll-like receptor binding. Toll-like receptors (TLRs) have recently emerged as regulators of neuronal survival and developmental neuroplasticity. And the TLR3-deficient mice could be anxious in EPM tasks [21]. The DEG3 which were mainly induced by the EPM in KO_EPM mice revealed the molecular function of actin filament binding. The actin signaling can be affected by the metabotropic glutamate receptors 7 (mGluR7). And the ablation of mGluR7 can lead to the reduced anxiety and depression behaviors related to PFC dysfunction [22]. We found that KO mice induced by EPM will cause the changes in the expression of a series of genes in the mPFC. And these genes were related to the anxiolytic behavior. Therefore, we surmised that knockout of the Efnb2 in PV+ neurons may affect the activity of PV neurons, which further affected the activity of the entire PFC, resulting in the changes in gene expression of KO_EPM.

The WGCNA was used to identify a distinct co-expression network module, which was associated with EPM stress caused by Efnb2 knockout, finding that MEgreen module genes had the most significant correlation with the KO_EPM group. The results of GO enrichment with genes of MEgreen module showed that it is associated with the distal axon, small GTPase binding, and synaptic vesicle membrane. KEGG pathway included Ras signaling and insulin signaling pathways. Previous studies have showed that axon terminals of cholecystokinin-immunoreactive basket cell expresses the type 1 cannabinoid receptor, and the interactions between cholecystokinin and the endogenous cannabinoid system in the basolateral nuclear complex of the amygdala could modulate anxiety-like behavior and fear learning/expression [23]. The Ras plays an essential role in various physiological processes consisting of several enzymes, peptides and receptors [24]. For example, the AT1 receptors could regulate activity of noradrenergic neurons in the brain. And the dysregulation in noradrenaline levels in the brain plays an important role in both anxiety and depression [25]. Insulin signaling accounts for the development of a variety of neuropsychiatric disorders through the interaction between the inflammatory condition and RAS, which impedes brain insulin pathway, resulting in neurobehavioral damage [26].

There were 17 genes which were shared among MEgreen module, DEG1 and DEG3. The main gene enriched by GO was Cnih3 among the 17 genes, which was related to transmembrane transporter complex, ion channel complex and regulation of AMPA receptor activity [27]. Previous studies have shown that AMPA receptors mediate fast excitatory synaptic transmission in the brain [28]. In order to further verify whether the anxiolytic behaviors in KO mice is accompanied by alterations in synaptic excitability, we performed whole-cell voltage clamp recordings in PFC brain region. The average sEPSCs frequencies of WT mice were significantly higher than those of KO mice. The finding is consistent with the results of previous studies. The dmPFC 5-HT6 receptor activation in mice decreased the anxiety. 5-HT6 receptor antagonist SB 271046 could increase sEPSC [29]. This undoubtedly provided synaptic-related evidence for the results of RNA-seq analysis. This study only focused on the changes of the average sEPSCs frequencies between KO and WT mice. The specific regulation mechanism is expected to be studied further in the future.

In conclusion, this article identified a distinct co-expression network module, which has the most significant correlation with KO_EPM by using WGCNA. We characterized the biological processes of this module. The module was related to the distal axon, Ras signaling pathway, and insulin signaling pathway. The whole-cell voltage clamp recordings also showed that Efnb2 gene knock-out could affect the average sEPSCs frequencies. This is the first study to identify the fear- and anxiety-associated genes through WGCNA in Efnb2 knockout mice. We confirmed that Efnb2 on PV+ neurons are important components of fear and anxiety behaviors, which suggested that the identification of the module associated with the fear and anxiety behaviors can provide a treatment strategy and therapeutic target to regulate fear and anxiety symptoms.

Materials and methods


All animal studies were performed in accordance with the National Institute of Health Guide for the Care and Use of Laboratory Animals and were approved by Animal Care and Use Committee in Shanghai Jiao Tong University. Mice were of a mixed CD1/B6 genetic background. All mice have been previously described as having loxP-flanked alleles in Efnb2 [30], which were combined with the parvalbumin (PV)-Cre mice [31].

All mice were housed and had access to food and water in a room with a standard 12 h dark/12 h light cycle. Lights were turned on at 7:00 AM and turned off at 7:00 PM. The room temperature was maintained at 24 ± 1 °C. Behavioral tests were performed between 9:00 AM and 4:00 PM. The behavioral experiments were conducted in the following order: marble burying, cold stress defecation, and EPM. A 3–5 day interval was set between tests. Before all experiments, mice were allowed to habituate in the testing room for 60 min.

Marble burying

Marble burying test was performed in a cage with the following dimensions: length, 48 cm; width, 35 cm; and height, 20 cm. The cage was filled with 7–8 cm of fresh, autoclaved wood chip bedding (Sharon et al., 2019). Before the test, each mouse was acclimatized to the cage for 10 min and then returned to the home cage. Then, the wood chips were spread evenly. Twenty glass marbles (4 × 5) were placed on top of the wood chips. Subsequently, each mouse was placed carefully back into the corner of cage, as far from marbles as possible. The mouse was allowed to explore for 10 min without being disturbed. Then, the number of marbles buried (50% or more covered) was counted. Before testing the next mouse, the wood chip spread was re-set, and the marbles were cleaned with 75% ethanol.

Cold stress defecation

Cold stress defecation experiment was performed in a refrigerator at 4 °C. First, the mouse was placed in a cage with a moist paper on the bottom and allowed to acclimate to the environment for 10 min. Then, the cage was placed in the refrigerator at 4 °C for 10 min. The amount of stool defecated during this period was quantified. Before the next mouse was tested, the moist paper in the cage was replaced.

Elevated plus maze

The EPM apparatus consisted of two open arms and two closed arms of equal dimensions (5 cm width × 35 cm length). Each arm was positioned at a 90° angle relative to the adjacent arms. The closed arms on the opposite sides were enclosed by 15 cm-high walls made of black Plexiglass, and the open arms were open without walls. The four arms were connected by a central area (5 cm width × 5 cm length). The entire EPM was lifted to 30 cm above the floor by the supports. Before the test, the mice were placed in the behavioral room to acclimate for 1 h. At the start of the test, each mouse was placed in the central area facing an open arm and allowed to freely explore the maze for 5 min. Before the next test, the EPM was wiped with a cleaning solution and dried with paper towels. The following variables were scored in the EPM: the percentage of time spent in the open arm (time spent on open arms/time spent in open and closed arms × 100) and the entries into open arms (entries into open arms/entries into open and closed arms × 100).

Total RNA extraction

The mice were sacrificed by cervical dislocation. Then, the brain was rapidly removed after decapitation. Subsequently, the mPFC region was dissected according to the atlas (Paxinos and Franklin, 4th Edition). The mPFC tissue was collected into the RNase-free tube and kept frozen at − 80 °C until analysis. Finally, 14 mPFC tissue samples were collected, including those from WT mice with loxP-flanked alleles in Efnb2 (WT; n = 4), Efnb2 gene knockout mice (KO; n = 4), WT mice after 1 h of the EPM test (WT_EPM; n = 3), and KO mice after 1 h of the EPM test (KO_EPM; n = 3). Total RNA was isolated and purified using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA amount and purity of each sample were assessed using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). RNA integrity was assessed by Bioanalyzer 2100 (Agilent, CA, USA) with RIN number > 7.0 and confirmed by electrophoresis with denaturing agarose gel.

RNA sequencing

The construction and sequencing of transcriptome libraries were performed at LC-BIO Bio-tech, Ltd. (Hangzhou, China) according to the protocol recommended by the vendor. Dynabeads Oligo (dT) 25-61005 (Thermo Fisher, CA, USA) was used to purify Poly (A) RNA from 1 μg total RNA. Then, the poly(A) RNA was fragmented into small pieces using Magnesium RNA Fragmentation Module (NEB, cat.e6150, USA) under 94 °C 5–7 min. Then, the cleaved RNA fragments were reverse-transcribed to create the cDNA by SuperScript™ II Reverse Transcriptase (Invitrogen, cat. 1896649, USA), which were next used to synthesize U-labeled second-stranded DNAs with Escherichia coli DNA polymerase I (NEB, cat.m0209, USA), RNase H (NEB, cat.m0297, USA) and dUTP Solution (Thermo Fisher, cat.R0133, USA). An A-base was added to the blunt ends of each strand, thereby preparing them for ligation to the indexed adapters. Each adapter contained a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single- or dual-index adapters were ligated to the fragments. Size selection was performed with AMPureXP beads. After the heat-labile UDG enzyme (NEB, cat.m0280, USA) treatment of the U-labeled second-stranded DNAs, the ligated products were amplified with PCR. The average insert size for the final cDNA library was 300 ± 50 bp. Finally, we performed the 2 × 150 bp paired-end sequencing (PE150) on an Illumina Novaseq™ 6000 (LC-Bio Technology Co., Ltd., Hangzhou, China).

Cutadapt software (version: cutadapt-1.9) was used to remove the reads that contained adaptor contamination. After removing the low-quality bases and undetermined bases, we used HISAT2 software (version: hisat2-2.0.4) to map reads to the genome. The mapped reads of each sample were assembled using StringTie (version: stringtie-1.3.4d. Linux_x86_64) with default parameters. Then, all transcriptomes from all samples were merged to reconstruct a comprehensive transcriptome using gffcompare software (version: gffcompare-0.9.8. Linux_x86_64). After the final transcriptome was generated, StringTie and ballgown were used to estimate the expression levels of all transcripts and to determine the expression level for mRNAs by calculating FPKM.

Weighted gene co-expression network analysis (WGCNA)

The WGCNA package was used to perform WGCNA in R studio based on FPKM expression data. The soft threshold (beta) value is equal to 14 which is calculated by the SoftThreshold() function of the WGCNA package. The process of automatic network construction and module detection were completed by the blockwiseModules() function of the WGCNA package. Specifically, Genes were clustered into branches of highly co-expression genes, and modules were identified with the tree cut algorithm.

The trait data contained four sets of binary variables, namely, WT, KO, WT_EPM, and KO_EPM, which were used to calculate the module trait relationships. The groups of interest were set to 1. The others were set to 0. Modules of interest were functionally annotated using annotation data package. The GO enrichment and the KEGG enrichment analysis were preformed using the clusterProfiler R package [32].

Ex vivo whole-cell patch-clamp recording

Whole-cell patch-clamp recordings were performed on WT and KO mice. The mice were first anesthetized by letting them inhale isoflurane and were sacrificed. The brains were rapidly removed and placed in the ice-cooled cutting solution (234 mM sucrose, 3.6 mM KCl, 1.2 mM MgCl2, 1.2 mM NaH2PO4, 12 mM glucose, 2.5 mM CaCl2, and 25 mM NaHCO3) for 90 s. The 300 µM coronal brain slices of PrL were cut on a vibroslicer (VT 1200S, Leica, Wetzlar, Germany) according to the mouse brain atlas. Then, the slices were immediately transferred into 32 °C artificial cerebrospinal fluid (ACSF, 124 NaCl, 2.5 KCl, 1.25 NaH2PO4, 1.3 MgSO4, 26 NaHCO3, 2 CaCl2, and 20 d-glucose, equilibrated with 95% O2 and 5% CO2) for 30 min and at room temperature for 1 h before recordings.

Electrodes (4–6 MΩ) were used for electrophysiological recording of neurons. The sEPSCs were recorded with cs-base interval solution (140 mM Cs-gluconate, 10 mM Hepes, 1.1 mM EGTA, 2 mM MgCl2, 3 mM MgATP, and 0.3 mM tris-guanosine triphosphate; pH 7.4 was achieved by adjustment with CsOH) with bath application of picrotoxin (PTX, 100 μM) and clamped at − 70 mV. The sIPSCs were recorded with high-Cl interval solution (140 mM CsCl, 1.1 mM EGTA, 10 mM Hepes, 2 mM MgCl2, 3 mM MgATP and 0.3 mM Tris-guanosine triphosphate; pH 7.4 was achieved by adjustment with CsOH) with application of 1 mM kynurenine and clamped at − 70 mV.

Signals were acquired by an Axon 700B amplifier (Axon Instruments, Foster City, CA) and fed into the computer through a Digidata-1550 interface (Axon Instruments) for data analysis (Clampfit 10.7). Access resistance (< 25 MΩ) was considered acceptable throughout all experiments. Moreover, the Signals in which resistance fluctuated in excess of 20% were excluded.

Statistical analysis

Statistical analyses were performed using Graph Pad Prism software (San Diego, CA, USA). In most cases, the unpaired two-tailed Student’s t test was used for normally distributed data. The Mann Whitney test was used for data not normally distributed. The D’Agostino-Pearson test was used for normality test. Data were given as means (± SEM), and differences in means were considered statistically significant at P ≤ 0.05. Significance levels [P < 0.05 (*), P < 0.01(**), P < 0.001(***) or P < 0.0001(#)] were indicated in the text and figures.

Availability of data and materials

The data supporting the conclusions of this article are included within its additional files.



Prefrontal cortex

PV+ :





Elevated plus maze


Weighted gene co-expression network analysis




Fragments Per Kilobase of transcript per Million mapped reads


Differentially expressed genes


Module Eigengene


Gene Ontology


Kyoto Encyclopedia of Genes and Genomes


Spontaneous excitatory post-synaptic currents


Spontaneous inhibitory post-synaptic currents


Medial prefrontal cortex


  1. 1.

    Osborn TL, Wasil AR, Weisz JR, Kleinman A, Ndetei DM. Where is the global in global mental health? A call for inclusive multicultural collaboration. Gen Psychiatr. 2020;33:e100351.

    Article  Google Scholar 

  2. 2.

    Ball TM, Ramsawh HJ, Campbell-Sills L, Paulus MP, Stein MB. Prefrontal dysfunction during emotion regulation in generalized anxiety and panic disorders. Psychol Med. 2013;43:1475–86.

    Article  Google Scholar 

  3. 3.

    Page CE, Coutellier L. Prefrontal excitatory/inhibitory balance in stress and emotional disorders: evidence for over-inhibition. Neurosci Biobehav Rev. 2019;105:39–51.

    CAS  Article  Google Scholar 

  4. 4.

    Kubota Y, Hattori R, Yui Y. Three distinct subpopulations of GABAergic neurons in rat frontal agranular cortex. Brain Res. 1994;649:159–73.

    CAS  Article  Google Scholar 

  5. 5.

    Ferguson BR, Gao WJ. PV interneurons: critical regulators of E/I balance for prefrontal cortex-dependent behavior and psychiatric disorders. Front Neural Circuits. 2018;12:37.

    Article  Google Scholar 

  6. 6.

    Pasquale EB. Eph-ephrin bidirectional signaling in physiology and disease. Cell. 2008;133:38–52.

    CAS  Article  Google Scholar 

  7. 7.

    Coulthard MG, Morgan M, Woodruff TM, Arumugam TV, Taylor SM, Carpenter TC, et al. Eph/ephrin signaling in injury and inflammation. Am J Pathol. 2012;181:1493–503.

    CAS  Article  Google Scholar 

  8. 8.

    Yamaguchi Y, Pasquale EB. Eph receptors in the adult brain. Curr Opin Neurobiol. 2004;14:288–96.

    CAS  Article  Google Scholar 

  9. 9.

    Dalva MB, McClelland AC, Kayser MS. Cell adhesion molecules: signalling functions at the synapse. Nat Rev Neurosci. 2007;8:206–20.

    CAS  Article  Google Scholar 

  10. 10.

    Talebian A, Britton R, Henkemeyer M. Abnormalities in cortical interneuron subtypes in ephrin-B mutant mice. Eur J Neurosci. 2018;48:1803–17.

    Article  Google Scholar 

  11. 11.

    Williams SE, Mann F, Erskine L, Sakurai T, Wei SN, Rossi DJ, Gale NW, Holt CE, Mason CA, Henkemeyer M. Ephrin-B2 and EphB1 mediate retinal axon divergence at the optic chiasm. Neuron. 2003;39:919–35.

    CAS  Article  Google Scholar 

  12. 12.

    Elowe S, Holland SJ, Kulkarni S, Pawson T. Downregulation of the ras-mitogen-activated protein kinase pathway by the EphB2 receptor tyrosine kinase is required for ephrin induced neurite retraction. Mol Cell Biol. 2001;21:7429–41.

    CAS  Article  Google Scholar 

  13. 13.

    Kania A, Klein R. Mechanisms of ephrin-Eph signalling in development, physiology and disease. Nat Rev Mol Cell Biol. 2016;17:240–56.

    CAS  Article  Google Scholar 

  14. 14.

    Henkemeyer M, Itkis OS, Ngo M, Hickmott PW, Ethell IM. Multiple EphB receptor tyrosine kinases shape dendritic spines in the hippocampus. J Cell Biol. 2003;163:1313–26.

    CAS  Article  Google Scholar 

  15. 15.

    Bouzioukh F, Wilkinson GA, Adelmann G, Frotscher M, Stein V, Klein R. Tyrosine phosphorylation sites in ephrinB2 are required for hippocampal long-term potentiation but not long-term depression. J Neurosci. 2007;27:11279–88.

    CAS  Article  Google Scholar 

  16. 16.

    Pellow S, Chopin P, File SE, Briley M. Validation of open: closed arm entries in an elevated plus maze as a measure of anxiety in the rat. J Neurosci Methods. 1985;14:149–67.

    CAS  Article  Google Scholar 

  17. 17.

    Pellow S, File SE. Anxiolytic and anxiogenic drug effects on exploratory activity in the elevated plus maze: a novel test of anxiety in the rat. Pharmacol Biochem Behav. 1986;24:525–9.

    CAS  Article  Google Scholar 

  18. 18.

    Degroot A, Wade M, Salhoff C, Davis RJ, Tzavara ET, Nomikos GG. Exposure to an elevated platform increases plasma corticosterone and hippocampal acetylcholine in the rat: reversal by chlordiazepoxide. Eur J Pharmacol. 2004;493:103–9.

    CAS  Article  Google Scholar 

  19. 19.

    Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9:559.

    Article  Google Scholar 

  20. 20.

    Kwon J, Kim YJ, Choi K, Seol S, Kang HJ. Identification of stress resilience module by weighted gene co-expression network analysis in Fkbp5-deficient mice. Mol Brain. 2019;12:99.

    CAS  Article  Google Scholar 

  21. 21.

    Okun E, Griffioen K, Barak B, Roberts NJ, Castro K, Pita MA, Cheng A, Mughal MR, Wan R, Ashery U, Mattson MP. Toll-like receptor 3 inhibits memory retention and constrains adult hippocampal neurogenesis. Proc Natl Acad Sci USA. 2010;107:15625–30.

    CAS  Article  Google Scholar 

  22. 22.

    Gu Z, Liu W, Wei J, Yan Z. Regulation of N-methyl-d-aspartic acid (NMDA) receptors by metabotropic glutamate receptor 7. J Biol Chem. 2012;287:10265–75.

    CAS  Article  Google Scholar 

  23. 23.

    Mcdonald AJ. Expression of the type 1 cannabinoid receptor (CB1R) in CCK-immunoreactive axon terminals in the basolateral amygdala of the rhesus monkey (Macaca mulatta). Neurosci Lett. 2020;745:135503.

    Article  Google Scholar 

  24. 24.

    Forrester SJ, Booz GW, Sigmund CD, Coffman TM, Kawai T, Rizzo V, Scalia R, Eguchi S. Angiotensin II signal transduction: an update on mechanisms of physiology and pathophysiology. Physiol Rev. 2018;98:1627–738.

    CAS  Article  Google Scholar 

  25. 25.

    Macova M, Pavel J, Saavedra JM. A peripherally administered, centrally acting angiotensin II AT(2) antagonist selectively increases brain AT(1) receptors and decreases brain tyrosine hydroxylase transcription, pituitary vasopressin and ACTH. Brain Res. 2009;1250:130–40.

    CAS  Article  Google Scholar 

  26. 26.

    Gong XX, Hu H, Qiao Y, Xu PF, Yang MQ, Dang RL, et al. The involvement of renin-angiotensin system in lipopolysaccharide-induced behavioral changes, neuroinflammation, and disturbed insulin signaling. Front Pharmacol. 2019;10:318.

    CAS  Article  Google Scholar 

  27. 27.

    Hawken NM, Zaika EI, Nakagawa T. Engineering defined membrane-embedded elements of AMPA receptor induces opposing gating modulation by cornichon 3 and stargazin. J Physiol. 2017;595:6517–39.

    CAS  Article  Google Scholar 

  28. 28.

    Santos SD, Carvalho AL, Caldeira MV, Duarte CB. Regulation of AMPA receptors and synaptic plasticity. Neuroscience. 2009;158:105–25.

    CAS  Article  Google Scholar 

  29. 29.

    Geng F, Tian J, Wu JL, Luo Y, Zou WJ, Peng C, Lu GF. Dorsomedial prefrontal cortex 5-HT6 receptors regulate anxiety-like behavior. Cogn Affect Behav Neurosci. 2018;18:58–67.

    Article  Google Scholar 

  30. 30.

    Gerety SS, Anderson DJ. Cardiovascular ephrinB2 function is essential for embryonic angiogenesis. Development. 2002;129:1397–410.

    CAS  Article  Google Scholar 

  31. 31.

    Liu J, Zhang MQ, Wu X, Lazarus M, Cherasse Y, Yuan MY, et al. Activation of parvalbumin neurons in the rostro-dorsal sector of the thalamic reticular nucleus promotes sensitivity to pain in mice. Neuroscience. 2017;366:113–23.

    CAS  Article  Google Scholar 

  32. 32.

    Yu GC, Wang LG, Han YY, He QY. Cluster Profiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7.

    CAS  Article  Google Scholar 

Download references


Not applicable.


This work was supported by the National Natural Science Foundation of China (Grant Numbers: 81501153, 81571326), Feixiang Personnel Training Program of Shanghai Mental Health Center (Grant Number: 2018-FX-02), Shanghai Key Laboratory of Psychotic Disorders (Grant Number: 13dz2260500).

Author information




YS designed the study, acquired data, interpreted data, and wrote the manuscript. LM acquired data, interpreted data, and wrote the manuscript. JC, SP, YC and YZ contributed to data acquisition. WW participated in data processing and answering questions during the revision process. JC and PJ designed this project and wrote the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Jinghong Chen or Peijun Ju.

Ethics declarations

Ethics approval and consent to participate

All experimental procedures and housing conditions for animals were approved by the National Institute of Health Guide for the Care and Use of Laboratory Animals and were approved by Animal Care and Use Committee in Shanghai Jiao Tong University (A2019102).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Summary of assembly statistics after Illumina sequencing.

Additional file 2.

DEGs in KO vs. WT, KO_EPM vs. WT_EPM, KO_EPM vs. KO and WT_EPM, vs. WT.

Additional file 3.

GO enrichment analysis for genes in DEG1 and DEG2.

Additional file 4.

Gene list of WGCNA.

Additional file 5.

List of genes within each WGCNA module.

Additional file 6.

GO enrichment analysis for genes in MEgreen module.

Additional file 7.

KEGG enrichment analysis for genes in MEgreen module.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Sun, Y., Ma, L., Chen, J. et al. RNA-seq co-expression network analysis reveals anxiolytic behavior of mice with Efnb2 knockout in parvalbumin+ neurons. Mol Brain 14, 118 (2021).

Download citation


  • Efnb2
  • Anxiety
  • RNA-seq