Schwann cell nodal membrane disruption triggers bystander axonal degeneration in a Guillain-Barré syndrome mouse model

In Guillain-Barré syndrome (GBS), both axonal and demyelinating variants can be mediated by complement-fixing anti–GM1 ganglioside autoantibodies that target peripheral nerve axonal and Schwann cell (SC) membranes, respectively. Critically, the extent of axonal degeneration in both variants dictates long-term outcome. The differing pathomechanisms underlying direct axonal injury and the secondary bystander axonal degeneration following SC injury are unresolved. To investigate this, we generated glycosyltransferase-disrupted transgenic mice that express GM1 ganglioside either exclusively in neurons [GalNAcT–/–-Tg(neuronal)] or glia [GalNAcT–/–-Tg(glial)], thereby allowing anti-GM1 antibodies to solely target GM1 in either axonal or SC membranes, respectively. Myelinated-axon integrity in distal motor nerves was studied in transgenic mice exposed to anti-GM1 antibody and complement in ex vivo and in vivo injury paradigms. Axonal targeting induced catastrophic acute axonal disruption, as expected. When mice with GM1 in SC membranes were targeted, acute disruption of perisynaptic glia and SC membranes at nodes of Ranvier (NoRs) occurred. Following glial injury, axonal disruption at NoRs also developed subacutely, progressing to secondary axonal degeneration. These models differentiate the distinctly different axonopathic pathways under axonal and glial membrane targeting conditions, and provide insights into primary and secondary axonal injury, currently a major unsolved area in GBS research.

Introduction
Guillain-Barré syndrome (GBS) is a postinfectious autoimmune disorder affecting the peripheral nervous system (PNS) resulting in acute-onset paralysis (1). GBS comprises a spectrum of axonal and demyelinating variants that exhibit distinct pathophysiological features. The 2 major variants are acute motor axonal neuropathy (AMAN) and acute inflammatory demyelinating polyneuropathy (AIDP), a classification derived from electrophysiological and morphological features and characterized by axonal injury and segmental demyelination, respectively. AMAN cases display either reversible conduction failure followed by rapid recovery, understood to be a consequence of axonal conduction block without axonal transection at the nodes of Ranvier (NoRs) and distal motor nerves (2,3), or a severe outcome associated with inexcitable nerves and extensive proximal and distal axon degeneration (4). In contrast, AIDP cases exhibit conduction slowing and/or conduction block, normally considered to be associated with segmental demyelination, and usually fully recover following regeneration of the myelinating Schwann cell (SC). Recent electrophysiological data have highlighted the limitations of the AMAN/AIDP classification by introducing the alternative term "nodo-paranodopathy" to describe clinical situations in which either axonal or SC nodal membrane injury (or both) account for the acute conduction failure, in the absence of segmental demyelination as an early feature (3). The direct pathological evidence of this electrophysiological conceptualization is incomplete, and the clinical outcome varied. Apart from primary axonal injury, secondary bystander axonal degeneration may occur following SC membrane injury, either locally at the nodal region or with segmental demyelination. The long-term severity of GBS is dictated by the extent of the primary and/or secondary bystander axonal degeneration, supported by studies indicating that high serum levels of the axonal structural protein, neurofilament-light chain, can be predictive of poorer patient outcome (5). As such it remains critical to prognostic modeling and therapy development to differentiate and understand the pathological and degenerative pathways in action at the axo-glial interface in GBS, especially the mechanisms underlying secondary axonal degeneration in AIDP-and SC-restricted forms of GBS.
Evidence from patient serology and human postmortem nerve material indicates that autoantibody and complement system involvement with macrophage recruitment are prominent mediators of GBS (6)(7)(8). The best-understood antigens that underlie the postinfectious autoimmune response in GBS are glycans borne by bacterial lipo-oligosaccharides (LOS), In Guillain-Barré syndrome (GBS), both axonal and demyelinating variants can be mediated by complement-fixing anti-GM1 ganglioside autoantibodies that target peripheral nerve axonal and Schwann cell (SC) membranes, respectively. Critically, the extent of axonal degeneration in both variants dictates long-term outcome. The differing pathomechanisms underlying direct axonal injury and the secondary bystander axonal degeneration following SC injury are unresolved. To investigate this, we generated glycosyltransferase-disrupted transgenic mice that express GM1 ganglioside either exclusively in neurons [GalNAcT -/--Tg(neuronal)] or glia [GalNAcT -/--Tg(glial)], thereby allowing anti-GM1 antibodies to solely target GM1 in either axonal or SC membranes, respectively. Myelinated-axon integrity in distal motor nerves was studied in transgenic mice exposed to anti-GM1 antibody and complement in ex vivo and in vivo injury paradigms. Axonal targeting induced catastrophic acute axonal disruption, as expected. When mice with GM1 in SC membranes were targeted, acute disruption of perisynaptic glia and SC membranes at nodes of Ranvier (NoRs) occurred. Following glial injury, axonal disruption at NoRs also developed subacutely, progressing to secondary axonal degeneration. These models differentiate the distinctly different axonopathic pathways under axonal and glial membrane targeting conditions, and provide insights into primary and secondary axonal injury, currently a major unsolved area in GBS research.
Schwann cell nodal membrane disruption triggers bystander axonal degeneration in a Guillain-Barré syndrome mouse model  (20). (B) Using a single anti-GM1 antibody (Ab, green), differential binding was observed at the distal motor nerves from triangularis sterni nerve-muscle preparations among genotypes. Open arrowheads indicate internodal Schwann cell (SC) abaxonal membrane anti-GM1 Ab deposition on WT and Glial nerves (absent along Neuronal nerves). Gliomedin (Gldn) immunostaining identifies the nodal gap. Boxed areas are enlarged underneath and represent differential anti-GM1 Ab binding at nodes of Ranvier (NoRs) among genotypes in relation to gliomedin (closed arrowheads). Dashed lines delineate the border of the axonal membrane determined by cytoplasmic CFP-positive axons. Scale bars: 10 μm (top panels) and 5 μm (lower panels). Asterisks indicate motor nerve terminals. (C) Caspr1 immunostaining (magenta) indicates the paranodes. Dashed lines delineate the border of the axonal membrane and arrowheads indicate anti-GM1 Ab binding beyond this membrane, suggesting binding on the glial membranes of the SC microvilli (open arrowheads) or paranodal loops (closed arrowheads) at WT and Glial NoRs. Scale bar: 2 μm. J Clin Invest. 2022;132(14):e158524 https://doi.org/10.1172/JCI158524 arrowheads, Figure 1B). Staining was enriched at the NoRs and colocalized with the SC microvilli marker gliomedin, and with paranodal Caspr1 at the node-proximal paranode. Binding on the glial SC microvilli and paranodal loop membranes was indicated by immunostaining beyond the border of axonally expressed cyan fluorescent protein (CFP) and Caspr1, demarcated by dashed lines and arrowheads in Figure 1, B and C. In Neuronal distal motor nerve, binding was restricted to the axolemma at the nerve terminal and NoRs; specifically, it colocalized with nodal gliomedin, was found between Caspr1 domains, was not detected beyond the demarcated axonal membrane, and was not observed anywhere along the SC abaxonal membrane. In Glial mice, internodal SC abaxonal membranes were positive, with the greatest enrichment at the paranodal loops most proximal to the nodal gap (node-proximal paranodal loops), flanking gliomedin. In some cases, staining could also be observed between Caspr1 domains and beyond the axonal membrane, which suggests SC microvilli positivity. Deposition on the axolemma was not clearly detected at the terminal or NoRs (also see Supplemental Figure 1A; supplemental material available online with this article; https://doi. org/10.1172/JCI158524DS1). A similar pattern for all genotypes was observed with cholera toxin B subunit labeling (a GM1 ligand, data not shown). As nerve terminal injury has previously been reported using Neuronal mice (21), the current research focuses on characterizing nerve terminal and NoR injury induced in Glial mice and the subsequent effects on axonal integrity.
Distal motor nerve axonal integrity remains intact in an acute ex vivo model following selective glial targeting. The TS ex vivo complement-mediated nerve-muscle injury model was used to assess axonal integrity when selectively targeting glial membranes in Glial mice compared with WT and Neuronal mice. Our previous research has shown that complement deposition shares a similar binding pattern to AGAb, indicating site-specific antibody-directed targeting of complement deposits (15). When a source of complement (provided in the form of normal human serum, NHS) was added to preparations pretreated with anti-GM1 antibody (injury), MAC deposition was observed in a neural membrane-selective pattern (Supplemental Figure 1B), as expected from the immunostaining patterns shown in Figure 1, B and C. Thus, in WT mice MAC was deposited at NoRs (arrowheads) and weakly over the MNT (asterisks) (Figure 2A). MAC staining was absent from anti-GM1 antibody-only control tissue. In Neuronal mice, MAC pore deposition was restricted to the axonal membrane and present at both the NoRs and MNTs in injured tissue. In contrast, in Glial tissue, MAC deposition was observed at the node-proximal paranodal loops, the perisynaptic SCs (pSCs) overlying the MNTs, and was more weakly seen along the internodal abaxonal SC membrane. Thus, the pattern of MAC deposition recapitulates the distribution of anti-GM1 antibody deposition, indicating that complement is being activated specifically at anti-GM1 antibody binding sites.
Controls from each genotype exposed only to anti-GM1 antibody (lacking NHS as a complement source) showed normal axonal integrity at NoRs and MNTs, judged by intact neurofilament immunostaining occupying these sites. In this acute injury (4 hours) paradigm, MAC pore deposition in Glial mice did not cause a significant loss of axonal integrity at either NoRs or MNTs compared to control ( Figure 2A). As expected, and in diametric con-at NoRs and motor nerve terminals (MNTs) (12), sites essential for saltatory conduction and synaptic transmission, respectively. The development of transgenic mice null for β4-GalNAc transferase that synthesizes complex gangliosides has demonstrated their essential role in the stability and function of the NoR (13). Significantly, the NoR is a major site of pathophysiology in GBS, first established from autopsy material that showed a lengthening of the nodal gap and axonal pathology with nodal autoantibody and complement deposition (6,7). Subsequently, an anti-GM1 antibody-driven rabbit model of AMAN replicated this nodal lengthening (14). Combined with rodent models, a causative role for anti-ganglioside antibody (AGAb) activation of the complement cascade leading to dysfunction of the nerve was thus established (14)(15)(16). These models demonstrated loss of axonal integrity, eversion of paranodal loops, and disruption of key axo-glial adhesion proteins: contactin-associated protein 1 (Caspr1) localization on the axonal membrane and neurofascin 155 (NF155) on the glial membrane, that together contribute to the formation of the normal paranodal axo-glial junction. Significantly, ankyrin G (AnkG), the axonal anchoring protein responsible for tethering and maintaining voltage-gated sodium (Nav) channels in conjunction with neurofascin 186 (NF186) to the actin-spectrin cytoskeleton at the NoR (17), was also disturbed in these models. Later, the direct axonopathic mechanism was shown to be mediated through membrane attack complex (MAC) pore formation, calcium influx (18), and activation of the calcium-dependent protease calpain (16,19) for which neurofilament, AnkG, actin, and spectrin are substrates.
Since gangliosides, including GM1, are expressed on both neuronal and glial membranes, it is impossible to separate primary and secondary injury from the consequences of a cell-specific (axonal or glial) membrane injury. Understanding these pathway differences is a crucial step in understanding the diverse degenerative mechanisms and informing the development of GBS variant-specific models and therapeutics. To unravel this, we developed transgenic mice that express gangliosides, including GM1, selectively on either glial or neuronal membranes (20). Here, we delineate axonal fate in distal motor nerve NoRs and nerve terminals, consequent to selective targeting of either SC or axonal membranes in these mice.

