Sex-dependent differences in the genomic profile of lingual sensory neurons in naïve and tongue-tumor bearing mice

Mechanisms of sex-dependent orofacial pain are widely understudied. A significant gap in knowledge exists about comprehensive regulation of tissue-specific trigeminal sensory neurons in diseased state of both sexes. Using RNA sequencing of FACS sorted retro-labeled sensory neurons innervating tongue tissue, we determined changes in transcriptomic profiles in males and female mice under naïve as well as tongue-tumor bearing conditions Our data revealed the following interesting findings: (1) FACS sorting obtained higher number of neurons from female trigeminal ganglia (TG) compared to males; (2) Naïve female neurons innervating the tongue expressed immune cell markers such as Csf1R, C1qa and others, that weren’t expressed in males. This was validated by Immunohistochemistry. (3) Accordingly, immune cell markers such as Csf1 exclusively sensitized TRPV1 responses in female TG neurons. (4) Male neurons were more tightly regulated than female neurons upon tumor growth and very few differentially expressed genes (DEGs) overlapped between the sexes, (5) Male DEGs contained higher number of transcription factors whereas female DEGs contained higher number of enzymes, cytokines and chemokines. Collectively, this is the first study to characterize the effect of sex as well as of tongue-tumor on global gene expression, pathways and molecular function of tongue-innervating sensory neurons.

Sex-dependent differences in orofacial pain has been clinically well-established with higher prevalence in women in several different pathological conditions [1][2][3][4][5][6] ; although the mechanism of sex-differences remain elusive.It is known that sensory neurons including trigeminal neurons, that regulate pain in the orofacial region, are genomically different in males and females even under naïve conditions [7][8][9][10] .These studies have primarily been conducted using the entire neuronal population of the ganglionic tissue.However, several reports reveal that sensory innervation can be different with each tissue type 7,[11][12][13][14] .Accordingly, we have previously identified subsets of sensory neurons expressed in mouse tongue 13 that are varied from those innervating the masseter muscle 15 .Therefore, it is vital in delineating tissue-specific sex-dependent differences in trigeminal sensory neurons.Moreover, a significant gap-in knowledge exists for sex-specific changes of trigeminal neurons that specifically innervate diseased-tissues.Such studies can provide crucial information about the regulation of trigeminal sensory neurons in tissue-specific pathologies.
Therefore, the current study identified the differences in tongue-innervating sensory neurons between males and females under naïve and diseased state.We used tongue cancer as our disease model as approximately 50% of oral cancer patients report pain throughout the course of the disease and of these, the prevalence of pain is highest in tongue cancer patients [23][24][25][26]49 . Pai from tongue cancer is extremely weakening and significantly deteriorates patient quality of life in addition to having cancer due to very limited treatment options available.Therefore, using the orthotopic tongue cancer xenograft model in mice, we performed bulk-RNA sequencing of isolated tongue-innervated sensory neurons from males and females, to identify changes in genes, biological processes and molecular function between sexes under normal and tumor-bearing conditions.

