ErbB Signaling Pathway Genes Are Differentially Expressed in Monozygotic Twins Discordant for Sports-Related Concussion

Abstract Transcriptional changes involved in neuronal recovery after sports-related concussion (SRC) may be obscured by inter-individual variation in mRNA expression and nonspecific changes related to physical exertion. Using a co-twin study, the objective of this study was to identify important differences in mRNA expression among a single pair of monozygotic (MZ) twins discordant for concussion. A pair of MZ twins were enrolled as part of a larger study of concussion biomarkers among collegiate athletes. During the study, Twin A sustained SRC, allowing comparison of mRNA expression to the nonconcussed Twin B. Twin A clinically recovered by Day 7. mRNA expression was measured pre-injury and at 6 h and 7 days postinjury using Affymetrix HG-U133 Plus 2.0 microarray. Changes in mRNA expression from pre-injury to each postinjury time point were compared between the twins; differences >1.5-fold were considered important. Kyoto Encyclopedia of Genes and Genomes identified biologic networks associated with important transcripts. Among 38,000 analyzed genes, important changes were identified in 153 genes. The ErbB (epidermal growth factor receptor) signaling pathway was identified as the top transcriptional network from pre-injury to 7 days postinjury. Genes in this pathway with important transcriptional changes included epidermal growth factor (2.41), epiregulin (1.73), neuregulin 1 (1.54) and mechanistic target of rapamycin (1.51). In conclusion, the ErbB signaling pathway was identified as a potential regulator of clinical recovery in a MZ twin pair discordant for SRC. A co-twin study design may be a useful method for identifying important gene pathways associated with concussion recovery.

Researchers have long relied on gene expression changes identified in animal models to develop targeted therapy for traumatic brain injury (TBI), but so far these therapies have not proved effective in human trials (Xiong et al., 2013). While these preclinical models have provided valuable mechanistic insights, they might not reflect the unique pathophysiology of concussion in humans (Boomsma et al., 2002). Identifying important and potentially modifiable gene expression changes in humans diagnosed with TBI may be essential for developing effective therapies.
The inability to study expressional changes in injured brain tissue has complicated efforts to identify important gene changes in humans with TBI and in particular with mild TBI (concussion). This obstacle has been largely overcome by the recognition that gene changes in peripheral blood monocyte cells (PBMCs) are an acceptable albeit imperfect reflection of brain transcriptional changes (Sullivan et al., 2006). However, an equally important obstacle, one that has not yet been overcome, is the large degree of genetic 'noise' that can obscure identification of a TBI-related transcriptional 'signal', when genetically unrelated individuals with and without injury are compared. This noise comes primarily from natural inter-individual variation in gene transcription as well as variation in environmental exposures affecting transcription.
Comparing mRNA expression before and after injury on a subject-specific level is one way to minimize interference from this genetic noise. Athletes are a population particularly well-suited to such studies, as they are readily identifiable at baseline and have a known, measurable risk for concussion. Another way to minimize genetic noise is to compare gene changes in genetically related individuals, such as twins. Monozygotic (MZ) twins have less variation in gene expression (0−1.76%) compared to genetically unrelated individuals (up to 14.13%; Sharma et al., 2005). Twin studies examining outcome differences within a twin pair discordant for exposure (co-twin control study) have been used to investigate a variety of potentially harmful brain exposures (e.g. combat stress, trauma, childhood sexual abuse; Brown et al., 2014;Kendler et al., 2000;Koenen et al., 2003;Mitchell et al., 2018), but not head injury. The combination of a co-twin control study design with a before-and-after study design would theoretically maximize the identification of gene changes unique to concussion. The serendipitous occurrence of a concussion occurring in one member of a pair of MZ twins participating in larger study of biomarkers in concussed collegiate athletes afforded us the unique opportunity to identify novel transcriptional changes with reduced genetic noise. Because the twins reported here were roommates and teammates, genetic noise may have been further reduced by their exposure to a very similar daily environment.
The objective of the current study was to identify important and potentially novel differences in messenger RNA (mRNA) expression among this single pair of MZ twins discordant for sportrelated concussion (SRC).