Results
Site-directed targeting of complex ganglioside expression to distal motor nerve membranes. In order to exclusively target axonal or glial membranes via anti-GM1 antibody-mediated injury, we developed transgenic mouse strains that have complex ganglioside expression limited to neurons (and thus also their axon projections) through the human Thy1.2 promoter [GalNAc-T -/--Tg(neuronal), Neuronal] or myelinating and non-myelinating glia through the mouse proteolipid protein (Plp) promoter [GalNAc-T -/--Tg(glial), Glial] ( Figure 1A), as previously reported (12,20,21). Firstly, we used triangularis sterni (TS) nerve-muscle preparations to study the differing binding patterns of a single monoclonal anti-GM1 antibody (DG2) (22) at distal nerves and NoRs from Neuronal and Glial mice compared to wild-type (WT) mice ( Figure 1, B and C). In WT nerve, anti-GM1 antibody binding was observed on the axonal membrane at the nerve terminal and weakly along the SC abaxonal membrane of the internode (open Given the alterations to NF186, we next studied its SC microvilli ligand gliomedin. Gliomedin immunostaining was present at the majority of NoRs across all genotypes and treatments, but with an abnormal distribution compared with control in injured WT and Glial mice ( Figure 3B) demonstrating a disruption to the SC microvilli glial membrane when directly targeted. The less dramatic changes are consistent with the lower levels of anti-GM1 antibody deposits on SC microvilli relative to paranodal loops.
We next considered that an indirect effect on Caspr1, the axonal partner of NF155, might be evident when targeting the glial membranes (Supplemental Figure 2B). Results show that Caspr1 immunostaining was disrupted or lost in injured Glial and WT tissue in which AnkB and NF155 were also absent. In contrast, Caspr1 immunostaining was relatively preserved at injured Neuronal NoRs, but with similar abnormalities as described for NF155. We attribute this distinction between Glial and Neuronal Caspr1 integrity to the preservation of AnkB and NF155 in Neuronal mice. Intriguingly, despite disturbances to the paranodal axo-glial junction, in preliminary studies we did not detect an invasion of voltage-gated potassium (Kv1.1) channels into the paranodal domain, as may be expected (Supplemental Figure 2C). Instead, there appeared to be a lengthening of the NoR between the domains in injured tissue from all genotypes. Myelin basic protein (MBP, a marker for compact myelin) intensity measurements across the NoRs (Supplemental Figure 3A) showed a reduction at the lateral node-proximal edges in injured Glial and WT mice relative to Neuronal mice. This indicates disruption to the myelin at this site or a widening of the NoR by a reduction in intensity at the margins between the domains. Differences in MBP immunostaining patterns can be observed at NoRs throughout Figure 3. There were no measurable differences in internodal MBP intensity (data not shown), but staining was qualitatively abnormal in Glial mice, being punctate and vesiculated in the co-presence of MAC pore deposition, suggesting some early myelin damage (Supplemental Figure 3B). Additionally, EthD-2-positive SC nuclei were observed after acute injury in Glial preparations, indicating initial SC injury (Supplemental Figure 3C).
To firmly establish that the observed disruptions to paranodal proteins were not partly a confounding consequence of structural dysfunction secondary to abnormal ganglioside profiles in the axonal compartment of transgenic Glial mice, a complement-fixing antibody against the myelin glycolipid sulfatide (27) was used to target the paranodal loop glial membrane in WT mice in order to look for similar perturbations. As seen with the anti-GM1 antibody in Glial mice, AnkB and Caspr1 immunostaining was lost/ abnormal compared with uninjured control (Supplemental Figure  4) at anti-sulfatide antibody-exposed paranodes.
Disruption of axo-glial adhesion proteins by glial membrane targeting induces functional consequences similar to those of axonal targeting. To assess the functional impact of disrupted anchoring proteins and CAMs ex vivo, we first studied Nav channel clustering by immunostaining ( Figure 3C). All genotypes with NHS-induced injury had significantly fewer NoRs with normal Nav channel clustering compared with controls (P < 0.001, two-way ANOVA). Notably, there were significantly more NoRs with absent Nav channel clusters in injured Neuronal tissue compared with Glial tissue (P < 0.05, two-way ANOVA), suggesting a greater severity trast to the Glial mice, complement deposition resulted in almost complete loss of axonal integrity at the NoRs and MNTs in injured WT and Neuronal tissue.
Previous studies have shown that AGAb can bind and injure pSCs overlying the nerve terminal but does not induce loss of terminal axon integrity in an acute injury model of 1-hour duration (23). We used the cell viability label, ethidium homodimer (EthD-2), to assess the health of pSCs ( Figure 2B). EthD-2-positive cell number significantly increased over MNTs from injured Glial and WT tissue compared with Neuronal tissue (P < 0.001, twoway ANOVA). Injury to pSCs alone, without direct axonal targeting, was insufficient to cause loss of axonal integrity at the MNTs at this acute time point. Taken together, these data suggest that direct targeting of the axonal membrane is necessary for acute loss of axonal integrity.
Disruption of intracellular cytoskeletal proteins at NoRs following selective glial and neuronal targeting ex vivo. The consequence of selective membrane targeting for glial and axonal cytoskeletal anchoring protein localization at NoRs was comparatively analyzed ex vivo in TS preparations exposed to anti-GM1 antibody and complement. The glial paranodal cytoskeletal anchoring protein, ankyrin B (AnkB) (24), and the axonal cytoskeletal anchoring protein, AnkG, were compared at the NoRs in Neuronal, Glial, and WT tissue ( Figure 2, C and D). When tissue was injured through the addition of NHS, AnkB immunostaining was lost in Glial and not Neuronal mice, although it appeared disrupted in the latter, with either weaker staining or a lengthened gap between domains, signifying either a lengthened nodal gap and/or a loss of NF155 from the innermost proximal-nodal border. Conversely, AnkG immunostaining was significantly absent from injured Neuronal compared with Glial NoRs ( Figure 2D, P < 0.05, two-way ANO-VA). A degree of disruption, as defined by an abnormal immunostaining pattern (outlined in Methods), to AnkB and AnkG occurred in Neuronal and Glial mice, respectively, indicating some subtle reciprocal disturbance of axoglial integrity when one or the other membrane is targeted. In WT mice in which both glial and axonal components of the NoR were targeted, both AnkB and AnkG were almost entirely absent.
Disruption of axo-glial adhesion proteins at NoRs following selective glial and neuronal targeting ex vivo. As AnkB and AnkG are key cytoskeletal proteins in the glial and axonal cytoplasmic compartments of NoRs, respectively, we next studied the disruption to the axo-glial cell adhesion molecules (CAMs: glial, NF155; axonal, NF186). A fraction of AnkB interacts with NF155 in the paranodal loops (24) and AnkG anchors NF186 to the cytoskeleton at the axonal nodal membrane (25). Caspr1 is the axonal partner of glial NF155 and together with contactin-1 these CAMs form the major component of the paranodal axo-glial junction.
A pan-neurofascin antibody (binds both NF155 and NF186) (26) was used to assess the integrity of glial NF155 and axonal NF186 in the different genotypes ex vivo ( Figure 3A and Supplemental Figure 2A). After injury, paranodal NF155 immunostaining was absent from NoRs in Glial mice, but present in Neuronal mice NoRs, albeit disrupted. In contrast, nodal NF186 immunostaining was absent from NoRs in injured Neuronal and moderately preserved in injured Glial mice. In injured WT mice, immunostaining for both NF155 and NF186 was severely disrupted or lost entirely. J Clin Invest. 2022;132(14):e158524 https://doi.org/10.1172/JCI158524 of injury to clustering when the axon is directly targeted; nevertheless, Glial injury also clearly affects Nav channel clustering, likely through an indirect pathway.
To study the consequent impact on function of losing Nav channel clustering and NoR architecture, we recorded perineural current waveforms from the distal nerves ex vivo ( Figure 3D). In control observations, we recorded Na + and K + waveforms in Glial and Neuronal tissue under exposure to either anti-GM1 antibody alone or NHS alone. Perineural currents showed normal Na + and K + waveforms, as expected, and there were no changes in the Na + /K + ratio among genotypes or treatments, showing no effect of antibody or NHS treatment alone (Supplemental Figure 3D). However, waveforms were completely absent from both injured Glial and Neuronal preparations, indicating severe functional impairment in both genotypes.
The above ex vivo results demonstrate that injury in WT leads to disruption of multiple proteins on both neuronal and glial membranes at the NoRs, as would be expected. However, selectively targeting the glial compartment principally causes SC nodal membrane disruption with mild axonal protein disruption, whereas targeting the axonal compartment principally leads to a major disruption of axonal structural and adhesion proteins, nodal lengthening, and leaves a partially intact axo-glial junction at the paranode. Significantly, either direct targeting to the axon or the SC nodal membrane has the same functional consequence: disruption to normal Nav channel clustering and acute conduction failure.
Distal motor nerve axon integrity remains intact following acute selective glial membrane targeting in vivo. Following the above characterization of the differential injury to the distal nerve in ex vivo models, we next sought evidence for selective injury in vivo. WT, Neuronal, and Glial mice were injected intraperitoneally (i.p.) with either anti-GM1 antibody and NHS (injury) or NHS only (con-trol). The diaphragm was principally affected; therefore, noninvasive whole-body plethysmography was used to monitor respiratory function in live mice and the diaphragm was subsequently removed for immunostaining or morphological analysis.
Injured Neuronal mice presented with the most severe respiratory phenotype consisting of a pinched, wasp-like abdomen indicative of a paralyzed diaphragm, as previously reported (21). In comparison, injured Glial and WT mice exhibited a milder respiratory dysfunction ( Figure 4A). Tidal volume, used as a measure of respiratory function, was reduced at 5 hours after injury in all injured groups compared with baseline and genotype-matched control, as illustrated in representative respiratory flow charts, but most severely reduced in the Neuronal mice ( Figure 4A). Additionally, the flow charts illustrate increased respiratory rate (tachypnoea) in all injured groups, although owing to a wide range in variability in individual mouse activity, this did not reach significance (data not shown). Serum from terminal bleeds was assessed for circulating anti-GM1 antibody and confirmed antibody presence in the Neuronal and Glial injury groups ( Figure 4A) but was undetectable in WT mice. This finding replicates a phenomenon we previously reported where circulating antibody removal by endocytosis was observed in WT mice, thus rendering them unsuitable for in vivo injury modeling, as insufficient antibody remains available to bind and injure the intraneural target tissue (28).
We next assessed complement deposition and axonal integrity following in vivo injury in Glial mice compared to Neuronal mice ( Figure 4B). MNTs with overlying complement deposits increased in injured Glial and Neuronal mice, reaching significance compared with their respective controls (P < 0.05 and P < 0.001, twoway ANOVA). Conversely, complement deposits were detected more readily along distal nerve SC membranes in injured Glial mice compared with all other treatment groups (P < 0.001, twoway ANOVA). Consistent with serum antibody results showing rapid AGAb sequestration, complement deposits were not detected to a significant degree compared to control in WT mice at either the MNTs or along the distal nerve. As seen in the acute ex vivo experiments, MNTs or distal nerves occupied by neurofilaments were not significantly lower in injured Glial mice compared to control, whereas injured Neuronal mice showed a significant reduction in neurofilament occupancy, and therefore axonal integrity (P < 0.01, two-way ANOVA). In WT mice, as expected from the absence of antibody and complement deposition, neurofilament occupancy was unchanged. These data thus show that complement deposition along the distal nerve in injured Glial mice does not affect axonal integrity in this acute injury paradigm, unlike injured Neuronal mice. Because of the lack of significant antibody and complement deposition in injured WT mice, and normal NoR conformation (Supplemental Figure 5A), they were not subjected to further in vivo studies.
Disruption of structural and axo-glial adhesion proteins following glial membrane targeting in vivo. We next investigated injury to the NoRs in vivo following injury. Complement only-treated controls from each genotype showed no significant abnormalities in protein distribution. Anti-GM1 antibody-directed complement-mediated injury to the glial membrane in Glial mice resulted in a major reduction in AnkB immunostaining compared with Neuronal mice ( Figure 5A, P < 0.001, two-way ANOVA). NF155 Figure 3. Differential disruption to the node of Ranvier when neuronal and glial membranes are injured selectively ex vivo. Triangularis sterni nervemuscle preparations from WT, Neuronal, and Glial mice were treated ex vivo with anti-GM1 Ab and a source of complement (injury, Inj) or anti-GM1 Ab alone (control, Con). Disruption to nodal protein (magenta) organization at the node of Ranvier (NoR) due to injury was assessed; the site of expected staining is indicated by arrowheads for each marker. Representative images demonstrate normal nodal protein localization in all control tissue and absent or abnormal staining in injury groups, which coincides with nodal complement deposition (A and C, green). (A) A pan-neurofascin (Nfasc) Ab was used to assess paranodal NF155 (closed arrowheads) and nodal NF186 (open arrowhead). (B) SC microvilli marker gliomedin (Gldn) immunostaining at NoRs was assessed compared to controls. Asterisks indicate motor nerve terminals. (C) Changes to normal (black bars) Nav1.6 labeling were observed in injured tissue from all genotypes compared with associated controls. Diamond defines statistical comparisons of absent immunostaining (white bars). (D) Perineural recordings from distal motor nerves were performed on tissue from Neuronal and Glial mice treated with anti-GM1 Ab only, a source of complement (normal human serum, NHS) only, or a combination of Ab and NHS (injured). Representative recordings from 1 mouse per treatment demonstrate that normal Na + and K + waveforms were lost when the tissue was injured. Scale bar: 5 μm. Results are represented as the mean ± SEM. n = 3/genotype/treatment: 13-36 NoRs/mouse (median = 24, pNFasc); 15-33 NoRs/mouse (median = 19, gliomedin); and 11-30 NoRs/mouse (median = 23, Nav1.6) were analyzed. *P < 0.05, **P < 0.01, ***P < 0.001 (for comparisons between normal immunostaining); ### P < 0.001 (for abnormal NF155 immunostaining in Neuronal injury group compared to WT or Glial imjury in A) compared with control by 2-way ANOVA with Tukey's post hoc test. deposits were found extensively along the glial membranes, with particular enrichment at the paranodal loops. There was a significant loss in both neurofilament staining and endogenous CFP signal along the entire distal motor nerve of injured tissue. Glial control tissue (anti-GM1 antibody only) was positive for antibody staining, negative for complement deposits, and possessed normal intact neurofilament immunostaining and cytoplasmic CFP. To demonstrate that neurofilament loss was secondary to paranodal disturbance and not the product of pSC injury and/or death, we conducted a control ex vivo experiment in WT mice exposed to anti-sulfatide antibody. As pSCs do not form myelin and thus do not express sulfatide, they were not targeted by anti-sulfatide antibody (Supplemental Figure 6), whereas the glial membranes at the paranode, formed by myelinating SCs, were injured (Supplemental Figure 4). As with anti-GM1-treated Glial tissue, loss of neurofilament staining and cytoplasmic CFP was also observed in anti-sulfatide antibody-treated mice ( Figure 7B), indicating that targeting paranodal glial membranes over 20 hours ex vivo is a sufficient time period to cause subsequent loss of axonal integrity in the distal motor nerve.
Next, we investigated axonal integrity in an extended in vivo injury model in Glial mice. At 24 hours after NHS delivery, anti-GM1 antibody and complement deposits were observed along the distal nerve glial membrane in injured diaphragms (data not shown). Unlike the acute injury, neurofilament staining intensity was reduced at the MNT (data not shown) and reached a significant reduction at the first distal NoR in injured mice ( Figure 8A, AnkB, NF155, NF186, and Caspr1 were studied again to determine whether they remain disorganized or recover their normal distribution following initial acute injury. The results demonstrate that AnkB, NF155, and Caspr1 immunostaining remained disrupted 24 hours after injury to the glial membrane ( Figure 8B). Moreover, the impaired integrity of the axon at the NoRs progressed over time. There was an enhanced reduction in the number of NF186-positive NoRs in injured Glial mice at 24 hours compared with 6 hours ( Figure 8B, P < 0.01, Student's t test). Additionally, Nav channel and AnkG clustering at the NoRs was significantly reduced in injured mice compared with control, which was not previously significant at the 6-hour time point (Figure 8, C and D, P < 0.01 and P < 0.05, Student's t test), showing that SC membrane disruption precedes axonal nodal protein disturbances. Overall, these results corroborate the ex vivo findings, confirming that injury to the glial membrane develops over time, and that a loss of paranodal integrity likely precedes and drives secondary bystander axonal loss.

Discussion
The efforts to distinguish key differentiating features of axonal (AMAN) and demyelinating (AIDP) forms of GBS have developed from widespread clinical observations emerging over the last 30 years, culminating in the recent introduction of the term nodo-paranodopathy to describe forms in which the pathology is centered on the nodal complex (1, 3). The NoR has long been recognized as a site of pathology (6), but the effect, if any, of injury to this site on bystander or secondary axonal degeneration is little understood. Herein, we have developed a model system to investigate this injury that allows us to mechanistically separate axonal immunostaining was also significantly reduced in injured Glial mice compared with all other groups ( Figure 5B, P < 0.01, two-way ANOVA). NF186-positive NoRs were reduced in both Glial and Neuronal injured mice in comparison with genotype control mice, but only reached significance in Glial mice ( Figure 5A, P < 0.05, two-way ANOVA, Supplemental Figure 5B).
NoRs with normal Caspr1 immunostaining were significantly decreased in injured Glial mice compared with Neuronal mice ( Figure 5C, P < 0.001, two-way ANOVA), where instead there was a greater number of NoRs with disrupted immunostaining; however, this did not reach significance compared to the other groups in vivo. These data indicate that the integrity of the axo-glial junction at the paranode is impaired following direct targeting of glial membranes in injured Glial mice in vivo. As observed ex vivo, there was no difference in the intensity of MBP staining along the distal internode among treatment groups (Supplemental Figure  5C), which shows that at this acute time point, injury from glial targeting is restricted to paranodal SC nodal membrane disruption and the myelin sheath remains intact.
We next studied the presence of Nav channel clustering after complement-mediated injury in vivo. Nav channel clustering was unchanged in Glial injured mice in this acute injury study; however, injured Neuronal mice had significantly fewer NoRs with normal Nav channel staining compared with all other groups ( Figure  5D, P < 0.001, two-way ANOVA). These results indicate that the localization or presence of axo-glial paranodal proteins are most severely disrupted in Glial mice acutely in vivo, while the axonal proteins at the NoRs are only mildly disrupted.
In ultrastructural analysis of tissue from in vivo experiments, NoRs from Glial control mice showed normal arrangement of paranodal loops forming septate-like junctions with the axon, neurofilaments, SC microvilli, and compact myelin ( Figure 6, A and B). In stark contrast, the paranodal loops at NoRs from injured Glial mice were highly disrupted, potentially owing to the influx of water and ions through MAC pores ( Figure 6, C and D). This arrangement was not observed in any control tissue. SC microvilli remained visible, reflecting gliomedin immunostaining results. At higher magnification, a loss of transverse bands between the node-proximal paranodal loops and axon was observed ( Figure  6D, indicated above white line). While NoRs randomly sampled from injured Neuronal mice appeared normal ( Figure 6E), the architecture of their MNTs showed disruption, with a depletion of neurofilament and synaptic vesicles, and the presence of dense vacuolated mitochondria ( Figure 6G) compared with control ( Figure 6F), and as seen previously (29). Collectively, these data indicate that in vivo autoimmune injury to the glial membrane can result in paranodal disruption and dysfunction prior to segmental demyelination; this pattern of glial injury differs substantially from the primarily axonal pathology seen in Neuronal mice.
Extended glial injury leads to secondary axonal loss of integrity. We next investigated extended injury models both ex and in vivo to look for subacute downstream events, notably loss of axonal integrity that could signify secondary axonal degeneration following selective glial insult ( Figure 7A). To this end, only Glial mice, including those crossed with axonal cytoplasmic CFP reporter mice, were studied. After 20 hours of anti-GM1 antibody and complement injury treatment ex vivo, antibody and complement lying primary and secondary axonal failure and degeneration. In ex vivo preparations this was compared to the outcomes in WT preparations where both neural membranes are targeted. While it would be ideal to model anti-GM1 antibody-mediated neuropathy in vivo in WT mice, the major handicap has been WT mouse resistance to developing GBS following passive immunization. We have attributed this to the low levels of complement activity in mice, the sequestration and depletion of circulating anti-GM1 antibody by non-neural tissues, and protection afforded by the blood-nerve barrier (16,28,30). The model described herein removes the non-neural expression of GM1, thereby raising circulating antibody to pathogenic levels, uses NHS as a potent source of heterologous complement, and studies distal motor nerves projecting out with the blood-nerve barrier. The combined application of these 3 steps allows us to observe striking clinical and morphological phenotypes in vivo. Topical application of reagents in our ex vivo paradigm reflect a severe disease pathology. Owing to the injury at the MNT resulting in synaptic dysfunction with severe respiratory distress, the in vivo experiment was restricted to 6 hours for Neuronal mice, meaning that this model represents a milder stage of injury and explains the differences in injury reported between ex and in vivo. In the future, a subacute injury paradigm with a milder clinical phenotype (31) could be applied to Neuronal mice to study nodal changes in vivo in more detail.
Delayed loss of axonal integrity following glial membrane injury. Our results show that anti-GM1 antibody binding can be selectively targeted to, and activate complement on, axonal or glial membranes from SC membrane dysfunction unambiguously, using a single monoclonal autoantibody against the prototypic GBS nodal antigen, GM1 ganglioside.
The segregation of GM1 to axonal or glial membranes (20) permits us to study selective injury to either membrane at the axo-glial interface, and thereby to distinguish mechanisms under-  findings previously observed in AMAN models and is likely driven by MAC pore formation, with intra-axonal calcium influx and calpain activation (14,16). This likely leads to the following sequence of events: calpain cleavage of axonal cytoskeletal adaptor protein AnkG; dispersion of NF186; SC microvilli gliomedin remains present; Caspr1 loss at the node-proximal paranode; consequent perturbation of the Caspr1 binding partner NF155, and glial AnkB; Nav channel cleavage and/or cluster dispersion; and conduction failure. Conversely, Glial injury to NoRs likely leads to the following sequence of events: MAC-induced paranodal loop distortion; calpain cleavage of glial paranodal scaffold protein, AnkB; disruption of paranodal NF155 through general loop distortion and loss of partial interaction with AnkB; consequent disturbance of axonal Caspr1 binding partner; mislocalization and partial loss of gliomedin; dispersion of axonal NF186; followed by Nav channel cluster and AnkG dispersion and conduction failure. Despite these differences, the effect on Nav channel clustering and function is the same, likely owing to the loss of 2 Nav channel clustering mechanisms: cytoskeletal anchors and an intact paranode.
Inducible ablation has shown that NF186 has a critical role in maintaining PNS Nav clusters at the NoR (37), and NF186 stability can underlie long-term Nav channel clustering and axon health in concert with the paranodes (38). In Neuronal mice, AnkG is likely directly cleaved, and Caspr1 is indirectly disturbed through calpain cleavage of the axonal cytoskeletal proteins 4.1B and βII spectrin that tether it to the cytoskeleton (39). These axonal paranodal cytoskeletal proteins are critical in combination with NF186 for Nav channel positioning; combined paranodal βII spectrin and NF186 conditional deficiency causes impaired AnkG and Nav channel clustering (40), reflecting our observations. Combined with Nav channel clustering loss in our Neuronal model, the presence of MAC pores in the axonal membrane drives ionic imbalance, culminating in loss of membrane potential and function, as previously observed (16). These disturbances culminate in conduction block and nodal lengthening likely through paranodal disruption from loss of node-proximal Caspr1, or membrane insertion as observed in patient autopsy tissue (6).
When the glial membranes are selectively targeted, we propose the loss of the axo-glial junction proteins at the node-proximal paranodal loops begins with distortion, swelling, and loss of paranodal loop integrity through the influx of water and ions via MAC pores, as observed in ultrastructural images. Loss of the calpain substrate AnkB immunostaining is an indicator of MAC-poreinduced paranodal loop disruption. While AnkB loss is not entirely causal, it could contribute to the disruption of NF155 through their partial interaction. Normally, NF186 maintains Nav channel clusters through interactions with SC microvilli extracellular matrix proteins, including gliomedin (17). Therefore, we propose that the combined disturbance of SC gliomedin and the axo-glial CAMs leads to Nav channel dispersion, which underlies, in part, the ensuing conduction failure. Chang et al. (24) found no change to NF155 and Caspr1 following conditional ablation of AnkB; however, they propose a compensation by substitution with other glial cytoskeletal spectrins and actins. These proteins would presumably also be cleaved in our acute injury setting, thereby preventing this substitution from becoming functional and exacerbating the disruption. It is true that mice deficient in axo-glial CAMs have at distal NoRs and MNTs. Loss of neurofilament immunostaining, indicating a loss of axonal integrity, follows acutely in Neuronal but not Glial tissue, suggesting direct targeting of the axonal membrane, as previously observed in other models of AMAN (16,21). However, most significantly, we find that following early paranodal disruption and nodal axonal protein disturbances in Glial mice, axonal integrity is compromised following extended glial injury, thus indicating that indirect axonopathic mechanisms must account for subsequent axonal loss in the extended model. These results, first seen ex vivo, were then recapitulated in vivo in an extended Glial injury model. Notably, our results also confirm that despite pSC injury in our Glial paradigm, acute injury to these specialized glial cells does not account for the acute loss of axonal integrity seen in the distal motor nerve, aligning with previous observations (15,32). A diverse array of mechanisms have been proposed to contribute to secondary bystander axonal degeneration including, but not limited to, demyelination, mitochondrial dysfunction, and disrupted trophic support (33). Because of the proximity of axonal and glial membranes in the nodal complex, it is possible that cytotoxic mediators are released from the injured glial membrane that could indirectly cause axonal damage. There is no evidence that MAC pores can be shed and reinserted intact into another adjacent membrane (34); however, MAC formation is a highly toxic injury resulting in a wide range of uncontrolled bidirectional solute fluxes. Based on our data, herein we will discuss the contribution of paranodal axo-glial and nodal protein disturbances to axonal degeneration.
Selective membrane targeting leads to loss of axo-glial integrity and function through different pathways. To account for the acute cascade of molecular events leading to collapse of the axo-glial interface at the paranodal junction, our results suggest a differential mislocalization of proteins dependent on 4 factors: (a) local anchoring of CAMs through interactions with cytoskeletal adaptor proteins; (b) axo-glial stabilization and tethering through axo-glial trans interactions; (c) calpain substrate status of intracellular cytoskeletal adaptor protein domains; and (d) distortion of the paranodal loops and loss of integrity. In this study we interpret a loss of immunostaining as the absence or dispersion of a protein. For the abundant cytoskeletal protein neurofilament, absence is an indicator of compromised axonal integrity, which could progress to axonal loss/degeneration. Loss or disruption of nodal protein immunostaining may be through direct calpain cleavage of known substrates (AnkB, AnkG, and Nav) (35,36), or indirectly through cleavage of known tethering proteins causing dispersion, detailed below. Injury to Neuronal tissue reflects Figure 7. Extended ex vivo injury selectively targeting glial membranes results in secondary axonal degeneration. Triangularis sterni nerve-muscle preparations from Glial (A) and WT (B) mice were treated ex vivo with anti-GM1 Ab or anti-sulfatide Ab, respectively, and a source of complement (injury, Inj) or Ab alone (control, Con) for 20 hours. (A) Anti-GM1 Ab (orange) and complement (green) deposition along the distal motor nerve was strongly enriched at the paranodes (arrowheads) in injured compared with control tissue. Loss of axonal integrity along the distal nerve was monitored by presence of neurofilament H immunostaining (NFH, magenta) and cytosolic CFP (blue). (B) The experiment was repeated in WT mice using an anti-sulfatide Ab; the results reflect those reported in A. Scale bars: 10 μm (A) and 20 μm (B). Results are represented as the mean ± SEM. n = 3/treatment: 25-54 NoRs/mouse (median = 39) were analyzed. **P < 0.01, ***P < 0.001 compared with control by 1-tailed Student's t test.
ary loss of axonal integrity, manifest as a porous membrane. In a recent model of experimental autoimmune encephalomyelitis, calcium reporter mice were used to monitor calcium entry into axons traveling through inflammatory lesions (45). The authors discovered that extracellular calcium ions were entering through nanoruptures and levels could ultimately predict axonal fate. They predict these nanoruptures could be caused by toxic mediators or mechanical forces, the latter neatly aligning with the proposed disturbances to the regular pattern of the actin/spectrin axonal cytoskeleton in our extended model. Cytoplasmic CFP loss in our extended model strongly indicates a loss of axonal membrane integrity. Furthermore, neurofilament is a known calpain substrate and its loss in the same paradigm points toward calcium entry and calpain activation as the final executioner.
Clinical significance. A recent clinical development has been the introduction of the term nodo-paranodopathy as an explanatory concept focused on this site that in part escapes from the binary AMAN/AIDP classification (3). In parallel with these advances in clinical and electrophysiological observations, little effort has been applied to developing animal models that might support them experimentally. Exploring the mechanisms underlying the nodo-paranodopathy concept will undoubtedly benefit from the generation of the Neuronal and Glial transgenic mice and other relevant models (46).
The pathological situation of acute anti-GM1 antibodyassociated GBS in humans, be that AMAN, AIDP, nodo-paranodopathy, or a hybrid form, would likely resemble the situation described herein in the WT mouse. There are, however, many unknown factors that make this assumption of human-murine parity overly simplistic. Firstly, the distribution and concentration of GM1 in human axonal and glial membranes relative to the mouse is unknown. Secondly, both human and mouse polyclonal and monoclonal anti-GM1 antibodies differ considerably in their ability to bind to GM1 in different membrane environments. Cis interactions between GM1 and other molecules in the plane of the plasma membrane shape the topographical presentation of GM1 epitopes for antibody binding that may be enhancing or inhibitory (47). Thus, certain anti-GM1 antibodies may preferentially favor binding to glial GM1 versus axonal GM1, or vice versa. In the context of anti-GM1 antibody-associated GBS this could affect the AMAN/AIDP relative balance. Notwithstanding these considerations, it seems likely that anti-GM1 antibody-associated GBS in humans would be driven by a summation of direct and indirect axo-glial injury mechanisms that may evolve at different rates over time according to antibody titer, specificity, and host-intrinsic inflammatory and axonal vulnerability factors.
The specialized nature of the distal nerve site studied over a relatively brief time period is a constraint to translating our current model into a human homolog, although this distal site is certainly affected in human GBS (2,3). Looking forward, further extending this model over time in vivo, subject to animal welfare constraints, will allow us to elucidate the mechanism of secondary bystander axonal degeneration, or reversible conduction failure, and thereby add mechanistic insights into the likely effects of novel treatments in humans. Characterizing SC membrane and myelin injury will allow us to study whether these early disease stage findings exemplify reversible axonal damage or can lead to more permanent minimal disturbances in Nav channel clustering (41)(42)(43), likely due to the presence of NF186. But here, similar to Neuronal injury, a disturbance of both NF186 and paranodal integrity causes Nav channel cluster loss in Glial mice. Additionally, paranodal loop disruption also likely leads to a leakage of driving current due to disruption of axo-glial junctional scaffolds causing conduction block, as previously described (3).
Conduction failure can be transitory and reversible, or irreversible if the axons ultimately undergo transection with subsequent Wallerian degeneration. The tipping point -a concept also referred to as the metastable state -that dictates the capability of injured axons to recover through local repair in favor of complete axonal transection is likely to be mechanistically multifactorial (44). Therefore, distinguishing the outcome of immune attack on axonal or glial membranes that have similar electrophysiological outcomes has clear implications for therapeutic strategies and recovery. Whether, over time, our current models would proceed to segmental demyelination, either directly, or as a result of axonal transection with Wallerian degeneration, or be restricted to focal membrane injury around the nodal region, remains to be determined. However, from existing evidence it is reasonable to predict that both segmental demyelination and/or Wallerian degeneration would be demonstrated in a prolonged-injury model.
Loss of axo-glial integrity and secondary bystander-axon degeneration. It has been shown that the axo-glial CAMs NF155, Caspr1, and contactin-1 assemble a specialized axonal cytoskeleton at the paranode (protein 4.1B, αII and βII spectrin), which promotes nodal domain formation and anchoring to the actin cytoskeleton (17). Additionally, it has been proposed that axonal actins and spectrins form a periodic cytoskeleton with submembranous rings, and coalignment of submembranous paranodal cytoskeletal proteins is critical to arranging the paranodal CAMs (17). Our Kv1.1 data from the acute injury paradigm suggests the localization of these cytoskeletal proteins remains intact, as there is no invasion into the paranodal domain. However, the delayed consequence of a distortion to this regimented arrangement by a disturbance in glial CAM localization could be second- Figure 8. Extended in vivo injury selectively targeting glial membrane results in secondary axonal degeneration. Glial mice were dosed i.p. with 50 mg/kg anti-GM1 Ab followed 16 hours later with 30 μL/g normal human serum (NHS) (injury, Inj) or NHS only (control, Con). The experiment was terminated 24 hours after NHS delivery. The site of expected nodal protein immunostaining is indicated by arrowheads. (A) At this time point there was loss of neurofilament H staining (NFH, orange) at the motor nerve terminal (MNT, asterisk) and the staining intensity was significantly reduced at the first distal node of Ranvier (NoR). (B) Normal ankyrin B (AnkB), NF155, NF186, and Caspr1 (magenta) immunostaining was assessed at distal paranodes after injury compared to control. (C) There was a further reduction in distal NoRs with normal voltage-gated sodium (Nav) channel staining (magenta) in injured Glial mice compared with control at this extended time point. (D) Additionally, the Nav channel-tethering protein AnkG was notably absent. Scale bar: 5 μm. Results are represented as the mean ± SEM. n = 3/genotype/ treatment: 5-15 NoRs/mouse (median = 11, NFH intensity); 4-25 NoRs/ mouse (median = 18, panNFasc); 9-23 NoRs/mouse (median = 12, Nav); and 10-28 NoRs/mouse (median = 18, AnkG) were analyzed. *P < 0.05, **P < 0.01, ***P < 0.001 by 2-tailed Student's t test (A and B) or 2-way ANOVA with Tukey's post hoc test (C).
[NGS] + 0.5% Triton X-100 in PBS). For diaphragm sections from in vivo experiments, antibodies were prepared in PBS plus 3% NGS and 0.1% Triton X-100. Secondary antibodies were diluted in PBS plus 1% NGS. NHS was collected from a single donor, rapidly frozen, and stored in multiple aliquots at -70°C to preserve complement activity.

Mice
WT mice (Harlan) and 2 previously described transgenic strains, GalNAc-T -/--Tg(neuronal) (21) and GalNAc-T -/--Tg(glial) (20), all on a C57BL/6 253 background backcrossed for 7 generations, were used. Briefly, GalNAc-T -/--Tg(neuronal) mice express the full-length cDNA encoding GalNAc-T under the control of the human Thy1.2 promoter (restricted to mature neurons), and GalNAc-T -/--Tg(glial) mice express the full-length cDNA encoding GalNAc-T under the control of the mouse Plp promoter (restricted to glia, including myelinating and nonmyelinating SCs). Additionally, mice were crossed with B6.Cg-Tg mice that endogenously express cytosolic cyan fluorescent protein (CFP) in their axons (21). Mice from the GalNAc-T -/--Tg(neuronal) or GalNAc-T -/--Tg(glial) background exhibit age-dependent neurodegeneration (20); therefore, we elected to use 4-to 6-week-old mice, both male and female, that have no identifiable phenotype (12). The number of mice per treatment are reported per experiment. Mice were maintained under a 12-hour light/dark cycle in controlled temperature and humidity with ad libitum access to food and water. For each study, mice were euthanized by increasing-CO 2 inhalation.

Ex vivo injury models
Nerve integrity. Ex vivo TS nerve-muscle preparations from WT, Neuronal, and Glial mice were used for each study as described previously with some modifications (21). Each mouse provides 2 TS: one half was used as an uninjured control (anti-GM1 antibody only) and the second half as the injured preparation (anti-GM1 antibody plus NHS). For the acute injury model, all TS were incubated for 4 hours at 32°C in Ringer's with 100 μg/mL anti-GM1 antibody and additionally 40% NHS for the injury half. For the extended injury model, TS were incubated for 20 hours at room temperature. The TS were washed 3 times in Ringer's, followed by a 1-hour incubation with anti-MAC antibody at 4°C. For the EhtD-2 homodimer live stain (2 mM), the label was applied for 30 minutes at room temperature prior to fixation. The tissue was washed 3 times in Ringer's followed by fixation with 4% paraformaldehyde (PFA) at 4°C. Washes with PBS, 0.1 M glycine, and PBS followed. TS were transferred to 100% ethanol for 10 minutes at -20°C, and then thoroughly washed in PBS. Tissue was incubated overnight at 4°C in blocking solution plus mouse anti-MAC antibody, rat-anti MBP antibody, and α-bungarotoxin with one of the neural markers under investigation. TS were rinsed 3 times in PBS followed by a 2-hour incubation with secondary antibodies at room temperature in the dark. TS received final washes in PBS and were mounted in Citifluor mounting medium (Citifluor Products). For each neural marker, n = 3/4 mice per genotype.
Perineural recordings. Perineural recordings were made from TS nerve-muscle preparations set up as described previously (16,48), with some modifications. Briefly, preparations were set up for electrophysiological recordings in Tyrode's solution following a 3-hour incubation of the TS with 100 μg/mL anti-GM1 antibody plus 40% NHS in Tyrode's at 32°C. Antibody only and NHS only controls were used. For each of the 3 treatments, n = 4 mice/genotype. axonal degeneration; shifting the metastable state in favor of local repair rather than axonal transection (44) could have major effects on clinical outcome. Indeed, we observed some mild changes to myelin (vesiculation) and EthD-2-positive internodal nuclei that show signs of early SC body injury, the likely forerunner of segmental demyelination. SC nodal membrane disruption is just one of many mechanisms that occur in response to anti-GM1 antibody-mediated injury and comprise the path to degeneration. While certainly affected in human GBS, the distal nerve is a site with specialized features that must be taken into consideration. Distal axons are small in diameter and vulnerable to damage, and the pSCs are capable of promoting rapid repair and regeneration after axonal injury (31). Exploring these factors in models remains highly relevant to understanding the highly complex axo-glial interplay at the different anatomical sites affected in human GBS.
For ex vivo injury preparations, tissue was maintained alive in oxygenated (95% O 2 /5% CO 2 ) Ringer's solution (116 mM NaCl, 4.5 mM KCl, 1 mM MgCl 2 , 2 mM CaCl 2 , 1 mM NaH 2 PO 4 , 23 mM NaHCO 3 , 11 mM glucose, pH 7.4). Ex vivo electrophysiological preparations were maintained alive in oxygenated Tyrode's solution (137 mM NaCl, 2.7 mM KCl, 0.5 mM MgCl 2 , 2.5 mM CaCl 2 , 0.4 mM KH 2 PO 4 , 11.9 mM NaHCO 3 , 10 mM HEPES, 11 mM glucose, pH 7.4). For ex vivo preparations, all primary antibodies were applied in a blocking solution (3% normal goat serum gap/weakened staining; (c) Absent. If only 1 paranodal protein domain remained then this was categorized as absent. FIJI software (Fiji) was used to quantify intensity of neurofascin and MBP across the NoRs. Fluorescence intensity centered on the nodal gap was plotted for each marker. For all ex vivo analysis of the effect of complement on neural structures, only NoRs positive for complement deposition in the injured group were included for quantification.

Experimental design
For in vivo experiments, mice were randomly allocated to treatment groups using a random number generator. Animals were excluded from the injury treatment group if the serum was absent for AGAb or complement assessed by ELISA and topical complement assay, respectively. A power analysis was performed using G*Power software v3.0.10 to determine group size; n = 3 (ex vivo) and n = 4 (in vivo) was selected for each treatment group. Effect size for behavioral output was based on previous experiments, calculations made on the basis of 80% power, and 0.05 significance criteria. All tissue was coded to prevent researcher bias.

Statistics
The number of independent animals are described in the Methods and indicated in the figure legends. Statistical differences among genotypes were determined by Student's t test, 2-way ANOVA for multiple factors, or repeated measures 2-way ANOVA followed by a Fisher's, Tukey's, or Bonferonni's post hoc test for multiple comparisons using GraphPad Prism 7 software. Parametric testing was used, and differences were considered significant when P values were less than 0.05. Data were plotted as the mean ± SEM using dot plots or contingency graphs.

Study approval
All experiments using mice conformed to University of Glasgow institutional guidelines and complied with relevant guidelines on the care and use of animals outlined in the revised Animals (Scientific Procedures) and were performed in accordance with a license approved and granted by the United Kingdom Home Office (POC6B3485).

Author contributions
RM and HJW designed the research and wrote the manuscript. RM, CIC, JAB, MEC, DY, and EGR conducted experiments and acquired and analyzed the data. CLC and SR performed pilot studies. RM, CIC, JAB, MEC, CLC, SR, EGR, and HJW reviewed the manuscript.


Metabolomic predictors of phenotypic traits can replace and complement measured clinical variables in population-scale expression profiling studies

Population-scale expression profiling studies can provide valuable insights into biological and disease-underlying mechanisms. The availability of phenotypic traits is essential for studying clinical effects. Therefore, missing, incomplete, or inaccurate phenotypic information can make analyses challenging and prevent RNA-seq or other omics data to be reused. A possible solution are predictors that infer clinical or behavioral phenotypic traits from molecular data. While such predictors have been developed based on different omics data types and are being applied in various studies, metabolomics-based surrogates are less commonly used than predictors based on DNA methylation profiles.In this study, we inferred 17 traits, including diabetes status and exposure to lipid medication, using previously trained metabolomic predictors. We evaluated whether these metabolomic surrogates can be used as an alternative to reported information for studying the respective phenotypes using expression profiling data of four population cohorts. For the majority of the 17 traits, the metabolomic surrogates performed similarly to the reported phenotypes in terms of effect sizes, number of significant associations, replication rates, and significantly enriched pathways.The application of metabolomics-derived surrogate outcomes opens new possibilities for reuse of multi-omics data sets. In studies where availability of clinical metadata is limited, missing or incomplete information can be complemented by these surrogates, thereby increasing the size of available data sets. Additionally, the availability of such surrogates could be used to correct for potential biological confounding. In the future, it would be interesting to further investigate the use of molecular predictors across different omics types and cohorts. Supplementary Information The online version contains supplementary material available at (10.1186/s12864-022-08771-7).

Background
Genome-wide association studies (GWAS) have proven to be valuable in uncovering links between genes and a wide range of phenotypic traits. Such findings have led to the discovery of new disease-related biomarkers and are often the basis for gaining a better understanding of biological processes or disease-mechanisms. Since the introduction of the first GWAS powered by the availability of genome-wide single-nucleotide polymorphism (SNP) profiling, numerous studies have identified thousands of SNP-trait associations [1]. Technological advancements allowing high-throughput profiling of other molecular features, such as transcripts, or DNA methylation sites, also enabled population-scale studies of transcriptomics, epigenomics, and other omics data types.
Such studies are susceptible to confounding by biological and technical factors that can influence omics profiles and phenotypic traits of interest. However, measured values to correct for such confounding are often not available. As a solution, differences in cell type composition are commonly accounted for using information contained in the DNA methylation profiles themselves, by either reference-based imputation [2] or reference-free methods such as surrogate variable analysis (methods reviewed in [3] and [4]). Other well-known examples of inferring values for possible confounding factors from DNA methylation profiles include sex [5] and smoking status prediction. Bollepalli et al. [6] trained a smoking status classifier using multinomial LASSO regression. Machine learning approaches have also been applied to other omics data types to predict environmental exposures [7].
The value of such predictors is not only evident when complementing missing data to account for technical or biological confounding, but also for using them as outcome variables. These molecular surrogates can be used in association studies in order to link molecular features to clinical phenotypes or exposures. Since identified associations in ome-wide association studies often have only moderate effect sizes, a common approach to detect relevant features are cross-cohort meta-analyses [8]. However, the applicability of meta-analyses can be limited by availability of the respective outcomes of interest. Specific clinical, environmental, or phenotypic traits might not be recorded in every cohort, or the data collection might be based on different protocols, making the reported values for these traits not directly comparable.
As more and more multi-omics data sets become available, it becomes possible to make use of molecular predictors to infer phenotypic traits from specific omics layers. Blood is a key specimen in clinical diagnostics reflecting on the health state of an individual. While blood metabolomics methods partially overlap with classical clinical diagnostics methods, they can measure a wide range of metabolites and have the potential to play an important role in personalized medicine approaches [9]. Recently, Bizzarri et al. [10] trained predictors on proton NMR-based metabolomics (Nightingale Health) data. The authors applied logistic regression using elastic net regularization to train models for various clinical variables, including physiological measures, environmental exposures, and clinical endpoints. They demonstrated the use of these surrogates in metabolome-wide association studies to complement missing clinical data and correct for confounding. They further showed that metabolomic surrogates can help explore independent risk factors of all-cause mortality in older individuals [10].
We here propose the use of metabolomics-derived surrogates in analyses of other omics layers, not only as covariates to account for confounding factors, but also as outcome variables. We investigated whether values derived from molecular predictors represent a viable alternative to measured or reported clinical or phenotypic traits to serve as outcome variables in population-scale gene expression profile association studies. To this end, we applied 17 metabolomic predictors to metabolomics data from four large population cohorts inferring values for phenotypic, exposure, and clinical traits. We performed association studies on corresponding RNA-seq data sets employing either reported/measured or inferred values as outcome variables, and systematically compared the respective results of these analyses. For five of the outcomes, where reported values were not available, we evaluated the performance of the metabolomic surrogates based on results reported in literature.

Results
In this study, we performed expression profiling studies to evaluate the performance of 17 surrogate outcomes that are based on molecular predictors. Metabolomics data used to infer surrogate outcome values and RNAseq data used for the evaluation are part of multi-omics data sets of four large Dutch population cohorts: LifeLines (LL) [11], Leiden Longevity Study (LLS) [12], Netherlands Twin Register (NTR) [13], and Rotterdam Study (RS) [14]. An overview of the study workflow is shown in Fig. 1. For the different outcomes (Table 1), we compared expression profile association study results based on metabolomic surrogate outcomes to those based on reported outcomes whenever possible. Additionally, results were compared to literature findings.

Metabolomic surrogate outcomes
We inferred values for clinical variables by applying molecular predictors to metabolomics data. Recently reported metabolomic predictors trained on up to 22 Dutch population cohorts [10] were applied to infer values for outcome variables (Table 1). All 17 metabolomic predictors had been shown to perform accurately with mean Fig. 1 Overview of study workflow. Gene-wise models are fitted for various outcome variables based on reported information or metabolomic surrogates, respectively AUC values > 0.7 in a 5-fold cross-validation approach [10]. In order to avoid emphasis on clinical extremes, the metabolomic predictors trained by Bizzarri et al. are based on binary representations using clinical thresholds for continuous variables. The values returned by the predictors are continuous posterior probability scores for belonging to one of the two groups. In most cases, this is the clinical risk group. We here used these predicted values as surrogate outcomes and compared their use in expression profiling studies to reported or measured outcome values.

Expression profile association studies
In the next step, the 17 metabolomic surrogates were used as outcomes of expression profile association studies. In addition to analyses employing these surrogate outcomes, analyses using measured or reported outcome values were performed. For five outcomes that had limited availability of reported values (eGFR, diabetes, metabolic syndrome, blood pressure lowering medication, and alcohol consumption) only surrogate outcomes were used. In each linear regression model, known biological (age, sex) and technical (flow cell number, white blood cell composition) confounding factors were included (formulas available in Additional file 1).
For an initial assessment of the performance of each model, we compared the numbers of significant associations and the effect sizes between outcomes and outcome variable types. Additionally, test statistic (t-statistic) bias and inflation were estimated as parameters (mean and standard deviation) of the empirical null distribution using a Bayesian method implemented in the R package bacon [15] (Fig. 2). Numbers of identified significant associations ( Fig. 2A) varied strongly across outcomes. Highest numbers were found for the outcomes triglycerides, metabolic syndrome, and white blood cells. For several outcomes, including eGFR, LDL cholesterol, total cholesterol, and alcohol consumption, no or only few significant gene-trait associations were found. For outcomes where association study results based on surrogate outcomes could be compared to results for reported variables, the numbers of identified significant associ-ations averaged across all cohorts were higher for the metabolomic surrogate outcomes in two cases, and lower in 10 cases. However, the variation across the four cohorts was generally higher than the difference between models employing either reported or surrogate outcome variables. Similarly, high variation across cohorts was observed for the other parameters assessed to evaluate the performance of the models. Absolute effect size averaged across all genes (Fig. 2B) were generally small, with the highest values observed for the outcomes triglycerides and sex. In 10 cases, the mean absolute effect size averaged across cohorts was lower when using metabolomic surrogate outcomes instead of reported variables; in two cases, it was higher. We observed relatively low test statistic bias (Fig. 2C) across all outcomes and types of outcome variables. The bias, i.e., the deviation of the  empirical null distribution's mean from zero, averaged across cohorts decreased in four cases, increased in three cases, and remained similar in five cases when employing metabolomic surrogate outcomes instead of reported variables. Bias in the RS cohort was often higher than in the other cohorts. This may be explained with the differences in population structure. The RS cohort has a higher average age than the other three cohorts [16], indicating higher bias for the studied clinical variables in older populations. Inflation (deviation of the empirical null distribution's standard deviation from one, Fig. 2D) was highest for the outcome sex. In most cases, inflation averaged across cohorts remained stable when using different outcome variable types. For the outcome total cholesterol, it slightly decreased when using metabolomic surrogate outcomes instead of reported variables; for the outcomes lipid medication and hsCRP, it slightly increased. Since the number of significant associations and average effect sizes do not allow drawing conclusions about the similarity of the association study results employing different types of variables as outcome, we next performed pairwise comparisons of models with different types of outcome variables, i.e., reported vs. surrogate. Figure 3 shows the correlation of regression coefficients from gene-wise fitted linear models between two different types of outcome variables. Correlation coefficients were generally high for surrogate outcomes based on bestperforming metabolomic predictors. For outcomes based on predictors with reported AUC > 0.9 (triglycerides, LDL cholesterol, total cholesterol, HDL cholesterol, and sex) [10], the correlation coefficients averaged across cohorts ranged between 0.71 and 0.96. We observed a modest trend for decreasing correlations with decreasing performance of predictors. However, there were exceptions, with sex having the highest correlation values, although the predictor's AUC was reported to be lower than those for triglycerides or cholesterol. Lowest similarity of results with average absolute Pearson r < 0.5 were observed for the outcomes age, white blood cells, smoking, and hemoglobin, the latter having the lowest correlation values. We often observed that correlations were lower for the NTR cohort. This could be explained by a technical difference in the metabolomic profiles, with NTR missing glutamine [10].

Meta-analyses and replication studies
When comparing two alternative outcome variables, a lower or higher number of found significant associa- tions does not necessarily imply that results are better or worse, since the values alone do not indicate if this is due to a reduction or increase of false positive (noise) or true positive findings, respectively. We observed that the expression profile association study results differ when surrogate values differ from reported values. However, we do not know which set of results is correct, as reported values could contain inaccuracies. Under the assumption that true positive findings, as opposed to false positive results, can be replicated in different cohorts (validating the results), we performed replication studies to determine which outcome variable type is more consistently reflected in the RNA-seq data. We performed leave-onecohort-out meta-analyses and replication studies for all comparisons (except for hsCRP where only two cohorts were available) using the approach described by van Rooij et al. [16]. For each comparison, four meta-analyses were performed leaving one cohort out each time, and using the left out cohort for a replication analysis. Figure 4 shows the numbers of significant associations found in each meta-analysis (number of meta-analyzed genes) and the respective percentage of replicated genes. Overall, we did not find substantial differences in the numbers of meta-analyzed genes (Fig. 4A) except for the outcomes lipid medication, (high) age, and current smoking. While more genes were meta-analyzed when using the metabolomic surrogate for lipid medication, the reported variable yielded more meta-analyzed genes for age and smoking. For the latter two outcomes, results based on metabolomic surrogates could not be replicated (Fig. 4B) while on average 30-38% of results based on reported outcomes could be replicated. For these two outcomes we had also observed the highest differences between number of significations associations ( Fig. 2A). For other outcome variables, the percentage of replicated genes was quite similar between outcome variable types, but the cohort which was left out for the meta-analysis had a strong impact on the results. Highest average replication rates were observed for triglycerides with 69% for the reported and 67% for the surrogate outcome. For a number of outcomes, associations could hardly be replicated: LDL-associated cholesterol, hemoglobin, and alcohol consumption. This is in line with the fact that almost no significant associations were found for these outcomes (except for reported hemoglobin) in the individual cohorts (Fig. 2).

Gene set enrichment analysis
To arrive at a biological interpretation of the association study results, we performed gene-set enrichment analyses (GSEA) using pathways from the Reactome database. GSEA was applied to both individual cohort results and results from a meta-analysis of all four cohorts (Fig. 5). For direct comparisons of reported variables and metabolomic surrogates, we observed a highest overlap of significantly enriched pathways for HDL cholesterol and triglycerides in all cohorts, with 70-84% (meta-analysis 76%) and 65-80% (meta-analysis 77%) of significantly enriched pathways found by both reported and surrogate outcome, respectively. The overlap for the outcomes total cholesterol, lipid medication, hsCRP, BMI/obesity, and white blood cells was more variable across cohorts. Metaanalyzed results had an overlap of 38-69% and the order of significantly enriched pathways was highly comparable (see Additional file 2). The results for high age, current smoking, hemoglobin, and LDL cholesterol demonstrated a lower overlap than other outcomes and showed higher variation across the four cohorts compared to other outcomes. This is partially in line with the comparison of gene-wise linear models from association studies (Fig. 3), which showed that results for the outcome hemoglobin based on reported and inferred values were not correlated, and high age and current smoking were only moderately correlated. Since hardly any significant associations were found for hemoglobin (see Fig. 2A), the observed signal for this outcome was generally very low in the studied cohorts, independent of the type of outcome variable. It is surprising to observe that almost no significantly enriched pathways were observed for the outcome sex, even though many gene-trait associations were found and could be replicated in the meta-analysis and replication approach.
In order to evaluate the performance of surrogates for which results could not be compared to results based on reported outcomes, we compared significantly enriched pathways (see Additional file 2) to the literature. The top-ranked enriched pathways for low eGFR are related to translation. Since low eGFR is an indicator of kidney disease, this is in line with studies reporting increased translational activity to several kidney diseases [17,18]. For diabetes and metabolic syndrome, which is a risk factor for diabetes, 16 out of the top 20 significantly enriched pathways were found for both outcomes. These include pathways related to translation, signaling, infection, and amino acid deficiency and metabolism. This is in agreement with previously reported results [19,20]. For alcohol consumption, although almost no significant associations were found in individual analyses, several significant gene-outcome associations were found when meta-analyzing multiple cohorts (compare Figs. 2 and 4). Top-ranked positively enriched Reactome pathways from gene-set enrichment analysis (Additional file 2), including, e.g., innate immune system, signal transduction, and infectious disease, have been linked previously to chronic alcohol drinking [21].

Discussion
While certain omics predictors especially based on DNA methylation profiles [2][3][4][5] are regularly applied in (multi) omics data analyses, metabolomics-based predictors are not commonly used in the analysis of other omics data types. Previously, Bizzarri et al. had shown that metabolomic surrogates can be used to correct for confounding in metabolome-wide association studies [10]. In this study, we investigated the use of these surrogates as outcome variables in the analysis of another omics level. We systematically compared results of populationscale gene expression profile association studies against outcome variables that were either reported or inferred by molecular predictors. The results generally showed good agreement (Fig. 2). Most similar association study results across all assessment parameters are those for the outcomes triglycerides and HDL cholesterol. Many significant gene-trait associations were found, of which many could be replicated, and the majority of significantly enriched pathways were found by reported and surrogate outcomes. Regression coefficients of the models including the reported and surrogate outcomes were strongly correlated, and the majority of pathways found by GSEA were obtained by both outcome types. The top-ranked pathways positively enriched for both high triglycerides and low HDL-cholesterol include "GTP hydrolysis and joining of the 60S ribosomal subunit", "L13a-mediated translational silencing of Ceruloplasmin expression" and "Formation of a pool of free 40S subunits" which participate in the Eukaryotic translation initiation [22] and were previously shown to be enriched in a high-cholesterol and high-fat diet induced hypercholesterolemic rat model [23]. The metabolomic predictors for these outcomes are directly related to metabolic markers measured on the Nightingale platform and had shown high performance (AUC ≥ 0.95) [10]. It is expected that results based on predicted outcomes will depend on the accuracy of the prediction. Accordingly, we observed a slightly lower correlation between reported outcomes and metabolomic surrogates for molecular predictors that were known to have lower accuracy (Fig. 3). In the association studies for lipid medication and BMI/obesity, the molecular predictors yielded even more significant gene-trait associations than the reported outcomes ( Fig. 2A) and a higher number of significantly enriched pathways were obtained when using the metabolomic surrogates (Fig. 5). For lipid medication, this may be related to inaccurate recording of this trait in the questionnaires used. For BMI, this may be explained by the more direct capturing of metabolic processes that are associated with obesity as a combined measure of BMI and waist circumference, because BMI alone is not a perfect indicator of metabolic health [24]. Similarity of results based on surrogate and reported value for sex was high across all assessment parameters, but GSEA did not yield significantly enriched pathways in most cases. It is possible that the genes significantly associated with sex belong to too many pathways and/or that some genes within a pathway have a positive association while others have a negative association resulting in a failure to identify positively or negatively enriched pathways. For the outcomes total cholesterol and hsCRP, regression coefficients were moderately correlated and GSEA results were very similar. For the outcome white blood cells overlap of significantly enriched pathways in GSEA was smaller. However, the order of top-ranked pathways was similar (see Additional file 2). Additionally, the application of metabolomic surrogates for outcomes that were not reported in the data and the comparison of association studies and GSEA results with literature show that these surrogate outcomes allow transferring information from one data set (the training data) to another, thereby facilitating to study phenotypes or exposures in data sets which would otherwise not be possible. These results suggest that metabolomic surrogates are a useful tool to complement phenotypic information of multi-omics data sets and enable analyses of clinical outcomes even when they are not reported. This is especially useful when reanalyzing existing data sets. Even though this study comprises four different large population cohorts and 17 metabolomic surrogates, it will be interesting to investigate in the future whether similar results can be observed in other cohorts, for other clinical variables, or for other omics data types.
For both outcome variable types, only few gene-trait associations were significantly associated with the outcome LDL cholesterol. Similarly, low numbers of significant associations were observed for hemoglobin for the surrogate outcome. In contrast to that, the reported outcome yielded more associations which, however, could not be replicated. For high age and current smoking slightly fewer associations were found when using the surrogate outcomes. They also performed worse in the metaanalyses and replication studies compared to reported outcome values. Here, differences between results based on different outcome variable types are reflected in lower correlations of regression coefficients of the gene-wise models and in a smaller overlap in enriched pathways from GSEA. Possible reasons for the differences between surrogates and reported values are lower performance of the metabolomic predictors for these outcomes, or a lower biological signal for the respective clinical outcomes and thus increased noise in the studied data. It is known that aging is reflected in transcriptomics data [25], but the metabolomic predictor for high age trained on binarized data (≥ 65 y.o.) might not be an ideal surrogate to study this. Alternatively, a metabolomics-based biological age predictor based on continuous data [26] might perform better. Differences between GSEA results of different outcome variable types could also arise from different biological information captured by the metabolomic surrogates and by reported or measured values. This phenomenon is known from epigenetic clocks whose age predictions can differ from chronological age, and different clocks can reflect different aspects of biological age [27]). While many of the top-ranked pathways (Additional file 2) for smoking were found by both outcome variable types, some pathways were solely found by using either reported smoking status ("smoking_current") or metabolomic surrogate("s_current_smoking"). Several pathways related to translation initiation ("Formation of the ternary complex, and subsequently, the 43S complex", "Translation initiation complex formation", "Ribosomal scanning and start codon recognition") were only significantly enriched when using the reported variable as outcome. Translation of mRNA is known to be dysregulated in cancers [28]. Pathways only enriched when using the metabolomic surrogate include "Platelet activation, signaling and aggregation" and "Hemostasis". Increased platelet aggregation has been reported in smokers [29] and platelet-dependent thrombin levels were shown to be increased in smokers and following smoking [30]. This possibly indicates that the reported smoking behavior captures effects of longterm exposure to smoking better, while the metabolomic surrogate captures effects of acute smoking. It would be interesting to further investigate which aspects of the clinical phenotypes are captured by the metabolomic surrogates. This requires additional phenotypic information. To understand which aspect of smoking behavior is reflected in the omics data current smoking status alone might not be sufficient. More information including pack years and years since smoking cessation could help better understand the information captured by the predictors. It is also possible that different omics types capture different effects better, e.g., short-term and long-term effects. In that case, combining reported outcome variables, and/or molecular surrogates from different omics layers could be very useful, not only to study the effect of a certain exposure, but also to better adjust for confounding factors.