Results
Isolation of lingual TG neurons revealed differences in number of innervating neurons between males and females.To assess transcriptomic changes in tongue-innervating TG neurons in males and females under naïve and tumor-bearing conditions, we performed bulk-RNA sequencing of lingual TG neurons.
To achieve this, we first isolated retro-labeled TG neurons innervating mouse tongue using WGA-488 and subsequently flow-sorted fluorescently labeled neurons from four groups: Male normal (MN), Male tumor (MT), Female Normal (FN) and Female Tumor (FT).Representative gating strategy for the sorting protocol is shown in Fig. 1A.Flow sorting obtained an average of approximately 10 K to 23 K neurons per sample for all four groups (Fig. 1B and Suppl Table 1).Surprisingly, we found that the number of neurons isolated from all of female samples from both groups (i.e.FN and FT,) were higher (approx.22,000 cells) compared to male groups (i.e.MN and MT, approx.12,000 cells) (Fig. 1B).In potting the percentage of WGA + neurons sorted over all live events in the samples, we found that percentage of WGA + neurons isolated from female samples were significantly higher (approx.threefold, one-way ANOVA, p < 0.05) than in males (Fig. 1C) suggesting that female mice may have increased number of sensory neurons innervating the tongue tissue compared to males in mice.To further investigate this finding, we employed immunohistochemistry to evaluate nociceptive and non-nociceptive WGA + neurons in naïve male and female mice.We used TRPV1 to distinguish between the two neuronal classes and found that males had a higher percentage of TRPV1 + tongue innervating neurons compared to females Fig. 1D and E) (25.3% males vs. 18.89% females, two-way ANOVA, p = 0.0197).Accordingly, females had a higher percentage of TRPV1 negative tongue innervating neurons than males (Fig. 1D and E) (74.39% males vs. 81.10%females, two-way ANOVA, p = 0.0168).
Sex-dependent differences in gene expression in lingual neurons of naïve mice.Bulk-RNA sequencing was conducted of flow-sorted WGA + TG neurons from all four groups: MN, MT, FN and FT.Data were analyzed to identify differentially expressed genes (DEGs) using the criteria: RPKM > 5, FC > 1.5 and p < 0.05.Our RNA sequencing data of markers of neuronal and non-neuronal cells confirmed that our samples from all four groups were enriched in neurons post-flow sorting (Supplementary Fig. 1).In comparing FN versus MN groups, we identified a total of 81 DEGs between both sexes after excluding all sex-linked genes.Of these, 30 genes were exclusively expressed in females and only 2 genes exclusively expressed in males (Fig. 2A and B.).The remaining 49 genes, while expressed in both sexes, were significantly upregulated in females compared to males (Fig. 2C and D).Details of RPKM and FC of the top 10 genes is listed in Supplementary Table 2. Two genes, that were higher in males than females were Fam23a (transmembrane protein 236) and Ddx3y (DEAD box helicase).In further assessing the expression of female-specific genes, we tested the expression of selected two genes: Csf1R and C1qa using immunohistochemistry and confirmed their higher expression in females than in males.Csf1R was expressed in ~ 42% of all WGA + neurons in females compared to 5.5% of WGA + neurons in males (Unpaired Student T test, p = 0.041) (Fig. 2E).Of all Csf1R positive neurons innervating the tongue in females, 18% were found to be TRPV1 positive (Fig. 2F and G) (Paired Student's T test, p = 0.035).Similarly, ~ 52% of all tongue innervating neurons expressed C1qa in females whereas almost no expression of this gene was found in males (Fig. 2H) (Unpaired Student T test, p = 0.025).Majority of C1qa positive neurons in the tongue of females were TRPV1 negative (94%) with only a very small proportion of C1qa expression in TRPV1 + nociceptors (6%) (Fig. 2I and J) (Paired Student's T Test, p = 0.0005).
To assess whether the genes exclusively expressed in females may have a role in nociception, we tested the function of Csf1R in sensitizing nociceptors due to its expression in TRPV1 + neurons.Accordingly, we explored whether its ligand Csf1 sex-selectively increases TRPV1 responses in TG neurons.As shown in Fig. 3A, proportion of WGA + and WGA-neurons that were responsive to CAP did not change with Csf1 treatment in females.However, Csf1 significantly increased CAP-evoked calcium accumulation of WGA + as well as WGAcells (Fig. 3B and C).Expectedly, Csf1 treatment neither altered the percentage of CAP responsive neurons (Fig. 3D) nor CAP-evoked calcium accumulation in male neurons (Fig. 3E and F).These data indicate that Csf1R in sensory neurons may sex-specifically contribute to nociceptive responses.
Between DEGs observed from MT versus MN and FT versus FN, we found 18 genes that were common to both comparisons whereas majority of the genes were regulated in a sex-specific manner (Fig. 4G).Interestingly though, out of the 18 common DEGs, 11 were upregulated in both comparisons whereas 3 genes that were upregulated in MT versus MN were significantly downregulated in FT versus FN (Suppl Table 3).In contrast 4 genes that were downregulated in MT versus MN were upregulated in FT versus FN (Suppl Table 3).www.nature.com/scientificreports/Collectively, these data suggests that the regulation of lingual neurons by the tongue tumor may be considerably different in males and females.

Biological processes and molecular function regulated upon tumor growth in males and females.
We next conducted gene ontology analyses to identify key biological processes in males and females upon tumor growth.Fourteen distinct biological processes were identified with upregulated DEGs in males upon tumor growth (Fig. 5A).These processes pertained to signaling mechanisms, immune process and inflammation, metabolic processes as well as some others including response to glucocorticoid stimulus, circadian regulation glial cell proliferation and apoptosis.Majority of DEGs in males were involved in signaling pathways and immune and inflammatory processes with predominant processes being second messenger signaling, response to growth factors, inflammatory response or leucocyte migration (Fig. 5A and Suppl Table 4).Additionally, considerable number of DEGs involved in cell differentiation and cell death were found to be upregulated in males upon tumor growth.Similar to males, sixteen processes were identified in females (Fig. 5B), although all were classified under regulation of signaling pathways, immune and inflammatory processes as well as other processes such as cell death, transport regulation and positive regulation of angiogenesis.DEGs contributing to metabolic processes and cell proliferation and differentiation were not found in females, unlike in males (Fig. 5B).Five specific processes were found to be common in both sexes.These included response to interferon-gamma, IL-1 signaling, leucocyte migration, inflammatory response and apoptosis (Fig. 5A and B).However, the number and type of DEGs involved in each of these pathways were different between males and females (Suppl Table 4 and  5).While no specific processes were identified for downregulated DEGs in males due to very few numbers of genes, downregulated DEGs in females were found to be associated with processes such as negative regulation of angiogenesis, locomotion, cell adhesion and cytoskeletal organization (Fig. 5C).List of downregulated DEGs associated with each of the BPs in females in given in Suppl.Table 6.Additionally, we also identified specific transcription factors (TFs), ligands, peptides, growth factors, receptors, channels, enzymes, chemokines and cytokines that were differentially regulated in males and females upon tumor growth.Of all DEGs between both comparisons, we found 11 TFs, 19 ligands, peptides and growth factors, 24 channels and receptors, 65 enzymes and 14 chemokines/cytokines (Fig. 5D).While the number of TFs were higher in male neurons upon tumor growth, female neurons had higher number of channels and receptors, enzymes and chemokines and cytokines, regulated post-tumor growth (Fig. 5E and F).Taken together, these data indicate significant differences in regulation of genes, processes and pathways by tongue tumor, between males and females.