Materials and Methods
A single pair of right-handed MZ twins was recruited as part of a National Institutes of Health funded study of National Collegiate Athletic Association (NCAA) athletes at the University of Rochester and Rochester Institute of Technology designed to validate bloodbased markers of concussion. During the study, one twin (Twin A) sustained a concussion while playing soccer, affording us the rare opportunity to compare mRNA expression to the uninjured twin (Twin B).
The study methods have been described in detail elsewhere (Merchant-Borna et al., 2016). In brief, athletes were eligible for inclusion if they were ≥18 years of age and played an NCAA varsity contact sport. A baseline, preseason (T 0 ) assessment included a venous blood sample, cognitive performance (ImPACT test) and postural stability (Balance Error Scoring System) on all participating athletes. Athletes were followed prospectively for development of SRC as defined by the Third International Conference on Concussion in Sport (McCrory et al., 2009). After Twin A sustained a SRC, a second blood sample was obtained from both twins within 6 h of the injury (T 1 or 'acutely' postinjury). White matter integrity was measured with diffusion tensor imaging (DTI) at 2 days and 1-year postinjury. Cognitive performance, balance and symptoms were remeasured 3 days postinjury. A third blood sample, repeat cognitive performance and balance measures were obtained at 7 days postinjury (T 2 or 'subacutely' postinjury).
The institutional review boards at University of Rochester and Rochester Institute of Technology approved this protocol; written informed consent was obtained from study participants prior to subject participation. Supplementary Material is available on the Cambridge Core website

PBMC and mRNA Isolation
The methods to isolate mRNA from PBMCs have been previously described (Merchant-Borna et al., 2016). In summary, PBMCs were isolated within 1 h of blood collection (Fuss et al., 2009), and the resulting pellets were suspended in RPMI-10 medium and stored at −190°C until analysis. RNA was isolated from PBMC pellets using TRIzol® Plus RNA Purification Kits (Life Technologies, Grand Island, NY) followed by DNase I-Amplification Grade Kits (Life Technologies). A NanoDrop DN-1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE) was used to verify the purity and concentration of the RNA samples. The integrity of the RNA samples was measured on an Agilent bioanalyzer 2100 using the RNA 6000 Nano Kit (Agilent, Sana Clara, CA). Only samples with an RNA integrity number (RIN) greater than 7.0 were used in analysis; samples below this threshold were excluded. Each RNA underwent reverse transcription using a GeneChip® 3' IVT Expression Kit, converting the RNA strands to biotinylated complementary RNA, which was hybridized to Affymetrix HG-U133 Plus 2.0 microarray (Affymetrix, Sana Clara, CA), stained with streptavidin-phycoerythrin and washed. Raw data were generated by laser scanning imaging.

DTI Methods
Both twins were scanned on a Siemens 3T Tim Trio scanner with a 32-channel head coil using a single-shot spin-echo echoplanar imaging (60 diffusion directions with b = 1200 s/mm2, 10 images with b = 0 s/mm2, TR = 8900 ms, TE = 86 ms, FOV = 256 Å∼ 256 mm2, matrix = 128 Å∼ 128, voxel size = 2 mm3 by 2 mm3 by 2 mm3, 70 axial slices). A double-echo gradient echo field map sequence (echo time difference = 2.46 ms, EPI dwell time = 0.75 ms) was acquired with the same resolution as the DTI sequence and was used to correct for distortion caused by B0 inhomogeneity. Diffusion-weighted images in conjunction with an atlas-based Markov random field representation were used to define 31 WM structures (Bazin et al., 2011). FA values were measured in these 31 regions. For each ROI, the FA difference between Twin A and Twin B was calculated and compared to the 95% CI for the mean FA difference among 65 uninjured female MZ twins (mean age 22.9 þ 2.1 yrs) scanned between 2007 and 2012 as part of the Queensland Twin Imaging Study (QTIM;De Zubicaray et al., 2008;Kochunov et al., 2014). ROIs in which FA differences between the twins were either above or below this 95% CI were considered indicative of axonal injury. The QTIM used the Johns Hopkins University white matter parcellation atlas (Mori et al., 2008) to define 41 ROIs, 15 of which were identical to the ROIs defined in the study twins using the atlas-based Markov method.