Conclusions
In our systematic comparison of expression profiling results using either reported variables as outcome or surrogate outcomes inferred from metabolomics profiles, we demonstrated that in many cases metabolomic sur-rogates yield similar results as reported variables. We showed that the availability of these surrogate outcomes extends the possibilities of studying various clinical outcomes in population cohorts. It can enable the reuse of existing multi-omics data with limited reported clinical (meta)data. This allows for inclusion of more cohorts in meta-analyses, even when outcomes of interest were not reported for all cohorts. This approach also increases possibilities to study clinical outcomes by allowing to infer important confounding factors. Especially investigations that rely on reuse of existing data, e.g., in the case of rare disease studies which often also suffer from low sample sizes, will benefit from this approach.

Data
In this study, we analyzed RNA-seq and metabolomics data from four large Dutch population cohorts: LifeLines (LL) [11], Leiden Longevity Study (LLS) [12], Netherlands Twin Register (NTR) [13], and Rotterdam Study (RS) [14,31]. The data is provided by the Dutch node of the European Biobanking and BioMolecular Resources and Research Infrastructure (BBMRI-NL).
RNA-seq data of all four cohorts was generated by the BBMRI-NL Biobank-based Integrative Omics Study (BIOS) Consortium at the Human Genotyping facility (HugeF) of ErasmusMC, the Netherlands. RNA sample processing and sequencing is described in detail by Zhernakova et al. [32]. Briefly, total RNA was extracted from whole blood, depleted of globin transcripts, and paired-end sequencing of 2x50-bp reads was conducted using the Illumina HiSeq 2000 platform. Read alignment to reference genome hg19 was performed using STAR (v2.3.0). We used the "Freeze2 unrelated data sets", which contain maximum sets of unrelated individuals and are available within the BIOS workspace at the SURF Research Cloud via the R package BBMRIomics v3.4.2 [33].
Metabolomics data was generated by the BBMRI-NL Metabolomics Consortium in 2014 as described by van den Akker et al. [26]. Briefly, metabolite concentrations were measured in EDTA plasma by proton nuclear magnetic resonance ( 1 H-NMR) spectroscopy on the platform of the Nightingale Health Group (Helsinki, Finland) [34].