Discussion
It is well established that chronic orofacial pain is sexually dimorphic with a higher prevalence in women than men for various pain conditions including temporomandibular joint disorders 50 , apical periodontitis 4 , oral mucositis 32,33,43 , burning mouth syndrome 31 as well as oral cancer 42,46,47 .Yet, there is a large gap in knowledge about the mechanisms that lead to sex-dependent differences in orofacial pain.Moreover, sex-dependent global genomic changes in trigeminal neurons during disease is entirely unexplored.One prior study that reported changes in transcriptomic profiles of trigeminal ganglia in conditions of masseter muscle inflammation in rats, only included males 51 .On the other hand, another study that investigated changes in gene expression of TG following neuropathic pain in males and females, did not evaluate the whole genome 52 .Furthermore, both these studies employed whole ganglionic tissues that primarily represent non-neuronal population as recently confirmed by Mecklenberg et al. 7 .To address this drawback, a third study employed single-cell sequencing of human and mouse TG tissues to implicate genes and cell types in migraine.However, sequencing of TG tissue in this study only evaluated sex-differences under naïve conditions and genomic data obtained from mouse migraine models was not separated by sex 10 .Therefore, no study till date has investigated changes in transcriptomic profiles of the trigeminal sensory neurons during disease in both sexes.Furthermore, it is critical to study  www.nature.com/scientificreports/neurons-innervating specific tissues over all ganglionic neurons as it is indicated that type of sensory innervation is tissue-specific.For example, it has been reported that dental pulp innervating neuronal soma primarily are large diameter myelinated neurons unlike other tissues in the orofacial region such as the skin or mucosa 53,54 .Accordingly, percentage of TRPA1 expressing neurons are higher in the oral mucosal tissues than in dental pulp 54 .Similarly, we have shown that about 20% of tongue innervating sensory neurons are CGRP positive 13 whereas almost 50% of masseter muscle innervating neurons are CGRP positive 15 .Therefore, in the current study, we explored changes in genomic profile of isolated sensory neurons innervating mouse tongue tissue and identified sex-dependent differences.We employed the approach of bulk-RNA sequencing of FACS sorted retro-labeled mouse tongue innervating TG neurons, to be able to concentrate neuronal population from whole ganglionic tissue as well as explore tissue-specific differences in neurons between sexes.Additionally, we utilized the preclinical mouse tongue cancer model as our disease model, as it has been shown that patients with tongue cancer are more commonly in pain than other oral cavity cancers and it is widely included in clinical and preclinical studies of oral cancer pain [35][36][37][38][39][40] .However, till date, there is only a couple of studies investigating sexually dimorphic mechanism for tongue cancer pain 42,46 .While these studies demonstrated the role of the immune response in sex-dependent differences in oral cancer pain, differential gene expression in sensory neurons upon tumor growth cannot be excluded, as we and others have reported that sensory neuronal activity is altered upon oral tumor growth 27,36,37,[55][56][57] .
FACS sorting of lingual TG neurons was performed from normal and tumor-bearing animals for both sexes and our data indicated that female mice may have higher number of tongue-innervating sensory neurons than males as observed by significant increased percentages of sorted cells from female samples than male samples.In further exploring this result, we found that female mice had higher percentage of TRPV1-and lower percentage of TRPV1 + neurons innervating the tongue than in male mice.However, recently Scheff et al., demonstrated no difference in capsaicin sensitivity in the tongue between males and females under naïve conditions 56 , indicating that while the percentage of TRPV1 + neurons were found to be different between the two sexes, the total number of TRPV1 + lingual neurons may not be different between males and females.Therefore, given that the total number of tongue-innervating neurons are higher in females, our data suggests that perhaps the total number of non-nociceptive sensory neurons in mouse tongue may be higher in females than in males.Whether or not www.nature.com/scientificreports/ the increased innervation of non-nociceptive neurons contributes to sex-dependent pain observed in any of the chronic lingual diseased states is yet to be determined.To our knowledge, this is the very first evidence of differences in number of neurons innervating orofacial tissues between sexes.Further validation of this finding using different tracer dyes, strain of mice and alternative approaches such as whole tissue imaging of retro-labeled TG tissues, is warranted.Besides, whether this finding is specific to mice or exists in tongue tissues of higher order species such as human and non-human primates is unknown.
In analyzing RNA sequencing data between naïve male and female neurons, we found that several genes (i.e. 30 genes) were exclusively expressed in females and not in males.This was in accordance to previous report showing increased DEGs in naïve female versus male neurons from whole mouse trigeminal ganglia 7 .Interestingly, many of these DEGs such as CSF1R (colony-stimulating factor receptor 1) 58,59 , C1qa (complement component 1 q) [60][61][62] , Sh2b3 (lymphocyte adapter protein) 63 ,Hhex (Hematopoetically expressed homeobox) 64,65 and Retnlg (Resistin-like gamma) 66 , are known to be primarily expressed in immune cells and have shown to play various roles in immune processes 58,62,67 , myeloid cell differentiation 64,65 , cytokine signaling 63 and inflammation 66,68 .In fact, some of these genes such as Csf1R and C1qa are considered specific markers of mononuclear phagocytic system such as microglia, monocytes, macrophages and dendritic cells 58,59,61,62 .Csf1R in microglia and macrophages has been reported to be activated in neuropathic pain via its ligand Csf1 expressed in immune cells or in sensory neurons 67 .Additionally, it has been shown that Csf1R exerts its action in a sex-specific manner by favoring a response in males compared to females 69 .However, for the first time, we report the expression of Csf1R in trigeminal sensory neurons, although its expression is specific to female neurons.Our immunohistochemical analyses of Csf1r and C1qa confirmed our RNA sequencing result and revealed that while majority of C1qa was expressed in non-TRPV1 expressing neurons, Csf1R was expressed in nociceptive and non-nociceptive neurons.Moreover, we also found that Csf1R signaling by its ligand Csf1 sensitizes TRPV1 responses in female sensory neurons but not in males indicating that at least some of the DEGs exclusively expressed in females may contribute to nociception.Interestingly, this effect was observed in WGA + and WGA-cells suggesting the role of Csf1/ Csf1R signaling in females may not be restricted to the tongue tissue.Similar to Csf1R, the contribution of C1q proteins in pain including orofacial pain as been reported to be in microglia [70][71][72] .Therefore, the sensory neuronal role of C1q proteins in pain is entirely unexplored possibly due to the fact that most studies are conducted in males that do not express neuronal C1q.Our data show that unlike Csf1R, C1q in females, is primarily expressed in TRPV1-neurons and not in TRPV1 + neurons.Therefore, further experiments are needed to fully characterize its expression in TRPV1-nociceptor and non-nociceptor subtypes, that will provide insights on the whether C1q may be involved in peripheral, secondary or central sensitization during tongue cancer or even other orofacial conditions.Nonetheless, our data from normal male and female animals indicate that perhaps female sensory neurons may be functionally different than male neurons innervating the mouse tongue.
Our analyses for gene expression changes post tongue-tumor growth revealed altered expression of several DEGs in both sexes.However, the number of DEGs in females was considerably higher than in males (83 DEGs in males vs. 382 DEGs in females) suggesting that male neurons are more tightly regulated in spite of the tumor compared to female neurons.It would be interesting to investigate whether this result is specific to tongue tumors or even other lingual pathologies.Notably, while we observed similar pain behaviors (as measured by previously reported feeding behavior and von Frey thresholds in the virbissal pad 27 ) and tumor growth at the selected time point between both sexes (Fig S2A , B and C), there were only 18 genes that were commonly altered between males and females, indicating that the mechanisms of neuronal regulation in oral cancer may be considerably different between sexes, even if pain manifestation is similar.The concept of differences in sex-specific mechanisms with similar pain behaviors is not novel and has been described in neuropathic pain models by others.For example, it has been reported that pain behavior responses are not different between sexes post spared nerve injury 73,74 ; however processing of pain through either microglial function 74 or the gene expression profile in sensory neurons is sex-specific 75 .Clinically, association of gender differences in oral cancer pain has been varied.While, several studies report no association of gender and oral cancer pain [76][77][78][79][80][81][82][83][84][85] , one study reported that initial presence of pain may be more common in men 35 and another study reported higher levels of function-related pain in men than women 47 .However, three other studies report worse pain scores for women than men 86,87 .Pre-clinically, differences in pain behaviors appears to be varied as well and may depend on the tumor model used.For example, Scheff et al. 46 , showed that females had worse gnawing behavior in the 4NQO tongue cancer model compared to males whereas we observed no difference in pain behaviors in the current study with HSC3 tongue tumors (Fig S1).Nonetheless, our data has clinical implications as it is crucial to understand the differences in mechanisms of oral cancer pain between males and females, to be able to understand if there may be potential differences in treatment responses between the sexes to ultimately be able to develop effective therapeutic strategies for both sexes.
Since the tumor developed in both sexes was from the same cell line (i.e.HSC3), perhaps the differential regulation of neurons between males and females may be due to stark differences in the tumor microenvironment or alterations occurring at the ganglionic level upon tumor growth.Intraganglionic activation of macrophages 88 and satellite glial cells 89,90 has been reported in pain.While no studies have yet probed the contribution of these cells for oral cancer pain, a similar response within the TG can be expected that in turn would lead to changes in the neuronal soma.On the other hand, the peripheral tumor microenvironment consists of various cell types including keratinocytes, fibroblasts, endothelial cells, Schwann cells, and immune cells.Schwann cells 91 and immune cells 42,46 have already been shown to contribute to tongue cancer pain and as mentioned above, immune cells even play a sex-dependent role in tongue cancer pain 42,46 .Therefore, examining the role of these cell types in mediating global changes in the neurons of males and females would be crucial in better understanding of the impact of the tumor microenvironment in sexually dimorphic tongue cancer pain.
It is noteworthy though that despite only 18 DEGS common to both males and females, 5 biological processes (BPs) were found to be common between males and females post tumor growth.These included apoptosis, response to interferon gamma, inflammatory response, leucocyte migration and IL-1 signaling.Of the 18 www.nature.com/scientificreports/common DEGs, four genes were associated with two of the common BPs; i.e. inflammatory process and apoptosis.These four genes were Timp 1(tissue inhibitor of metallopeptidase 1) Gal (galanin), Chi3l3 (chitinase-like 3) and Cxcl10 (chemokine-ligand 10).Timp1 is a glycoprotein that is known to promote cell proliferation and anti-apoptosis; is implicated in cancer progression 92 ; and negatively regulates matrix metalloproteinases and disintegrin-metalloproteinases (ADAMs) 93 .Interestingly, ADAM17 is implicated in oral cancer pain 94 .Accordingly, Timp1 has been shown to attenuate inflammatory pain in preclinical models 95 .Galanin is a neuropeptide and is considered a potent modulator of inflammation by promoting cytokine production in immune cells 96,97 .It also has been shown to induce cell death in pheochromocytoma cells 98 .Importantly, the role of galanin in pain has been reported to have pro-and anti-nociceptive functions and whether or not galanin plays a role in peripheral nociceptive mechanisms is yet to be confirmed 99 .Interestingly, galanin release from sensory neurons have been shown to promote oral cancer progression 100 .
Chitinase-like proteins belong to the family of glycoside hydrolase and are involved in the regulation of the innate immune response 101 .Chi3l3 has been specifically reported to orchestrate recruitment of eosinophils in meningitis and autoimmune neuroinflammation 102,103 .However, its role in pain is not yet defined.
Cxcl10 is an important chemokine for inflammatory processes and its function in pain has been studied in several pain models including neuropathic pain and inflammatory pain [104][105][106] .
Aside from the above-mentioned four genes, many other genes were selectively altered in each sex, yet associated with the common BPs.For example, the cytokine, interleukin-1beta (IL1β) was specifically shown to be induced in female post tumor growth and was not expressed in males.The role of IL1β has been widely reported in pain, with the cytokine mostly produced by immune cells and other non-neuronal cells during injury 107,108 .However, because most of the studies have been conducted using male animals, neuronal IL1β has not been reported previously.To this end, a sex-specific role of IL1β in pain has not been studied till date.Interestingly, while induction of neuronal IL1β was only observed in females, downstream signaling of IL1 pathway was observed in both sexes post tumor growth indicating that perhaps IL1 β in males may be increased in the periphery or within the ganglia in non-neuronal cells as demonstrated in other injury models [107][108][109][110] .
Another example was expression of Hspa1a which encodes for heat-shock protein 70 (Hsp70).Our RNA sequencing data showed that while this gene was expressed in both sexes, it was specifically upregulated only in males post tumor growth.Interestingly, hsp70 has been demonstrated to have a protective role in pain during nerve damage 111,112 , migraine 113 and opioid-induced hyperalgesia 114 .Therefore, upregulation of Hsp70 in males post-tumor growth might indicate that male neurons might express an endogenous feedback mechanism to suppress pain that may be lacking in females.
In addition to the common BPs, sex-selective BPs were identified post-tumor growth.Noteworthy femaleselective BPs including cytokine/chemokine production and signaling as well as angiogenesis.The contribution of various cytokine/chemokines in nociception and inflammation [115][116][117] is widely established and our data showed that the number of cytokine/chemokines altered in female neurons was higher than in males (i.e. 14 in females versus 3 in males).This is an intriguing yet not a surprising observation as it aligns with the finding that female neurons express immune cell markers as described above.
In contrast, male-specific BPs included growth factor signaling and ion transport.One of the top genes upregulated in males upon tumor growth; fibroblast growth factor 23 (Fgf23) is reported to be associated to both of these processes.It is not only a growth factor that induces downstream signaling via its receptors 118 but also controls phosphate homeostasis 119 .Furthermore, it is also considered a bone-derived hormone that is regulated by inflammation and associated with bone pain 120,121 .
Taken together, our data points to significant differences in the regulation of lingual sensory neurons in males and females upon tongue tumor growth.The current study is significant as it is the first to comprehensively characterize the genomic profile of tongue-innervating neurons under naïve and tumor-bearing conditions.Our data lays the foundation for future investigations to identify potential sex-specific targets in sensory neurons and study their functional and mechanistic contribution in tumor-induced pain.Of interest are the DEGs belonging to the class of transcription factors in males and enzymes and cytokine/chemokines in females as the total number of DEGs in these classes were higher in the respective sexes.One limitation of the study include use of one cell-line to induce the tongue tumor and at this point, the study cannot confirm which DEGs and processes are common across tumors induced by different OSCC cell lines.We also note that the study was conducted using the xenograft model with athymic mice that lack the response of T lymphocytes.While this model allows to explore the regulation of sensory neurons by human tumors, it excludes the involvement of T cells in tumor-induced pain.Therefore, identifying DEGs in both sexes that are common to the xenograft as well as the syngeneic model will be important in selecting relevant targets for functional studies.Besides, it would be useful to delineate the temporal effect of the tumor on neuronal genes and processes to further gain insight into sex-specific regulation of sensory neurons in oral cancer.