Analysis
Partek Genomics Suite Software, version 6.6 (Partek Inc., St. Louis, MO), was used to analyze the microarray data. The robust multiarray average method was used to correct for background signaling once interrogating probes were imported. Additional corrections were applied to the GC content probes. Quantile normalization was used to standardize the probesets, and the data distribution was normalized by transforming gene expression on a log-2 scale.
In order to isolate the effect of concussion on the gene expression in Twin A, we first calculated the change in mRNA expression for each twin across each time point; pre-injury to 6 h postinjury (T 1 -T 0 ), pre-injury to 7 days postinjury (T 2 -T 0 ) and 6 h postinjury to 7 days postinjury (T 2 -T 1 ) (Supplementary Table S1). The difference in the longitudinal change of mRNA expression was then determined between the twins. Taking the difference in the longitudinal change maximizes detection of gene changes unique to brain injury. Important genes were identified based on a foldchange threshold value of greater than 1.5 or less than −1.5 for all three time point comparisons.

Pathway Enrichment Analysis
Gene set enrichment analysis by hypergeometric test was performed in R (R Core Team, 2013), version 1.2.1335, using Biocarta (Nishimura, 2001), KEGG (Ogata et al., 1999), and Reactome (Vastrik et al., 2007) gene sets from msigdb (Liberzon et al., 2015) and gmt utilities from qusage (Liberzon et al., 2015) with a conservative estimate of 20,000 for size of universe. P values were adjusted using the Benjamini-Hochberg method (Benjamini & Hochberg, 1995) and adjusted p values < .2 output. Cytoscape (Shannon et al., 2003) was used to visualize ErbB KEGG pathway. The gene set enrichment analysis was confined to the baseline-tosubacute time interval (T 2 -T 0 ) in order to understand the genetic changes that may be associated with clinical recovery.

Clinical and White Matter Changes
Twin A was diagnosed with an SRC while participating in collegiate soccer. She was struck in the face by a ball kicked at close range and was witnessed to have a loss of consciousness of several seconds. She was amnestic for the event and reported headache, dizziness, lightheadedness. She was immediately removed from play. Both twins underwent DTI scanning 2 days postinjury. Of 31 ROIs, we were able to compare 15 ROIs to 65 female MZ twin pairs from the QTIM study. Of these 15 ROIs, FA differences between the study twins were outside the regional 95% CI in 5 ROIs (33%). Clinical testing 3 days postinjury revealed a marked worsening in verbal memory, visual motor speed and reaction time but not postural stabilityand increase in concussive symptoms, compared to her personal baseline and compared to testing done simultaneously on Twin B (Figure 1).
By 7 days postinjury, Twin A's symptoms, verbal memory, visual motor speed and reaction time had returned to baseline levels suggesting clinical recovery. The improvement in Twin A's visual memory score between 3 and 7 days postinjury is also consistent with concussion recovery. However, why Twin A's baseline visual memory score was lower than that of Twin B and why Twin B's score declined between 3 and 7 days in the absence of concussive injury are unclear. Possibilities include sleep deprivation, pain and emotional stress. By Day 7, Twin A was cleared to begin the gradual return to play progression, which she successfully completed on day 13 postinjury. Twin B continued sport participation throughout the duration of this study. By 1-year postinjury, FA differences between the study twins returned to within the 95% CI in 4 of the 5 ROIs that had been outside the 95% CI at 2 days (not shown).