Data analysis
All analyses were implemented in an R v4.0.3 [35] workflow employing R packages renv v0.14.0 for package management and drake v7.13.2 for workflow management. The analyses were run in the BIOS workspace of the SURF Research Cloud which is part of the multiomics analysis platform of BBMRI-NL. The code to run the analyses in available in GitHub [36] and archived in Zenodo [37].

Data preprocessing Normalization of values for clinical traits
Numeric values were used for all reported clinical variables. In case of categorical variables, they were binarized as follows. For smoking status, "current smoker" was coded as 1, and "former-smoker" and "non-smoker" were coded as 0; for sex, "male" and "female" were coded as 1 and 0, respectively; for lipid medication, "statins" were coded as 1, and "no" and "yes, but no statins" were coded as 0. In order to be able to compare effect sizes in association studies, all clinical variables were standardized to zero-mean and unit-variance (z-score normalization).

RNA-seq data preprocessing
Samples with more than 10% missing values in the RNAseq data were excluded from the analysis. Additionally, for comparisons of models employing either reported or inferred values as outcome, samples missing reported values were excluded. Subsequently, features (transcripts) missing in more than 10% of the samples were removed from the data set. Number of retained samples and features are given in Additional file 1.
The RNA-seq read counts as provided by the BBMRIomics R package were then normalized and transformed based on a previous evaluation of analysis strategies [16] as follows. Scaling factors for library sizes were calculated using the trimmed mean of M-values (TMM) method [38] implemented in the R/Bioconductor package edgeR v3.32.1 [39]. Using these scaling factors to adjust for sequencing depth, counts were transformed to log 2 counts-per-million (CPM) reads, their meanvariance relationships were estimated using voom [40] implemented in the R/Bioconductor package limma v3.46.0, and the associated observation-level weights were used in the subsequent linear modeling to adjust for heteroscedasticity.