Animals.
Six-to eight-week-old adult Balb/c male and female athymic nude mice (Jackson Labs, Bar Harbor, ME, USA) were used for all experiments.All animal handling and procedures were performed according to approved UTHSCSA IACUC protocols and conformed to the guidelines of International Association for the Study of Pain (IASP).The study is reported in accordance to the ARRIVE guidelines.Animals were housed in the UTHSCSA laboratory of Animal Resources (LAR) for at least 4 days prior to start of experiments.line: HSC3, tongue tumors were induced in mice as described by us previously 27,36,37 .HSC3 cells were purchased from Health Science Research Resources Bank, Japan.Animals were anesthetized with isoflurane inhalation and 50 ul of 3.5 × 10^5 HSC3 cells were injected unilaterally in the ventral side of the tongues using insulin syringes.Normal group similarly received 50 ul of 3.5 × 10^5 human normal oral keratinocyte (NOK) cell line; OKF6-TERT2 as described previously 27,36,37 .NOK cell line was kindly provided by Dr. Cara Gonzales at UTHSCSA.Animals were then allowed to recover in their cages.Animals were used for experiments at day 15 post-cell inoculation.
Retro-labeling of tongue-innervating TG neurons.Tongue-innervated sensory neurons were labeled as described by us previously 13 .Briefly, animals were anesthetized with isoflurane inhalation and 10 uls of 1% wheat germ agglutinin (WGA)-AF488 (Promega), diluted in 1%DMSO, was injected bilaterally twice in tongue of each animal.The first injection was given in the superficial epithelial layer and 4 h later, a second injection was given in the deeper muscular layers of the tongue.Animals were allowed to recover in their cages for 2 days before harvesting trigeminal ganglia (TG) tissues.
Preparation of single-cell suspension and flow sorting.Tongue-innervating sensory neurons were isolated from 4 groups: Male normal (MN), Female normal (FN), Male tumor (MT) and Female tumor (FT).Neurons were isolated by preparation of single-cell suspensions of TG tissues followed by flow sorting of WGA + cells.Each sample was prepared by pooling 4 TGs from 2 animals for normal groups and 4 animals for tumor-bearing groups as only ipsilateral TGs were collected from this group.A total of 3 samples per group was prepared.TG tissues were dissected and collected in cold 1X HBSS buffer, washed three times with 1X HBSS and incubated with 5 ul of 50 ng/ml dispase (type 2, Sigma) and 75 ul of 2.5 mg/ml liberase (Roche) for 60 min at 37 °C for enzymatic digestion.Following this, tissues were centrifuged at 2 min at 1000 rpm and washed with 5mls of DMEM containing 5% FBS and resuspended in 1.5 ml DMEM with 5% FBS and triturated with a Pasteur pipette to breakdown the tissues and prepare a homogenous solution.The solution was then strained with a 100um strainer to remove all debris, supernatant collected in an eppendorf tube and subjected to flow sorting.Flow sorting was performed using FACSAria III (BD Biosciences; San Jose, CA) using 130 μm nozzle.Consecutive gates were used to isolate WGA-AF488 labeled TG neurons.First, debris was excluded by forward scatter area (FSC-A) and side scatter area (SSC-A) gating.Second, duplets and clumps were excluded by side scatter width (SSC-W) and side scatter area gate (SSC-A) gate.Third, WGA-AF488 + bright cells were gated compared to unstained TG control and sorted directly to RLT buffer (Qiagen) containing 1% 2-mercaptoethanol (Sigma) to be able to use for RNA extraction.
RNA isolation.RNA was extracted using Qiagen RNeasy micro kit (Qiagen) with on column DNase I digestion according to manufacturer's instructions.RNA quantity was evaluated using Agilent 2100 Bioanalyzer RNA 6000 Nano chip (Agilent Technologies, Santa Clara, CA).
Bulk RNA sequencing.Bulk-RNA sequencing of samples was conducted at the UTHSCSA Genome Sequencing Facility.RNA integrity was determined using Fragment Analyzer (Agilent, Santa Clara, CA) prior to library preparation.All samples were ensured to have RIN values > 6.5 to proceed with library preparations.RIN values for each sample is listed in Supplementary Table 1.RNA-seq libraries were prepared according to SMART-seq2 protocol 122 , with the following modifications: PCR preamplification to 10-12 cycles, two rounds beads cleanup with 1:1 ratio after cDNA synthesis, and 0.6-0.8dual beads cleanup for Nextera XT DNA-seq library purification.RNA-seq libraries were sequenced using Illumina HiSeq 3000 system (Illumina, San Diego, CA) with 50 bp single-read sequencing module.Upon sequencing completion, short read sequences from RNAseq were first aligned to UCSC mm9 genome build using TopHat2 aligner and then quantified for gene expression by HTSeq to obtain raw read counts per gene, and then converted to RPKM (Read Per Kilobase of gene length per Million reads of the library).
Analyses of RNA sequencing data.Sequencing data was analyzed as previously described 123 .Briefly, differential expression analysis was performed using DESeq 124 algorithm (R/Bioconductor) to estimate the differential expression in read counts and their statistical significance.Significantly differentially expressed genes were selected based on following criterion: 1) RPKM > 5 in one of the comparison group, 2) fold-change > 1.5, and 3) differential expression p-value < 0.05 along with False-discovery rate (FDR) for each gene.When comparing female tumor samples with male tumor samples, all genes coded in sex-chromosomes (chrX and chrY) were excluded.Functional assessment of DEGs were performed by using over-representation statistic (Fisher's Exact test) for Gene Ontology using Panther platform.
Immunohistochemistry. Protocol for immunostaining is described previously by us 13,36  www.nature.com/scientificreports/Z-stack images were acquired of the V3 region of TGs from at least 2 animals per group and a total of at least 6-10 images were taken for each antibody combination per group.All images were obtained with a 20 × objective at fixed acquisition parameters across all groups and were unaltered from that initially taken.Laser gain settings were determined such that no-primary control did not show any positive staining.Quantitation was achieved using Adobe Photoshop 2023 for number of neurons above threshold, in the V3 region of the TG tissue in each image.A total of 200-500 neurons per analyses were counted.
Primary culture of TG neurons.Mouse TG neurons were cultured for calcium imaging experiments.
Naïve male and female mice were injected with WGA as described above, bilaterally in the tongue and 2 days later, TGs were dissected for neuronal cultures as described previously 125 .Tissues were washed in HBSS and then dissociated with 1 mg/ml of collagenase and dispase (Roche, Indianapolis, IN, USA) for 45 min at 37 degrees.Following washing out the enzymes, cultures were plated on poly-D-lysine/laminin coated glass coverslips (BD Biosciences, San Jose, CA, USA) and grown overnight in either DMEM media supplemented with glutamine, penicillin/streptomycin and 2% Fetal bovine serum (for vehicle group) or the same media with 100 ng/ml recombinant colony-stimulating factor 1 (Csf1, Cell Signaling, Boston, MA, USA).Csf1 stock was diluted in saline.
Calcium imaging.Calcium accumulations of cultured sensory neurons was determined by fluorescence imaging as described previously 125 .TG cultured were incubated with calcium-sensitive dye, Fura-2 AM (2 μm; Molecular Probes, Carlsbad, CA, USA) in Hanks modified buffer.Imaging was performed with a Nikon TE 2000U microscope fitted with a × 40/1.35NA Fluor objective.Following selection of small-to-medium neurons and baseline measurements, cells were applied with 15 nM capsaicin (CAP) to record TRPV1 responses.The net changes in calcium influx were calculated by subtracting the average baseline intracellular calcium [Ca2 +]i level from the peak [Ca2 +]i value achieved after exposure to CAP.
Feeding behavior.Feeding behavior was performed as previously described by us 27,36 .Briefly, animals were food deprived for 18 h in their cages with bedding, following which, animals were placed in individual cages without bedding to acclimate for 30 min to an hour.Premeasured standard chow was then provided to the animals to feed freely for 1 h.At the end of the hour, food intake was measured by calculating the difference (in g) from the initial amount and the remaining amount of food.Baseline measurements and food intake at day 15 post HSC3 cell inoculation was determined for male and female mice.
vonFrey threshold in the face.Mechanical thresholds in the vibrissal pad were determined using von Frey filaments as previously described by us 27,36 .Briefly, animals were first acclimatized to von Frey filaments a day prior to experimentation by exposing the animals to the lowest force of the filament.Following this, mechanical responses at baseline and at day 15 post-cell inoculation were recorded by applying sequential series of calibrated von Frey filaments (bending forces ranging from 0.008 to 6 g of force) to the vibrissal pad of each mouse.Each filament was applied 5 times and the percent response was recorded for each filament.Head withdrawal, face swipe, or grimace were considered a positive response.Stimulus-response curve was plotted and analyzed via non-linear regression analysis.EF 50 was calculated as the force at which animals show a 50% response rate.
Tumor volume measurement.Tumor volumes were measured at day 15 post HSC3 cell inoculation in male and female mice as described previously 36 .Animals were anesthetized with 2% isoflurane inhalation and tongue tumors were measured by determining length, width and height of the tumors using calipers.Tumor volumes were calculated by the formula π/6 (length × width × height) and presented as mm3.
Statistical analyses.Sample sizes were calculated using G-power application to obtain 80% power at a two-sided tail with α error probability of 0.05.All animals were allocated to groups using simple randomization.All statistical analyses and graphical representations were performed in GraphPad Prism 9.0.Data are presented as mean ± standard error of the mean (SEM) or as heatmaps and volcano plots for RNA sequencing data.Statistical significance was determined using either Paired Student T-test, Unpaired Student's T -test, One-way ANOVA or Two-way ANOVA with Sidak's post-hoc test and p < 0.05.

Figure 1 .
Figure 1. Isolation and Estimation of Tongue-Innervating Sensory Neurons.(A, B and C) Male and Female mice were injected with 3 × 10^5 HSC3 cells in the tongue and at day 13 post-cell inoculation, tongue tissues were bilaterally injected with 1% WGA-488.Normal group received no HSC3 cells.Two days later, TG tissues were dissected to make single-cell suspension and subjected to flow sorting.Animals were grouped into male normal (MN), male tumor (MT), female normal (FN) and female tumor (FT).N = 3 samples per group.(A) Gating Strategy for flow sorting is shown.(B) Number of cells sorted for each group.Data points represent numbers of cells in each sample.(C) Percentage of WGA + neurons of total events in each group.Data points represent percentage of cells in each sample.Data are represented as mean ± SEM and analyzed by oneway ANOVA with Sidak's post-hoc test.p< 0.05.(D and E) Naïve male and female mice were injected with 1%WGA-488 and 2 days later, TG tissues were harvested for immunohistochemistry. Images were taken at 20 × magnification using the Nikon C1 confocal microscope.N = 2 mice per group.(D) Percentage of TRPV1 + / WGA + and TRPV1-/WGA + TG neurons in naïve males and females.Each data point represents average percentage of neurons per mouse.A total of 10 images were taken per mouse.Data in bar graphs are represented as mean ± SEM and analyzed by one-way ANOVA with Sidak's post-hoc test.(E) Representative images of WGA and TRPV1 staining of TG tissues in male and females are shown.Arrows indicate colocalization of TRPV1 and WGA.

Figure 2 .Figure 3 .
Figure 2. Differentially Expressed Genes in Male versus Female Lingual Neurons.(A-D) Flow sorted neurons were subjected to RNA sequencing.DEGs were identified in normal male and female mice based on RPKM > 5, FC > 1.5 and p < 0.05.N = 3 per group.(A) Heatmap of genes that were expressed exclusively in tongueinnervating neurons of normal male and female mice.(B) Number of genes in males and females is depicted as bar graphs.(C) Heatmap of genes differentially expressed in female normal (FN) versus male normal (MN).(D) Bar graphs shows number of genes upregulated and downregulated in FN versus MN. (E-J) Validation of RNA sequencing data in normal male and female by Immunohistochemistry. N = 2 mice per group.Images taken with C1 Nikon Confocal Microscope at 20 × magnification.(E) Percentage of WGA + neurons expressing Csf1R in males and females.Each data point represents average percentage of neurons per mouse.6-7 images were taken per mouse.Data are represented as mean ± SEM and analyzed by unpaired Student's T Test at p < 0.05.(F) Percentage of Csf1R in lingual TRPV1 + and TRPV1-neurons.Each data point represents average percentage of neurons per mouse.Data are represented as mean ± SEM and analyzed by paired Student's T Test at p < 0.05 (G) Representative images of immunostaining of Csf1R and TRPV1 in WGA + neurons in males and females.White arrows indicate colocalization of Csf1R and WGA whereas orange arrows indicate colocalization of Csf1R, WGA and TRPV1.(H) Percentage of WGA + neurons expressing C1qa in males and females.Each data point represents average percentage of neurons per mouse.6-7 images were taken per mouse.Data are represented as mean ± SEM and analyzed by unpaired Student's T Test at p < 0.05 (F) Percentage of C1qa in lingual TRPV1 + and TRPV1-neurons.Each data point represents average percentage of neurons per mouse.Data are represented as mean ± SEM and analyzed by paired Student's T Test at p < 0.05 (G) Representative images of immunostaining of C1qa and TRPV1 in WGA + neurons in males and females.White arrows indicate colocalization of C1qa and WGA whereas orange arrows indicate colocalization of C1qa, WGA and TRPV1.◂

Figure 4 .
Figure 4. Effect of tongue tumor on transcriptomic profile of lingual neurons in male and female mice.DEGs were identified by conducting two comparisons.(A).Volcanic Plots for all genes for MT versus MN comparison.DEGs identified are colored showing downregulated genes on the left and upregulated genes on the right.(B).Number of upregulated and downregulated genes are plotted as bar graph.(C).Top three upregulated and downregulated genes in MT versus MN are plotted as heatmap as well as tabulated for values of RPKM, Fold change (FC) and p-value.Data for heatmap plotted as fold change.Similarly, (D).Volcanic plots for FT versus FN. (E).Bar graph for upregulated and downregulated genes for FT versus FN. (F).Top 3 three upregulated and downregulated DEGs in FT versus FN as heatmaps and tabulated for RPKM, FC and p-values.Data in heatmap plotted as fold change.(G).Venn Diagram indicating number of overlapping and non-overlapping between MT versus MN and FT versus FN.

Figure 5 .
Figure 5. Biological Processes and Function of tongue-tumor controlled genes in males and females.Panther Pathway Analysis platform was used to elucidate biological processes from (A) DEGs upregulated in MT versus MN, (B) DEGs upregulated in FT versus FN and (C).DEGs downregulated in FT versus FN.For each analysis, the number of DEGs associated with each biological process is plotted as bar graphs.(D) Additional analyses was conducted to identify molecular function (MF) of all DEGs in males and females.Pie chart of the number of DEGs identified as transcription factors, ligands/receptors/growth factors (GFs), Channels/Receptors, Enzymes and Cytokines/Chemokines is shown.(E) Heatmap of select DEGs for each of the molecular function indicate the differences in expression between MT versus MN and FT versus FN.Data plotted as fold change.(F) Number of total DEGs obtained for each MF with MT versus MN and FT versus FN. https://doi.org/10.1038/s41598-023-40380-6 https://doi.org/10.1038/s41598-023-40380-6www.nature.com/scientificreports/In vivo orthotopic xenograft tongue tumor model.Using human oral squamous cell carcinoma cell

Table 1 .
List of Primary Antibodies used for Immunohistochemistry.