Messenger RNA Expression Changes
54,675 probesets representing 38,000 genes were compared for important differences in longitudinal mRNA expression between the twins. Changes in 158 probes representing 153 unique genes had either a greater than 1.5 or less than −1.5 difference in the longitudinal change between the twins and were thus considered important (Figure 2). The majority (144)    Heat map of the difference in longitudinal change in mRNA expression between single MZ twin pair discordant for concussion across the three time intervals (Xaxis): baseline to acute (T 0 to T 1 ), baseline to subacute (T 0 to T 2 ) and acute to subacute (T 1 to T 2 ). The Y-axis represents individual genes, color-coded by upregulation (red) or downregulation (blue) were protein-coding; however, 2 antisense RNA transcripts and 7 noncoding RNA transcripts also were changed. Considering all time points, more probesets were upregulated (111 probesets) than downregulated (68 probesets).

Discussion
In the current study, we identified the ErbB signaling pathway as a potential regulator of recovery in single pair of MZ twins discordant for concussion. The chance occurrence of concussion in one member of a pair of MZ twins provided the rare and unique opportunity to identify novel concussion-related transcriptional changes in the setting of reduced genetic noise. An important study strength was the acquisition of PBMCs from each twin before injury, due to their participation in a larger NIH-funded study. This allowed us to further isolate the effects of concussion by comparing pre-and postinjury transcriptional changes in each twin.
Another key strength was the use of structural neuroimaging to verify the diagnosis of concussion. Although this diagnosis was based on Twin A's symptoms at the time of injury, the finding of abnormal FA difference between the twins in five white matter brain regions acutely postinjury increased confidence in the validity of a concussion diagnosis and suggests her injury resulted in reduced white matter integrity. The resolution of the FA differences in four of these five regions 1 year after clinical recovery further supports that these ROIs were injured acutely.
Although we observed important changes in single transcripts at both the acute and subacute time points, the pathway enrichment analysis of changes at the subacute time point, when Twin B was back to baseline clinically, is likely a more robust indicator of biologic pathways involved in concussion recovery. Pathway enrichment analysis has several advantages over analyzing single genes. Statistical power is improved by reducing the dimensionality of large gene lists into a smaller number of systems or pathways (Reimand et al., 2019). Moreover, pathway analysis simplifies data interpretation by constructing familiar concepts (e.g. biological processes) from the gene set, facilitating identification of putative causal mechanisms and drug targets (Reimand et al., 2019). Previous concussion research has taken advantage of pathway analysis to better interpret and understand genetic data (Johnson et al., 2018;LaRocca et al., 2019;Merchant-Borna et al., 2016).
Our pathway analysis of mRNA expression between pre-injury baseline and 7 days postinjury suggested the ErbB signaling pathway to be a potential regulator of the biologic processes involved in recovery following SRC. Specific genes in this pathway with important longitudinal changes in the twins included EREG, NRG1, MTOR and EGFR. Although insufficient ErbB signaling is known to be associated with neurodegenerative diseases, such as Alzheimer's disease and multiple sclerosis, its role in TBI recovery has been less frequently reported (Bublil & Yarden, 2007). In an animal model of TBI, Erlich et al. (2000) demonstrated significant upregulation of ErbB4 receptor (which binds NRG1) in neuronal and glial cells at the site of injury. LaRocca et al. (2019) reported increased expression of salivary microRNAs targeting transcripts of the ErbB pathway among 50 amateur mixed martial artists exposed to subconcussive head impacts during a fight (Atif & Hicks, 2019).
ErbB signaling consists of four homodimerized, or heterodimerized, receptors that regulate a variety of pathways including PI3K/AKT/mTOR, p38/MAPK, calcium signaling and others (Yarden & Shilo, 2007;Kataria et al., 2019). Previous TBI studies have targeted the PI3K/AKT/mTOR and p38/MAPK pathways post-TBI in order to improve functional outcome (Bachstetter et al., 2013;Erlich et al., 2007). Postmortem brain tissues of CTE patients have demonstrated a number of genes related to the MAPK and calcium signaling pathway to be significantly downregulated (Seo et al., 2017). Yet, studies have not extensively investigated the role of the ErbB signaling in regulating these downstream pathways postconcussion.
We identified NRG1, a growth factor that binds ErbB receptors, to be upregulated from pre-injury baseline to 7 days postinjury. The exact role this gene plays in concussion recovery cannot be determined from our results, but evidence from other disease processes (spinal cord injury, brain ischemia, Alzheimer's disease, Parkinson's disease) suggests that NRG1 controls several critical functions that could potentially promote concussion recovery including myelination, blood-brain barrier (BBB) function, synaptic function and immune response.  NRG1 has been shown to protect microvasculature and reduce permeability of the BBB via signaling in brain endothelial cells (Lok et al., 2009;Lok et al., 2012), which could mitigate disruption of the BBB following TBI (Blyth et al., 2009). In the postconcussion period, white matter tracts are susceptible to injury and demyelination, which is characterized by loss of myelin sheath and oligodendrocyte cell death (Mierzwa et al., 2015;Shi et al., 2015). NRG1 plays an important role in the development, maturation and survival of oligodendrocytes and seems to affect the quality of myelination in the CNS (Shahriary et al., 2019). Additionally, NRG1 serves a neuroprotective role through anti-inflammatory actions, including the regulation of microglia, astrocyte and immune cell response after in vitro CNS injury (Alizadeh et al., 2017;Shahriary et al., 2019). These studies provide evidence for a possible role of NRG1 in the brain's recovery process following concussion through the ErbB signaling pathway.