Metabolomics data preprocessing
The Nightingale Health metabolomics features inquired are the 56 variables selected by van den Akker et al. [26]. Outliers were identified as the samples having more than 1 missing observation (301 removed), more than 1 data point under the detection limit (49 removed) and having a value more than 5 standard deviations away from the overall mean observed within BBMRI-NL (0 removed). Remaining with a total of 12926 samples (LLS = 2343, LL = 1475, RS = 5136, NTR = 3972). The remaining 4210 missing values (0.58% of the entire dataset) were imputed as zero and the metabolomic features were zscaled using the mean and standard deviations observed in BBMRI-NL. The samples were further filtered based on availability of corresponding RNA-seq data. Sample sizes per cohort and model are listed in Additional file 1.

Metabolomic surrogates
Seventeen logistic regression models trained on up to 22 cohorts, including LL and RS, were applied to the dataset as described in Bizzarri et al. [10]. The surrogates used in this study are the posterior probabilities which represent how likely an individual is at risk for each of the inquired common clinical variables.

Expression profile association studies
To determine association of profiles with clinical outcomes, gene-wise linear regression models were fitted using limma [41]. Known potential biological (sex, age) and technical confounders (flow cell number, white blood cell composition) were included in the models. Association studies were performed separately for the respective type of outcome variable, i.e., reported variable or metabolomic surrogate. Parameters for each linear model are summarized in Additional file 1. We adjusted p-values and effect sizes for statistical bias and inflation using the Bayesian method bacon [15], which estimates bias and inflation as parameters from the empirical null distribution of test statistics (t-statistic). Additionally, p-values were adjusted for multiple testing using the false discovery rate (FDR) [42,43]. Results for two different variable types for the same outcome, were compared by calculating Pearson correlation coefficients of regression coefficients from gene-wise fitted models.

Meta-analysis
Leave-one-cohort-out meta-analyses and replication studies in left out cohort were performed as described in [16].

GSEA
Gene-set enrichment analyses (GSEA) were performed using the R/Bioconductor package fgsea v1.16.0 [44] and gene sets retrieved from the Reactome Pathway Database [22]. Genes were ranked by − log 10 (p b ) * |β b | with p b = bacon-corrected p-value and beta b = baconcorrected effect size. The number of permutations for initial estimation of p-values was set to 1 × 10 4 ; the boundary for calculating p-values was set to 1 × 10 −50 .
Additional file 1: TWAS model parameters. This .csv file contains variable names, sample and feature numbers, and covariates per TWAS model and cohort.
Additional file 2: Significantly enriched pathways from GSEA. This .html file contains plots showing significantly enriched pathways for each outcome variable. GSEAs are based on meta-analyzed TWAS results.