Transcriptome Analysis of the Human Tibial Nerve Identifies Sexually Dimorphic Expression of Genes Involved in Pain, Inflammation, and Neuro-Immunity

Sex differences in gene expression are important contributors to normal physiology and mechanisms of disease. This is increasingly apparent in understanding and potentially treating chronic pain where molecular mechanisms driving sex differences in neuronal plasticity are giving new insight into why certain chronic pain disorders preferentially affect women vs. men. Large transcriptomic resources are now available and can be used to mine for sex differences to gather insight from molecular profiles using donor cohorts. We performed in-depth analysis of 248 human tibial nerve (hTN) transcriptomes from the GTEx Consortium project to gain insight into sex-dependent gene expression in the peripheral nervous system (PNS). We discover 149 genes with sex differential gene expression. Many of the more abundant genes in men are associated with inflammation and appear to be primarily expressed by glia or immune cells, with some genes downstream of Notch signaling. In women, we find the differentially expressed transcription factor SP4 that is known to drive a regulatory program, and may impact sex differences in PNS physiology. Many of these 149 differentially expressed (DE) genes have some previous association with chronic pain but few of them have been explored thoroughly. Additionally, using clinical data in the GTEx database, we identify a subset of DE, sexually dimorphic genes in diseases associated with chronic pain: arthritis and Type II diabetes. Our work creates a unique resource that identifies sexually dimorphic gene expression in the human PNS with implications for discovery of sex-specific pain mechanisms.

Introduction

Sex-differential gene regulation and resultant changes in transcriptome, proteome, and metabolome shape sexually dimorphic physiology and behavior in animals. Sex-differential molecular profiles in human tissues have a profound effect on health, resulting in disease susceptibility, prevalence and pathophysiology differences between sexes. Acute and chronic pain have a staggering global disease burden, and prevalence of many chronic pain conditions like fibromyalgia and neuralgia have been shown to be higher in women (De Toledo et al., 2016), and sex-differential molecular changes in the peripheral and central nervous system (CNS) have been implicated in preclinical models (Mogil, 2012). Transcriptome profiles of human Dorsal Root Ganglia (DRG) and Trigeminal Ganglia have been characterized previously (Flegel et al., 2015; LaPaglia et al., 2018; Ray et al., 2018) but these studies are underpowered for capturing subtle transcriptional changes between sexes. Human tibial nerve (hTN) transcriptomes, which contains axons from DRG neurons, along with a panel of other harvested tissues, have been profiled using RNA-sequencing (RNA-seq) in hundreds of male and female donors as part of the GTEx project (Lonsdale et al., 2013). Some studies have characterized sex-differential gene expression changes (Chen et al., 2016; Gershoni and Pietrokovski, 2017) and investigated the evolutionary and regulatory basis of such changes across the repertoire of GTEx tissues, but none have focused on the hTN.