Limitations
An obvious limitation of the study is sample size, with only one twin pair evaluated, which increases the uncertainty that the important genetic changes we identified are unique to concussion. Additionally, as the MZ twins were athletes, the generalizability of results must be carefully considered, as physical exertion itself has been found to alter gene expression (Connolly et al., 2004). While  (Ogata et al., 1999) Genes with important differential change between the twins discordant for concussion from baseline to subacute (T 0 to T 2 ) are shown in either blue (downregulation) or red (upregulation). Genes of the ErbB pathway that did not have important upregulation or downregulation are indicated in green. Receptors of the ErbB pathway are either homodimerized or heterodimerized, demonstrated in the figure by the blocks containing two gene symbols (e.g. ErbB4/ErbB4). Important upregulated genes (red) included epiregulin (EREG), neuregulin 1 (NRG1) and mechanistic target of rapamycin (MTOR). Epidermal growth factor (EGFR) was downregulated (blue). The þP symbol indicates a gene in which the encode protein is responsible for a phosphorylation event in the ErbB pathway the concussed twin was prohibited from returning to play until recovered, the control twin continued to engage in physical activity. Therefore, some of the important gene changes we identified may reflect not only processes related to brain injury but also to physical exertion. Furthermore, epigenetic divergence of MZ twins could contribute to the differences in gene expression observed (Fraga et al., 2005), with additional genetic variability in female twin pairs (such as the pair in the present study) caused by random x-inactivation and/or reproductive hormone cycling.

Conclusions
By evaluating a set of MZ twins discordant for concussion and reduced white matter integrity, we identified the ErbB signaling pathway as the top differentially expressed gene network from pre-injury to 7 days postinjury. Since the concussed twin appeared clinically recovered by 7 days, we believe that ErbB signaling may be important for concussion recovery. Understanding genetic regulation of concussion recovery can help identify biologic processes that may be modifiable and suitable targets for future therapies. Because MZ twins have less inter-individual difference in gene expression than the general population, focusing on them as study subjects may facilitate detection of gene networks specific to concussion. A co-twin study design may thus be a new avenue for identifying modifiable targets with an improved chance of leading to effective concussion treatments. Data. The full gene data set used in this study has been publicly available on Gene Expression Omnibus (GEO) repository since December 31, 2021, under the following code/name: GSE176388 -Expression data from MZ Twins discordant for concussion