We cataloged sex differences in the hTN transcriptome under baseline conditions focusing on potential functional impact in the context of pain, inflammation and neuro-immunity. Previous analysis of the human peripheral nerve transcriptome has focused on changes in disease pathologies like nerve sheath tumors (Lee et al., 2014), or diabetic neuropathy (Hur et al., 2011; Luo et al., 2017). To our knowledge, sexual dimorphism studies on mammalian peripheral nerves is limited to a single microarray study in a rodent model of Type II diabetes (O'Brien et al., 2016). Our work thus fills an important gap in transcriptome studies of the human PNS and we created an online database (https://www.utdallas.edu/bbs/painneurosciencelab/sensoryomics/sexdiffnerve/) identifying sex differentially expressed genes in the hTN, and identified potentially functionally important genes by intersectional analysis with relevant databases. We further analyzed the cell type of expression of these genes, putatively identifying axonally localized transcripts. We also investigated the possible regulatory role of sex hormones in the sexual dimorphic genes by analyzing correlation of age with gene expression.

Materials and Methods

Data Use Policy

We only analyzed anonymized samples for which the corresponding donor consent information was available in the GTEx dataset (dbGAP phs000424.v7.p2) at the time of analysis. All of the datasets we analyzed were approved for General Research Use (GRU) and thus have no further limitations outside of those in the NIH model Data Use Certification Agreement. Our work was further approved by the UT Dallas IRB (protocol number 15-237).

Donor Selection by GTEx Consortium

The Genotype Tissue Expression (GTEx) consortium performed RNA-seq on a panel of harvested tissues, including the hTN, in a large cohort of consented tissue, organ or post-mortem donors. Donors of both sexes, from all ancestry groups between the ages of 21–70, and with Body Mass Index (BMI) between 18.5 and 35 were eligible. Post-mortem donors were also constrained to have had no whole blood transfusion within 48 h prior to death, and time between death and tissue collection was constrained to be <24 h. Donors were further restricted to have no history of metastatic cancer, no chemotherapy or radiation therapy within 2 years prior to death, and no communicable diseases that disqualify people for organ or tissue donation.

Sample Processing by GTEx Consortium

This provided a diverse dataset that is uniquely suitable for studying differences across sub-cohorts, such as sex differences in our study. Post-excision, aliquots were stabilized in a solution containing ethanol and methanol, acetic acid and Paxgene Tissue fixative (Qiagen). A portion of each tissue was then subject to RNA extraction and quantification. Next Generation sequencing (RNA-seq) was performed using the Illumina TruSeq platform, averaging ~50 million reads per sample. Sequenced reads were mapped to the human reference genome and transcriptome using the STAR toolkit (Dobin et al., 2013) and relative abundance quantified as Transcripts per Million (TPM) using the RSEM tool (Li and Dewey, 2011). RNA quality was checked using the RNA Integrity Number, and sequencing quality was analyzed using the RNA-SeqC toolkit (DeLuca et al., 2012). Donor level data included basic demographics, use of medication, medical history, results of laboratory tests, and circumstances of death were collected from the donor or next of kin and compared with the medical record. Donor level data and sample level data (which includes ischemic time, and comments from prosector and pathology reviewer) were analyzed for completeness and quality before release by the consortium. Further details of bench protocols and computational pipelines used by the GTEx consortium can be obtained from Lonsdale et al. (2013).

GTEx Data Requantification

PAXgene preserved hTN RNA-seq samples (dbGAP phs000424.v7.p2) with total RNA sequenced on the Illumina Truseq platform (and available donor consent information at the time of analysis) were identified. Samples noted to have sepsis, HIV infection, Type I diabetes or having both chronic joint pain and Type II diabetes were not used. The GTEx uniform processing pipeline provided relative abundance of genes in the form of normalized read counts as Transcripts per Million (TPM). GTEx RNA-seq assays used rRNA-depleted total RNA libraries containing reads from non-polyA transcripts with the proportion of such reads potentially varying between samples (Cui et al., 2010). We thus limited our analysis to validated coding genes by re-constraining the TPMs of coding genes [based on GENCODE annotation (Harrow et al., 2012)] to sum to a million.

Separation Into Cohorts

We mined the associated clinical information for samples to classify donors into three cohorts based on well-understood phenotypic changes in peripheral nerves: those noted to have arthritis or rheumatoid arthritis (chronic joint pain cohort, CJP), those noted to suffer from Type II diabetes (Type II diabetes cohort, T2D, which often causes diabetic neuropathic pain), and those without either of these diseases (baseline cohort, BSL) based on studies in the literature that suggested phenotypic and molecular profile changes in peripheral nerves in arthritic and diabetic patients (Hur et al., 2011; Pongratz and Straub, 2013). We then studied whether these cohorts could be analyzed together, or needed to be analyzed separately.

For the entire dataset including all three cohorts, only stably expressed genes (filtering out genes with median TPM across cohorts <0.5 or maximum TPM across cohorts <1.0) were analyzed, and each gene's normalized entropy score e_i (based on Ray et al., 2018) using samples from all cohorts was calculated as a measure of its variability:

$\begin{array}{l}{e}_i=-\frac{1}{lo{g}_2(N)}{\displaystyle \sum _j}\frac{t_{i,j}}{{\sum }_k{t}_i,k}lo{g}_2\frac{t_{i,j}}{{\sum }_k{t}_i,k}\mathrm{\text{.}}\end{array}$

where t_{i, j} is the TPM of the ith gene in the jth sample, with i indexing genes and j and k indexing samples, N being the total number of samples, and log 0 being defined as 0. Principal Components Analysis (PCA) was performed for samples (Figure 1A) using stably expressed, highly variable genes (with normalized entropy using samples in all three cohorts in the 90th percentile or above). The top two principal components were found to account for 72% of the total variance in the dataset but did not spatially segregate the samples on the basis of sex or membership in the three cohorts, which is expected since sex differences or disease pathologies like arthritis or diabetes are not expected to cause transcriptome-wide changes (like cancer). Based on the first two PCA dimensions, we performed outlier analysis to identify samples that were two standard deviations or further away from the origin (around which the principal component values are centered). We note an ~50% increase in the proportion of outlier samples in the CJP and T2D cohorts (both >28%) as compared to the BSL cohort (18.1%). Additionally, a preliminary analysis finds that over 250 genes are differentially expressed at the effect size of 1.5-fold change or above for both BSL–CJP and BSL–T2D comparisons (Figure 1B). These include genes well-known to be implicated in inflammation (for BSL–CJP comparison) and diabetes (for BSL–T2D comparison). These findings suggest very distinct phenotypes in the three cohorts, and any sex difference studies in the cohorts are thus more suited for analysis separately due to the large inter-cohort differences.

FIGURE 1

Figure 1. Population analysis of GTEx human tibial nerve samples. (A) PCA of the gene expression profile across all three cohorts of highly variable genes and subsequent plotting of the top 2 principal components show a larger fraction of outlier samples in CJP and T2D with respect to BSL. (B) The number of genes that are DE at 1.2-fold, 1.5-fold, and 2-fold when comparing the BSL with CJP and T2D are shown, along with some of the strongest DE genes. (C) The histogram of age in male and female BSL subcohorts show a similar distribution.

However, while rigorous statistical hypotheses testing to identify DE genes was performed on BSL (168:80 Male: Female Ratio/MFR) to characterize sex differences in healthy hTN (Supplementary Table 1), hypotheses testing was not performed in CJP (11:10 MFR) or T2D (44:16 MFR) due to these cohorts being underpowered for DE gene identification. Expression values, median fold change and Strictly Standardized Mean Differences for male–female gene expression changes in the two cohorts are shown in Supplementary Tables 2, 3 as a starting point for future studies. Strictly standardized mean difference (s_i) (Zhang, 2010) is a sample statistic that is often used to look at differences in effect size, while controlling for dispersion of data, under limited replication.

$\begin{array}{l}{s}_i=\frac{{\mu }_{F,i}-{\mu }_{M,i}+\in }{\sqrt{{\sigma }_{F,i}^2+{\sigma }_{M,i}^2}+\in }\end{array}$

where μ_{F, i} and μ_{M, i} are mean TPMs for female and male member of a cohort for the ith gene, σ_{F, i} and σ_{M, I} are the corresponding standard deviations with covariance assumed to be 0, with i indexing genes, and ε being a small smoothing factor (0.001) added to both the numerator and denominator.

Analysis of Potential Confounding Factors in Male and Female BSL Subcohorts

Based on donor level data, we performed a thorough analysis of documented variables that could potentially be a confounding factor in our analysis (results shown in Table 1). We identified several variables that were potentially confounding factors, including demographic (age, race), anatomical (BMI), and medical history based (presence of autoimmune disease like lupus or scleroderma, and inflammatory diseases like sarcoidosis or cellulites) variables. We find no appreciable difference in the breakup of BSL male and female subcohorts by race, and distribution of age (Figure 1C) and BMI are also comparable. Examined disease conditions were either completely absent in both subcohorts, or present in a small number of males (sarcoidosis). Height and weight distributions for the male and female subcohorts are noted to be different, as expected.

TABLE 1

Table 1. Statistics for BSL cohort.

DE Gene Identification

We performed non-parametric statistical hypotheses testing to identify DE genes. To minimize effects of multiple testing, we filtered out lowly expressed or undetectable genes (based on the conservative filtering criterion of median gene TPM <0.5 or maximum gene TPM <1.0 in both males and females in a cohort). We also filtered out genes that were ubiquitously expressed in BSL by calculating the normalized entropy (defined previously) across all samples in the BSL cohort. Since higher normalized entropy signifies more ubiquitous expression, we retained only those expressed genes in our analysis whose normalized entropy was less than the 75th percentile of normalized entropies of expressed coding genes, thus performing gene filtering in a manner agnostic to the sex of the samples. Wilcoxon rank-sum test (Wilcoxon, 1945) was used to calculate p-values for differences in the male and female sub-cohorts in BSL for the median (50th percentile) and the upper quartile (75th percentile), which can be robustly estimated given the cohort size. To test for differences at the 50th percentile (median), the entire male and female BSL sub-cohorts were used for comparison. In order to identify differences at the 75th percentile (upper quartile), only values greater than or equal to the medians of the male and female BSL sub-cohorts were used. To account for multiple testing, the Benjamini–Hochberg Procedure (BHP) (Benjamini and Hochberg, 1995) was used for both the 50th and 75th percentile tests with a False Discovery Rate (FDR) threshold of 0.05 for both tests (suggesting a combined FDR of 0.1). Since test statistics for the two tests are well-correlated due to overlap of data, the empirical FDR is effectively in the range (0.05, 0.1).

Age–TPM Correlational Analysis

In order to understand how age affects gene abundances on the sex differentially expressed gene set, we first performed PCA on the BSL samples using only the 149 sex-DE genes. Correlation [Pearsons' Correlational Coefficient R (Pearson, 1895)] was calculated between the principal components and age. We also calculated Pearson's R separately for the male and female subcohorts of BSL for the sex-DE genes. For the male subcohort, BHP was performed on the corresponding p-values to identify R > 0.20 or R < −0.20 as a statistically significant correlation. For the female subcohort, while the number of samples were less, and the number of comparisons also less (Y chromosomal genes not tested), we used the same thresholds for consistency.

The list of sex differentially expressed genes are present in Table 2, and a searchable public database with gene expression values for all genes, as well as all code and data for the analyses can be found online at the companion website https://www.utdallas.edu/bbs/painneurosciencelab/sensoryomics/sexdiffnerve/.

TABLE 2

Table 2. Statistically significant DE genes for BSL male vs. female subcohort analysis (statistically significant outcomes in boldface).

Results

We identified 29 genes in males, and 19 in females that were statistically significantly DE at ≥1.2-fold change (in their relative abundance/TPM) between the sexes in BSL at the 50th percentile. Additionally, we identified an additional 74 genes in males (Figure 2; Table 2), and 27 genes in females that were DE at the 75th percentile (Figure 3; Table 2). While several of these genes (including the 24 DE sex chromosomal genes) have been previously identified in the literature (Chen et al., 2016; Gershoni and Pietrokovski, 2017), our analysis identified new DE genes not identified in previous studies (Chen et al., 2016; Gershoni and Pietrokovski, 2017) including several that are potentially relevant to PNS function and pain mechanisms, including ISLR2, SP4, and TPPP2 (Chu et al., 2011; Aoki et al., 2014; Panza et al., 2015).

FIGURE 2

Figure 2. Differentially expressed genes in male tibial nerves. Quantile values for TPMs of sex-differentially expressed genes in BSL that are more abundant in males. Sex linkage, signal transduction and regulatory potential, and association with pain are also shown. Genes are shown in decreasing order of the area between the curves.

FIGURE 3

Figure 3. Differentially expressed genes in female tibial nerves. Quantile values for TPMs of sex-differentially expressed genes in BSL that are more abundant in females. Sex linkage, signal transduction and regulatory potential, and association with pain are also shown. Genes are shown in increasing order of the area between the curves.

Neural Tissue Enriched Genes

We found several genes that are sex-differential in expression and are enriched in expression in neural tissues. To identify neural tissue enriched genes, we identified genes in our DE gene list with high (≥0.5) Neural Proportion Scores (quantified in Ray et al., 2018) suggesting that they are primarily expressed in neurons or glia (in the CNS and PNS). NTRK1, which is known to be enriched in mammalian sensory neurons, was more abundant in males at both the 50th and 75th percentiles (Figure 2), suggesting sex-differential axonal mRNA trafficking. RNF165, involved in axonal growth, was also more abundant in male samples at the 75th percentile (Figure 2), and is another putative axonally transported mRNA. Well-known Schwann cell genes like Sonic Hedgehog (SHH) (Babcock et al., 2011) and Artemin (ARTN) (Lippoldt et al., 2013), that have been shown to be involved in pain, were also DE and higher in males at the 75th percentile (Figure 2). ISLR2 and MEGF11, known to be involved in axonal pathfinding, were differentially expressed in subsets of males and females, respectively (Figures 2, 3), and are potentially expressed in neurons and glia, respectively. Other neurally enriched genes that were DE in our datasets include PGBD5 and RNA helicase MOV10L1 in females [whose ortholog MOV10 has been shown to be implicated in nerve injury response in rodents (Melemedjian et al., 2011)], and AGAP2, NXPH3, and NRGN in males (Table 2). A small number of genes (SHH, EFCAB5, FABP7), which were neural tissue enriched, showed stronger or comparable sex-dependent expression in the hTN with respect to other tissues we examined (Supplementary Table 4).

To potentially identify the cell type of expression of neural tissue enriched genes, we analyzed the mouse nervous system single cell gene expression database at mousebrain.org (Zeisel et al., 2018) as well as the literature. Despite species' divergence of transcriptional regulation, we identified 15 genes whose murine orthologs were distinctly expressed in neuronal or glial cells in the mouse PNS (Figure 4A). Several of the neuronally expressed genes (NTRK1, ATP8A2, SRRM3, TMEM255A, RNF165, ISLR2) are expressed in the human DRG and have been shown in the literature to be important for axonal viability and routing or neurite outgrowth, suggesting that these are likely candidates for axonal mRNA transport and local translation.

FIGURE 4

Figure 4. Cellular enrichment and functional roles of differentially expressed genes. (A) For DE genes that are expressed in the human DRG, enriched in neural tissues, and detected in murine PNS single cell RNA-sequencing, median expression in TPMs in the human DRG, enrichment of expression in human neural tissues [based on the human neural proportion score (Ray et al., 2018)], and murine cell types of expression in the DRG (based on Zeisel et al., 2018) are shown. (B) LUM and WISP2 show similar sex-differential expression trends in BSL and CJP, and ARTN shows similar trends in BSL and T2D. (C) Protein Interaction Network for pro-inflammatory genes that are more abundant in the male cohort, suggesting a central role of NOTCH1. (D) Quantile plot for SP4 male and female sub-cohorts, and gene set enrichment analysis showing overlap of known SP4 targets from TRANSFAC (Matys et al., 2006) expressed in the hTN and present in Comparative Toxicogenomics Database (Davis et al., 2012) based pain associated gene list.

Pro-inflammatory Gene Signatures

The DE gene set with higher expression in subsets of male BSL sub-cohort showed a different pro-inflammatory gene signature (ARTN, HP, NOTCH1, CCL2, DIO2, and others), with respect to the female BSL sub-cohort (SULF1, GPR64, KRT4), suggesting sex-differential expression patterns of inflammatory gene markers under normal conditions (Table 2). This agrees with studies showing sex differences in clinical markers of inflammation (Casimir et al., 2010). Interestingly, the ratio of DE pro-inflammatory to anti-inflammatory genes in males were higher than in females, also in agreement with studies on human endothelial cells showing greater pro-inflammatory effects of androgens over estrogens (Annibalini et al., 2014). We analyzed the BSL male sub-cohort pro-inflammatory genes in T2D and CJP, and found only 3 inflammatory genes (LUM, WISP2, ARTN) with sex-dimorphic expression at comparable effect sizes in either T2D or CJP with respect to BSL (Figure 4B) suggesting that the pro-inflammatory gene signature is unlikely to be caused solely by a subset of donors in BSL that are suffering from inflammatory conditions unreported in the clinical record. The larger number of genes detected as more abundant in the male BSL sub-cohort can potentially be attributed to a combination of Y-chromosomal gene expression and downstream regulatory effects, and a larger set of DE pro-inflammatory genes.

We mined the set of pro-inflammatory protein that were more abundant in males for protein interaction based on both the literature and the StringDB database (Szklarczyk et al., 2016). We identified a Protein Interaction Network (PIN) connected through NOTCH1 (Figure 4C), known to be involved in neuropathic pain (Xie et al., 2015). This suggests that some of the pro-inflammatory genes more abundant in males (FABP7, FLT1, CCL2, PDGFRA, LUM, COMP) are potentially transcriptionally regulated by Notch signaling, or are co-regulated with NOTCH1. This is consistent with our knowledge of Notch signaling which is understood to drive a pro-inflammatory phenotype (Quillard and Charreau, 2013) in multiple tissues.

DE Regulatory Genes

Underlying causes for such sex-differential gene expression include Y chromosome gene expression, incomplete X inactivation in females (Gershoni and Pietrokovski, 2017), differential androgen and estrogen receptor driven regulation between sexes (Chen et al., 2016) and transcription regulatory programs controlled directly or indirectly by sex chromosomal gene products. Additionally, we find DE transcriptional regulators on autosomes including DMRTC1B and MED13 in males, and SP4 in females. SP4 is known to regulate pain-related genes, including several expressed in hTN (Figure 4D). While fold changes in these autosomal genes under baseline conditions are unlikely to affect transcription in their regulatory targets, the presence of sex-differential expression in a subset of samples suggests that under pathological conditions like inflammation and pain, more prominent sex-differential expression of these transcriptional regulators could potentially drive sex-specific regulatory programs.

Correlation of Gene Abundance With Age

We performed a principal components analysis of only the 149 sex differentially expressed genes, with the top two components clearly spatially separating male and female samples (Figure 5A). We also find that the top six principal components accounted for >90% of the variance (Figure 5B), but more interestingly four of the top six components are correlated with age (Figure 5C), both positively and negatively, showing that different parts of sex differential gene expression are both positively and negatively correlated with age. At the level of individual genes, we find that several genes have very different (weak) correlational patterns in males and females (Table 2). One possible explanation for this is regulation (direct or indirect) by sex hormones, especially via ESRA and ESRB, which have been studied in GTEx datasets (Chen et al., 2016). In fact, we find genes in our datasets (DDX3X, SULF1, TMEM125, and others more abundant in females; FLT1, MUSTN1, COMP, ARTN, REEP1, SEL1L2, and others more abundant in males) with negative correlation between female BSL gene abundance values and age, and a smaller correlative effect (or positive correlation) in male BSL samples. Several genes in our dataset like MUSTN1, REEP1, and FLT1 have been shown to be regulated by estrogen in the literature.

FIGURE 5

Figure 5. Association of differentially expressed genes with age. (A) PCA of the 149 sex-differentially expressed genes in BSL yields top two principal components that spatially segregate male and female samples. (B) The top six principal components explain >90% of the variance, and (C) 4 out of the top 6 principal components are weakly correlated with age of the samples.

Potential Protein Interaction Networks (PINs)

We investigated whether the gene products of the DE gene sets potentially interacted with genes whose expression is enriched in mammalian DRGs, which would help identify candidate sex-differential PINs in the PNS (Figures 6A,B). We used DRG-enriched genes from Ray et al. (2018) and our DE gene sets to identify putative PINs. The largest connected components from PINs generated using the StringDB database (Szklarczyk et al., 2016) show multiple DRG-enriched and male-preferentially expressed genes known to be expressed in glia (DUSP15, PRX, EGR2, DHH, FOXD3, ARTN), and involved in pain and inflammation, which points to a potential role for glia in sex differential pain processing in human peripheral nerves [shown in preclinical models (Averitt et al., 2018)]. Additionally, the presence of several neuronally expressed genes in the PIN among the set of DRG-enriched genes (GFRA3, NTRK1, NGFR, RET, and PPM1J) also suggests sex-differential glia-neuron crosstalk, which in turn can affect neuronal plasticity and excitability differently between the sexes.

FIGURE 6

Figure 6. Protein interaction networks of sex differentially expressed genes and human DRG enriched genes. Largest connected components of PINs based on known StringDB interactions between mammalian DRG-enriched genes (based on Ray et al., 2018), and DE genes that are more abundant in the (A) female BSL sub-cohort, and (B) male sub-cohort.

Roles in Pain and Neuro-Immunity

Out of 149 DE genes we identified, 56 (>33%) genes have known roles in pain, inflammation, innate immunity and other neuro-immune functions (Table 3). Sex differences in mammalian neuro-immune systems have been studied and linked to disease prevalence and susceptibility (Osborne et al., 2018), with PNS and CNS immune response sex differences studied in rodent pain models (Sorge et al., 2015; Mapplebeck et al., 2016; Lopes et al., 2017). In hTN, we find that DE genes involved in the innate immune system (C4A/B, CPAMD8 in males; DDX3X, TF in females) as well as genes known to be expressed in infiltrating or resident immune cells (ADA, PDGFRA in males; CD2, CD8B in females). We also found 34 genes that were associated with pain based on either the Comparative Toxicogenomics Database (Davis et al., 2012) or the Human Pain Genetics Database (Meloto et al., 2018), including genes that are central to pain pathways like NTRK1, and CCL2 in males and TMEM97 in females (Sahn et al., 2017; Ray et al., 2018).

TABLE 3

Table 3. Functionally relevant genes that are sex-differentially expressed in the human tibial nerve.

Discussion

Our work provides a biostatistical framework, and thoroughly catalogs sex-differential gene expression in hTN. The public resource (https://www.utdallas.edu/bbs/painneurosciencelab/sensoryomics/sexdiffnerve/) we generated provides a starting point for sex difference studies in human peripheral nerve drug target discovery, gene regulation and pathophysiology. We cataloged sex differences in gene expression in the PNS finding that many sexually dimorphic genes have known associations with pain, inflammation and neuro-immune interactions. This resource can be used for hypothesis-driven work to identify cell type and sub-compartment specificity of gene expression (e.g., NTRK1, LUM) by in situ hybridization or other techniques. These types of experiments will be important in identifying whether neuron-specific transcripts found in the tibial nerve are bona-fide axonally-localized transcripts or if they represent low level expression in other cell types that have not been recognized in single cell datasets to date. For genes like NTRK1 we favor the hypothesis that this is an axonally localized mRNA. From that perspective, it is likewise interesting that some of these genes are DE suggesting that sex hormones may drive mRNA localization in neurons. We are not aware of any previous studies that describe sex hormone-mediated control of axonal mRNA localization.

Our study can also be used to back-translate into transgenic preclinical models to identify potential sex differences in mechanisms of PNS physiology that may be relevant to pain therapeutics. While studies of molecular mechanism in human cohorts are difficult, sex differences in regulatory programs, signaling cascades and pathophysiological and behavioral changes can be studied in rodent-based perturbation models. Both SP4 and NOTCH1 are known to play key regulatory roles with respect to gene products that drive pain pathways (Chu et al., 2011; Xie et al., 2015), and are good candidates for further study.

While there are no rodent studies of baseline transcriptome sex differences in peripheral nerves, O'Brien et al. (2016) found sex differences in transcriptomes of sciatic nerves of diabetic mouse models show clear differences in the inflammasome. It is possible that the baseline differences in molecular profiles that we identified drive greater sexual dimorphisms under neuropathic conditions. While our study was underpowered for identifying sex-DE genes in the CJP and T2D cohorts, it can be expanded as GTEx cohort sizes continue to grow to potentially identify clinically relevant sex differences in CJP or T2D. Moreover, co-expressed (and putatively co-regulated) gene modules based on Whole Genome Correlation Network Analysis (Langfelder and Horvath, 2008) can be identified by finding correlated expression changes across cohorts and sexes. Given that sex differences in gene expression likely contribute to sexual dimorphism in neurological disease, such as chronic pain, exploiting these transcriptomic resources will be increasingly important for mechanism and drug discovery.

There are limitations to our work. As mentioned in the previous paragraph, our study was not able to identify sex-differential expression in Type II diabetes or arthritis, largely due to the small sample size. However, our work does identify some putative targets for future prospective studies. Most prominently, our findings suggest differential expression of immune-regulatory genes in tibial nerve of males and females. It may be reasonable to hypothesize that these differences will become more pronounced in chronic pain disorders. It is striking that we observed higher expression of some macrophage-expressed and/or macrophage recruiting genes (CCL2, NOTCH1) at higher levels in males given that many recent studies have identified macrophages and microglia as important contributors to pain specifically in male rodents (Mapplebeck et al., 2016; Lopes et al., 2017; Averitt et al., 2018; Inyang et al., 2018). We also observed higher expression of a well-known T cell gene in female TN, CD8B. This is interesting because T cells have been associated both with pain promotion (Sorge et al., 2015) and pain resolution (Krukowski et al., 2016) in female mice. Another shortcoming is that we have not addressed any specific hypotheses with interventional, prospective studies. This will be a future goal using rodent models. Finally, there is a potential confound of age in our analysis of sexual dimorphisms in gene expression. Our detailed analysis of this confound shows that age cannot explain all of the differences observed, and even when it does, it may be explained by changes in sex hormones across the lifespan.

In conclusion, our work demonstrates robust sexual dimorphisms in gene expression in the human tibial nerve. Some of these dimorphisms are likely important for understanding how the immune system interacts with the peripheral nervous system. Insofar as this is a major focus of current research into acute and chronic pain mechanisms (Mapplebeck et al., 2016; Lopes et al., 2017; Averitt et al., 2018; Inyang et al., 2018), these findings give new insight into translating findings from animal models into humans, as well as potential back-translation for testing of new hypotheses in animal models based on human molecular data.

Author Contributions

PR, AA, GD, and TP conceived the project. PR designed and performed experiments, analyzed and interpreted the data. JK performed exon-level analysis. DT-F performed pilot experiments. AW performed experiments for Supplementary Table 4. PR drew all figures. PR and TP wrote the paper, with input from GD.

Funding

NIH/NINDS grants R01NS065926 (to TP), R01NS102161 (to TP and AA) and R01NS098826 (to TP and GD).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

A preprint of this work is deposited in biorXiv at https://doi.org/10.1101/450197.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnmol.2019.00037/full#supplementary-material

Abbreviations

hTN, human tibial nerve; CJP, chronic joint pain cohort; T2D, Type II diabetes cohort; BSL, baseline cohort; DE, differentially expressed; BHP, Benjamini-Hochberg procedure;MFR, male:female ratio; PIN, protein interaction network; DRG, dorsal root ganglion; PNS, peripheral nervous system; PAR: pseudo-autosomal region; TPM, transcripts per million; BMI, body mass index; PCA, principal components analysis.

References

Annibalini, G., Agostini, D., Calcabrini, C., Martinelli, C., Colombo, E., Guescini, M., et al. (2014). Effects of sex hormones on inflammatory response in male and female vascular endothelial cells. J. Endocrinol. Invest. 37, 861–869. doi: 10.1007/s40618-014-0118-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Aoki, M., Segawa, H., Naito, M., and Okamoto, H. (2014). Identification of possible downstream genes required for the extension of peripheral axons in primary sensory neurons. Biochem. Biophys. Res. Commun. 445, 357–362. doi: 10.1016/j.bbrc.2014.01.193

PubMed Abstract | CrossRef Full Text | Google Scholar

Averitt, D. L., Eidson, L. N., Doyle, H. H., and Murphy, A. Z. (2018). Neuronal and glial factors contributing to sex differences in opioid modulation of pain. Neuropsychopharmacology 44, 155–165. doi: 10.1038/s41386-018-0127-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Babcock, D. T., Shi, S., Jo, J., Shaw, M., Gutstein, H. B., and Galko, M. J. (2011). Hedgehog signaling regulates nociceptive sensitization. Curr. Biol. 21, 1525–1533. doi: 10.1016/j.cub.2011.08.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B Methodol. 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

Casimir, G. J., Mulier, S., Hanssens, L., Zylberberg, K., and Duchateau, J. (2010). Gender differences in inflammatory markers in children. Shock 33, 258–262. doi: 10.1097/SHK.0b013e3181b2b36b

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C.-Y., Lopes-Ramos, C. M., Kuijjer, M. L., Paulson, J. N., Sonawane, A. R., Fagny, M., et al. (2016). Sexual dimorphism in gene expression and regulatory networks across human tissues. bioRxiv [preprint]. doi: 10.1101/082289

CrossRef Full Text | Google Scholar

Chu, C., Zavala, K., Fahimi, A., Lee, J., Xue, Q., Eilers, H., et al. (2011). Transcription factors Sp1 and Sp4 regulate TRPV1 gene expression in rat sensory neurons. Mol. Pain 7:44. doi: 10.1186/1744-8069-7-44

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, P., Lin, Q., Ding, F., Xin, C., Gong, W., Zhang, L., et al. (2010). A comparison between ribo-minus RNA-sequencing and polyA-selected RNA-sequencing. Genomics 96, 259–265. doi: 10.1016/j.ygeno.2010.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Davis, A. P., Murphy, C. G., Johnson, R., Lay, J. M., Lennon-Hopkins, K., Saraceni-Richards, C., et al. (2012). The comparative toxicogenomics database: update 2013. Nucleic Acids Res. 41, D1104–D1114. doi: 10.1093/nar/gks994

PubMed Abstract | CrossRef Full Text | Google Scholar

De Toledo, I. P., Réus, J. C., Fernandes, M., Porporatti, A. L., Peres, M. A., Takaschima, A., et al. (2016). Prevalence of trigeminal neuralgia: a systematic review. J. Am. Dental Assoc. 147, 570–576.e2. doi: 10.1016/j.adaj.2016.02.014

PubMed Abstract | CrossRef Full Text | Google Scholar

DeLuca, D. S., Levin, J. Z., Sivachenko, A., Fennell, T., Nazaire, M.-D., Williams, C., et al. (2012). RNA-SeQC: RNA-seq metrics for quality control and process optimization. Bioinformatics 28, 1530–1532. doi: 10.1093/bioinformatics/bts196

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Flegel, C., Schöbel, N., Altmuller, J., Becker, C., Tannapfel, A., Hatt, H., et al. (2015). RNA-seq analysis of human trigeminal and dorsal root ganglia with a focus on chemoreceptors. PLoS ONE 10:e0128951. doi: 10.1371/journal.pone.0128951

PubMed Abstract | CrossRef Full Text | Google Scholar

Gershoni, M., and Pietrokovski, S. (2017). The landscape of sex-differential transcriptome and its consequent selection in human adults. BMC Biol. 15:7. doi: 10.1186/s12915-017-0352-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrow, J., Frankish, A., Gonzalez, J. M., Tapanari, E., Diekhans, M., Kokocinski, F., et al. (2012). GENCODE: the reference human genome annotation for The ENCODE Project. Genome Res. 22, 1760–1774. doi: 10.1101/gr.135350.111

PubMed Abstract | CrossRef Full Text | Google Scholar

Hur, J., Sullivan, K. A., Pande, M., Hong, Y., Sima, A. A., Jagadish, H. V., et al. (2011). The identification of gene expression profiles associated with progression of human diabetic neuropathy. Brain 134, 3222–3235. doi: 10.1093/brain/awr228

PubMed Abstract | CrossRef Full Text | Google Scholar

Inyang, K. E., Szabo-Pardi, T., Wentworth, E., McDougal, T. A., Dussor, G., Burton, M. D., et al. (2018). The antidiabetic drug metformin prevents and reverses neuropathic pain and spinal cord microglial activation in male but not female mice. Pharmacol. Res. 139, 1–16. doi: 10.1016/j.phrs.2018.10.027

CrossRef Full Text | Google Scholar

Krukowski, K., Eijkelkamp, N., Laumet, G., Hack, C. E., Li, Y., Dougherty, P. M., et al. (2016). CD8⁺ T cells and endogenous IL-10 are required for resolution of chemotherapy-induced neuropathic pain. J. Neurosci. 36, 11074–11083. doi: 10.1523/JNEUROSCI.3708-15.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

LaPaglia, D. M., Sapio, M. R., Burbelo, P. D., Thierry-Mieg, J., Thierry-Mieg, D., Raithel, S. J., et al. (2018). RNA-Seq investigations of human post-mortem trigeminal ganglia. Cephalalgia 38, 912–932. doi: 10.1177/0333102417720216

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, W., Teckie, S., Wiesner, T., Ran, L., Granada, C. N. P., Lin, M., et al. (2014). PRC2 is recurrently inactivated through EED or SUZ12 loss in malignant peripheral nerve sheath tumors. Nat. Genet. 46:1227. doi: 10.1038/ng.3095

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., and Dewey, C. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12:323. doi: 10.1186/1471-2105-12-323

CrossRef Full Text | Google Scholar

Lippoldt, E. K., Elmes, R. R., McCoy, D. D., Knowlton, W. M., and McKemy, D. D. (2013). Artemin, a glial cell line-derived neurotrophic factor family member, induces TRPM8-dependent cold pain. J. Neurosci. 33, 12543–12552. doi: 10.1523/JNEUROSCI.5765-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Lonsdale, J., Thomas, J., Salvatore, M., Phillips, R., Lo, E., Shad, S., et al. (2013). The genotype-tissue expression (GTEx) project. Nat. Genet. 45:580. doi: 10.1038/ng.2653

CrossRef Full Text | Google Scholar

Lopes, D. M., Malek, N., Edye, M., Jager, S. B., McMurray, S., McMahon, S. B., et al. (2017). Sex differences in peripheral not central immune responses to pain-inducing injury. Sci. Rep. 7:16460. doi: 10.1038/s41598-017-16664-z

CrossRef Full Text | Google Scholar

Luo, L., Zhou, W. H., Cai, J.-J., Feng, M., Zhou, M., Hu, S.-P., et al. (2017). Gene expression profiling identifies downregulation of the neurotrophin-MAPK signaling pathway in female diabetic peripheral neuropathy patients. J. Diabetes Res. 2017:8103904. doi: 10.1155/2017/8103904

PubMed Abstract | CrossRef Full Text | Google Scholar

Mapplebeck, J. C., Beggs, S., and Salter, M. W. (2016). Sex differences in pain: a tale of two immune cells. Pain 157, S2–S6. doi: 10.1097/j.pain.0000000000000389

PubMed Abstract | CrossRef Full Text | Google Scholar

Matys, V., Kel-Margoulis, O. V., Fricke, E., Liebich, I., Land, S., Barre-Dirrie, A., et al. (2006). TRANSFAC^® and its module TRANSCompel^®: transcriptional gene regulation in eukaryotes. Nucleic Acids Res. 34, D108–D110. doi: 10.1093/nar/gkj143

CrossRef Full Text | Google Scholar

Melemedjian, O. K., Asiedu, M. N., Tillu, D. V., Sanoja, R., Yan, J., Lark, A., et al. (2011). Targeting adenosine monophosphate-activated protein kinase (AMPK) in preclinical models reveals a potential mechanism for the treatment of neuropathic pain. Mol. Pain 7:70. doi: 10.1186/1744-8069-7-70

PubMed Abstract | CrossRef Full Text | Google Scholar

Meloto, C. B., Benavides, R., Lichtenwalter, R. N., Wen, X., Tugarinov, N., Zorina-Lichtenwalter, K., et al. (2018). Human pain genetics database: a resource dedicated to human pain genetics research. Pain 159, 749–763. doi: 10.1097/j.pain.0000000000001135

PubMed Abstract | CrossRef Full Text | Google Scholar

Mogil, J. S. (2012). Sex differences in pain and pain inhibition: multiple explanations of a controversial phenomenon. Nat. Rev. Neurosci. 13:859. doi: 10.1038/nrn3360

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Brien, P. D., Hur, J., Robell, N. J., Hayes, J. M., Sakowski, S. A., and Feldman, E. L. (2016). Gender-specific differences in diabetic neuropathy in BTBR ob/ob mice. J. Diabetes Complicat. 30, 30–37. doi: 10.1016/j.jdiacomp.2015.09.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Osborne, B. F., Turano, A., and Schwarz, J. M. (2018). Sex differences in the neuroimmune system. Curr. Opin. Behav. Sci. 23, 118–123. doi: 10.1016/j.cobeha.2018.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Panza, P., Sitko, A. A., Maischein, H.-M., Koch, I., Flötenmeyer, M., Wright, G. J., et al. (2015). The LRR receptor Islr2 is required for retinal axon routing at the vertebrate optic chiasm. Neural Dev. 10, 23. doi: 10.1186/s13064-015-0050-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pearson, K. (1895). Note on regression and inheritance in the case of two parents. Proc. R. Soc. Lond. 58, 240–242. doi: 10.1098/rspl.1895.0041

CrossRef Full Text | Google Scholar

Pongratz, G., and Straub, R. H. (2013). Role of peripheral nerve fibres in acute and chronic inflammation in arthritis. Nat. Rev. Rheumatol. 9:117. doi: 10.1038/nrrheum.2012.181

PubMed Abstract | CrossRef Full Text | Google Scholar

Quillard, T., and Charreau, B. (2013). Impact of notch signaling on inflammatory responses in cardiovascular disorders. Int. J. Mol. Sci. 14, 6863–6888. doi: 10.3390/ijms14046863

PubMed Abstract | CrossRef Full Text | Google Scholar

Ray, P., Torck, A., Quigley, L., Wangzhou, A., Neiman, M., Rao, C., et al. (2018). Comparative transcriptome profiling of the human and mouse dorsal root ganglia: an RNA-seq–based resource for pain and sensory neuroscience research. Pain 159, 1325–1345. doi: 10.1097/j.pain.0000000000001217

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahn, J. J., Mejia, G. L., Ray, P. R., Martin, S. F., and Price, T. J. (2017). Sigma 2 receptor/Tmem97 agonists produce long lasting antineuropathic pain effects in mice. ACS Chem. Neurosci. 8, 1801–1811. doi: 10.1021/acschemneuro.7b00200

PubMed Abstract | CrossRef Full Text | Google Scholar

Sorge, R. E., Mapplebeck, J. C., Rosen, S., Beggs, S., Taves, S., Alexander, J. K., et al. (2015). Different immune cells mediate mechanical pain hypersensitivity in male and female mice. Nat. Neurosci. 18, 1081. doi: 10.1038/nn.4053

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2016). The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible. Nucleic acids research, gkw937.

PubMed Abstract | Google Scholar

Wilcoxon, F. (1945). Individual comparisons by ranking methods. Biometr. Bull. 1, 80–83. doi: 10.2307/3001968

CrossRef Full Text | Google Scholar

Xie, K., Qiao, F., Sun, Y., Wang, G., and Hou, L. (2015). Notch signaling activation is critical to the development of neuropathic pain. BMC Anesthesiol. 15:41. doi: 10.1186/s12871-015-0021-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeisel, A., Hochgerner, H., Lönnerberg, P., Johnsson, A., Memic, F., van der Zwan, J., et al. (2018). Molecular architecture of the mouse nervous system. Cell 174, 999–1014.e22. doi: 10.1016/j.cell.2018.06.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X. D. (2010). Strictly standardized mean difference, standardized mean difference and classical t-test for the comparison of two groups. Stat. Biopharm. Res. 2, 292–299. doi: 10.1198/sbr.2009.0074

CrossRef Full Text | Google Scholar