Analysing an allelic series of rare missense variants of CACNA1I in a Swedish schizophrenia cohort.

CACNA1I is implicated in the susceptibility to schizophrenia by large-scale genetic association studies of single nucleotide polymorphisms. However, the channelopathy of CACNA1I in schizophrenia is unknown. CACNA1I encodes CaV3.3, a neuronal voltage-gated calcium channel that underlies a subtype of T-type current that is important for neuronal excitability in the thalamic reticular nucleus and other regions of the brain. Here, we present an extensive functional characterization of 57 naturally occurring rare and common missense variants of CACNA1I derived from a Swedish schizophrenia cohort of more than 10 000 individuals. Our analysis of this allelic series of coding CACNA1I variants revealed that reduced CaV3.3 channel current density was the dominant phenotype associated with rare CACNA1I coding alleles derived from control subjects, whereas rare CACNA1I alleles from schizophrenia patients encoded CaV3.3 channels with altered responses to voltages. CACNA1I variants associated with altered current density primarily impact the ionic channel pore and those associated with altered responses to voltage impact the voltage-sensing domain. CaV3.3 variants associated with altered voltage dependence of the CaV3.3 channel and those associated with peak current density deficits were significantly segregated across affected and unaffected groups (Fisher's exact test, P = 0.034). Our results, together with recent data from the SCHEMA (Schizophrenia Exome Sequencing Meta-Analysis) cohort, suggest that reduced CaV3.3 function may protect against schizophrenia risk in rare cases. We subsequently modelled the effect of the biophysical properties of CaV3.3 channel variants on thalamic reticular nucleus excitability and found that compared with common variants, ultrarare CaV3.3-coding variants derived from control subjects significantly decreased thalamic reticular nucleus excitability (P = 0.011). When all rare variants were analysed, there was a non-significant trend between variants that reduced thalamic reticular nucleus excitability and variants that either had no effect or increased thalamic reticular nucleus excitability across disease status. Taken together, the results of our functional analysis of an allelic series of >50 CACNA1I variants in a schizophrenia cohort reveal that loss of function of CaV3.3 is a molecular phenotype associated with reduced disease risk burden, and our approach may serve as a template strategy for channelopathies in polygenic disorders.

Introduction

Schizophrenia is among the most disabling and chronic mental disorders[1,2] and has a high heritability.[3,4] However, the molecular processes that govern the aetiology of the disease are not understood. Recent large-scale schizophrenia genome-wide association studies (GWAS) have identified many genomic loci associated with schizophrenia susceptibility[5] that may converge on common biological pathways underlying its pathophysiology.[6] Deciphering how disease-associated common variants confer risk for schizophrenia is challenging, in part because the causative risk alleles within each GWAS locus are difficult to pinpoint and they typically map to non-coding regions of unknown function. Similarly, genetic studies powered by whole-exome sequencing revealed that the combined contribution of rare coding variations in neuronal genes increases schizophrenia risk.[7-9] Recently, exome sequencing revealed a handful of genes as risk factors for schizophrenia.[10] Interestingly, rare coding variants are enriched in genes mapped within schizophreniaGWAS risk loci,[7,10,11] suggesting that both rare and common variants within certain genetic loci may account for the risk variation in schizophrenia penetrance across individuals.

Among the 108 genomic loci implicated in schizophrenia risk by GWAS,[5] one region resides within the CACNA1I gene (Supplementary Fig. 1A). Interestingly, CACNA1I is implicated in human cognitive function[12] and impulsive traits,[13] and the deficits of these two phenotypes have also been associated with schizophrenia patients.[14,15] As the combined contribution of many risk variants is thought to be necessary for schizophrenia pathophysiology, a systematic strategy to characterize the functional impact of multiple risk alleles within a risk gene is necessary to gain insights into the mechanisms of the disease. CACNA1I encodes the voltage-gated CaV3.3 calcium ion channel, which is enriched in neurons of the thalamic reticular nucleus (TRN) as well as a subset of cortical and hippocampal neurons.[16-19] While the specific function of CaV3.3 in hippocampal or cortical neurons is not well understood, its role in TRN neurons is well established. A hallmark feature of TRN neuronal excitability is rebound bursting, which is the firing of action potentials as the membrane potential recovers rapidly from a period of hyperpolarization.[20-22] In TRN neurons, the CaV3.3 channels mediate rebound bursting which is critical for generating sleep spindle oscillations during sleep.[23] Sleep spindles are significantly impaired in schizophrenia patients,[24] providing physiological relevance of CaV3.3 in the aetiology of the disease. Previously, we showed that a de novo coding variant of CACNA1I (R1346H), identified in one schizophrenia proband,[25] impaired CaV3.3 channel function.[26] Mice harbouring R1346H exhibit impaired TRN excitability and a reduction in the production of sleep spindles during non-rapid eye movement sleep.[27] In regard to monogenic disorders, animal models carrying disease-causing mutations can provide important clues about the role of the gene in the pathological mechanisms of the disease.[28,29] However, in polygenic disorders such as schizophrenia, where the disease risk is distributed across many genes, it is necessary to obtain a comprehensive analysis of the variants present not only in patients but also unaffected individuals and ideally to explore variants with a broad range of allele frequencies of the gene of interest to determine the full landscape of the disease risk burden.

In this study, we analysed an allelic series of 57 missense CACNA1I variants derived from a Swedish schizophrenia case and control cohort of more than 10 000 individuals. While the total number of missense variants in the schizophrenia and unaffected individuals’ groups were not different, we hypothesized that the variants derived from schizophrenia patients encode CaV3.3 channels that may display distinct biophysical properties from those encoded by common coding variants or from unaffected individuals in this cohort. We used an isogenic inducible expression system to assess the functional impact of individual CACNA1I coding variants as previously reported.[26] The essential biophysical properties of each variant-encoded CaV3.3 channel were recorded and analysed using an automated planar patch-clamp instrument (SyncroPatch 384PE) and compared to wild-type channels. We first set out to identify which biophysical properties of the CaV3.3 channel were impacted by common and rare variants and to what extent they were impacted. Subsequently, and more relevant to the disease pathophysiology, we asked whether there was functional segregation between the missense variants from schizophrenia patients and those from unaffected individuals. These biophysical properties of CaV3.3 channel variants were then used to model the influence on TRN excitability[30] as a measure of physiological impact. Such functional analyses of the allelic series of variant CaV3.3 channels, in the context of human genetics, may provide important insight into the role of CACNA1I in schizophrenia risk.

Materials and methods

Methods for molecular and cellular biology and biochemistry, electrophysiological solutions and automated patch-clamp experimental details can be found in the Supplementary material.

Electrophysiology

Electrophysiological protocols and biophysical analyses

The holding potential was set to −100 mV for all voltage protocols. Peak current and the activation steady-state parameters were obtained from CaV3.3 currents elicited by 1 s depolarization steps from −120 to +20 mV at 10 mV increments per sweep. The current–voltage relationship (I-V) was constructed by plotting the maximum peak current magnitude, normalized by the cell capacitance, as a function of the voltage applied. The voltage–conductance relationship [G(V)] was calculated by dividing the maximum current magnitude at each voltage by the corresponding driving force:The reversal potential (Vrev) was calculated from the intersection with the x-axis of the linear extrapolation of the last four points of the I-V curve. The voltage-dependent activation parameters were obtained by fitting a single Boltzmann function to the normalized conductance (G(V)/Gmax):V1/2 ACT is the voltage at which half of Gmax is reached; zδ corresponds to the slope of the Boltzmann function (voltage sensitivity); R corresponds to the universal gas constant; T is the absolute temperature in Kelvin; and F is the Faraday constant. The mean midpoint values were represented as described by Horrigan and Aldrich.[31]

Steady-state inactivation parameters were obtained by a 1 s prepulse from −120 mV to +20 mV at 10 mV increments per sweep, followed by a 300 ms test pulse at −20 mV (Supplementary Fig. 2). The current amplitudes (I) at the test pulse were normalized to the maximum current amplitude at −120 mV (Imax) and represents the fraction of channels able to open at a given prepulse. The voltage-dependent inactivation relationship was fitted by a Boltzmann function according to the following equation:where V1/2 INACT is the voltage at which half of Imax is reached; zδ corresponds to the slope of the fitted Boltzmann function (voltage sensitivity); and R, T, and F are the same as in equation (2).

The recovery from inactivation time constant (τ) was obtained by applying a paired-pulse protocol (Supplementary Fig. 2) with a preconditioning pulse followed by a test pulse at −20 mV for 200 ms. The recovery interpulse had variable length, ranging from 10 to 1600 ms. The ratio between the current magnitude of the test pulse and the preconditioning pulse was plotted as a function of the interpulse time length. The data were fitted by a monoexponential function as follows:where I1 and I2 represent the current magnitudes elicited by the preconditioning and test pulse, respectively. A0 and Aend indicate the initial and plateau values of the exponential function, t is the interpulse length, and τ corresponds to the time constant of the exponential fit.

Biophysical parameter normalization

To calculate the change in the midpoint of the activation and inactivation voltages (ΔV1/2ACT and ΔV1/2INACT) for the analysed CaV3.3 channels, the values for each cell expressing wild-type (WT) hCaV3.3 [equation (5)] or a variant channel [equation (6)] were normalized to the mean wild-type hCaV3.3 values for V1/2ACT and V1/2INACT and presented as the mean shift in mV from the mean value of wild-type hCaV3.3 channels.The peak current density and recovery time from inactivation of the variant channels are presented as the fold change with respect to the mean values of WT hCaV3.3 channels. We used equation (7) for WT and equation (8) for variant channels.To calculate the z-scores, the specific biophysical parameters (voltage dependence activation and inactivation, current density, recovery from inactivation) of the different variants were transformed into dimensionless distance from the WT properties measured by the standard deviation according to the following equation:where μ, σ and n stand for the mean, standard deviation and size of the sample, respectively.

Statistics

For the electrophysiology recordings, in total, we analysed 3312 cells that each expressed one missense variant CaV3.3 channel, along with 924 cells expressing the wild-type channels. Each variant hCaV3.3 channel that carried a missense mutation was recorded from at least two independent induction experiments, and at least 30 cells expressing the same variant channel were used in the analyses. The quality control parameters for the Syncropatch recording and analyses were previously described.[32] We ran eight batches of experiments to characterize 57 variants. We recorded five variant channels and the wild-type channel on each 384-well plate and recorded 10 variant channels in each batch. Every parameter obtained for each cell carrying a variant was normalized to the mean wild-type hCaV3.3 signal on the same plate. Then, we conducted ANOVA for the normalized values from each variant and the normalized values for the wild-type channel within the same batch, followed by Dunnett’s post hoc correction for multiple comparisons for each batch (10 variants and the wild-type). The significance threshold was defined as P ≤ 0.001 after the post hoc correction. Data are presented as the mean ± SEM, unless stated otherwise. For other analyses, ANOVA was performed for group comparisons followed by Dunnett’s post hoc comparison unless otherwise noted. R was used to calculate the P-values and odds ratios (ORs) in Fisher’s exact tests.

Homology model of hCaV3.3

After a BLASTP search against the PDB,[33] human CaV3.1 (6KZO) was selected as a template (62–68% identities in each transmembrane domain repeat).[34] Unresolved cytosolic regions in the template were predicted by the employment of the I-TASSER[35] and JPred servers.[36] All fragments were combined, and energy minimization was performed to reduce amino acid clashes by using UCSF Chimera.[37]

NEURON simulations

We utilized a published model (ModelDB, #17663) as previously described.[26] We simulated the rebound bursting of TRN neurons by considering the significantly altered biophysical properties of each variant CaV3.3 channel. As some variants modified more than one biophysical property, we considered all significantly altered properties in the modelling. For this purpose, we simulated TRN excitability in heterozygosity: 50% of the simulated channels presented significant changes in the biophysical properties recorded for that particular variant and the other 50% were simulated according to the wild-type parameters. The model was used to categorize the different variants as gain-of-function or loss-of-function according to the number of spikes during the rebound bursting compared to wild-type at a −0.25 nA current injection (Fig. 4). Notably, the C1498R variant produced a channel with an altered inactivation profile at more negative potentials (Supplementary Fig. 3). While the current NEURON model does not consider the closed-state inactivation of CaV3.3, this variant produced only ∼50% of the current density compared to wild-type, and is categorized as loss-of-function in TRN excitability. The C1335S variant is included in the modelling and is a clear loss-of-function.

Figure 4

Modelling the impact of biophysical properties of hCa (A) Simulated sample traces of rebound bursting from a model TRN neuron responding to hyperpolarizing current injections of −0.1, −0.25 and −0.4 nA using NEURON. Middle: A wild-type TRN neuron. Left: A TRN neuron expressing hCaV3.3 channels with ΔV1/2ACT = −5 mV. Right: A TRN neuron expressing hCaV3.3 channels with ΔV1/2ACT = +1 mV. (B) The relationship between the ΔV1/2ACT of hCaV3.3 channels and the resultant properties of rebound bursting, including the maximum number of spikes in the rebound firing (top), threshold of the amount of hyperpolarization to rebound (middle) and latency to rebound (bottom) induced by a −0.2 nA current injection. (C) Systematic analysis of the impact of four biophysical parameters of hCaV3.3 on the number of action potentials (spikes) of TRN rebound bursting. Heat maps of the number of spikes in the rebound burst as a function of the injected current (from −0.1 to −0.4 nA at a step size of 0.05 nA) were produced by perturbing each of the four properties in the wild-type hCaV3.3 channels. The observed differences in the hCaV3.3 properties are simulated and vary on the x-axis. The rebound bursting output in the presence of wild-type hCaV3.3 channels in each graph is indicated with a thin red line.

Data availability

Raw data were generated at Stanley Center for Psychiatric Research at Broad Institute. Derived data supporting the findings of this study are available from the corresponding author on request.

Results

CACNA1I URVs are preferentially located in conserved transmembrane domains that lack common variants

We analysed three types of missense variants of CACNA1I derived from a Swedish schizophrenia cohort[7] based on minor allele frequencies. Ultrarare variants (URV) occur less than one in 100 000 individuals in the general population according to gnomAD,[38] while rare variants occur in less than one in 1000 individuals (Supplementary Table 1A). Common variants refer to those alleles with frequencies greater than 0.1%. We identified 23 missense URVs of CACNA1I to be analysed in this cohort, of which 10 were found exclusively in schizophrenia patients and 13 in unaffected individuals. We identified an additional 27 rare variants, together with seven common coding variants of CACNA1I. The most common mutation of CACNA1I in this cohort was I1040V, which represented 63% of the total alleles (Table 1). In this report, we excluded missense variants located in the C-terminal cytosolic region of channel CaV3.3 (Supplementary Table 1B), as they are all located after the proximal functional motif of the C-terminus and deletion beyond this motif has little functional impact on CaV3.3 channel properties in the heterologous expression system.[39] In addition, none of the variants in the C-terminus were predicted to have a significant functional impact according to a machine learning model[40] trained on the sequence and structure features of the ion channels (Supplementary Fig. 1B). The carboxyl cytosolic region is implicated in protein–protein interactions that are unlikely to be reconstituted in the heterologous expression system. In contrast, amino acid variations within the transmembrane domains of the ion channels, which affect biophysical properties, are reliably assessed in a number of heterologous expression systems.[41,42]

Table 1

Functional impact of ultrarare and common variants of CACNA1I in this study

Classification[a]	Variant	CTL[b]	SCZ[c]	Type	Allele freq.[d]	ACT	INACT	PEAK	REC
Schizophrenia	C1498R[e]	0	1	URV	0.000004061	0.54	3.01	−0.54	0.28
Unaffected	M1605I	1	0	URV	0.000004061	−0.04	−4.23	−0.28	0.29
Unaffected	L214M	1	0	URV	0.000004061	3.58	0.55	−0.24	0.05
Unaffected	R196Q	1	0	URV	0.000004061	0.83	−2.92	−0.23	0.06
Schizophrenia	A832G	0	1	URV	0.000004061	−4.31	−0.84	0.00	0.02
Unaffected	C1335S	1	0	URV	0.000004062	–	–	−0.96	–
Unaffected	N1404K	1	0	URV	0.000004062	1.80	−0.44	−0.50	0.11
Schizophrenia	F1547L	0	1	URV	0.000004062	1.06	1.39	0.13	0.02
Schizophrenia	D43G	0	1	URV	0.000004065	−0.22	1.18	−0.02	0.06
Schizophrenia	N643D	0	1	URV	0.000004067	3.94	1.15	−0.11	0.03
Schizophrenia	E39K	0	1	URV	0.000004069	2.12	−0.64	−0.07	0.07
Unaffected	G495E	1	0	URV	0.000004104	0.69	1.92	0.20	0.06
Unaffected	S6Y	1	0	URV	0.000004369	0.20	−0.85	−0.04	0.11
Schizophrenia	V1705L	0	2	URV	0.000004828	−2.43	0.18	0.04	0.06
Unaffected	S990R	1	0	URV	0.000006209	1.32	2.59	−0.06	0.12
Schizophrenia	A1352T	0	1	URV	0.000007219	3.20	2.20	−0.41	0.07
Schizophrenia	V180M	0	1	URV	0.000007316	−5.57	−0.47	0.10	0.11
Unaffected	V1418I	1	0	URV	0.000008122	2.94	1.50	−0.58	0.13
Unaffected	S1482Y	2	0	URV	0.000008122	−2.50	−2.60	−0.27	0.13
Mix	S839C	1	1	URV	0.000008122	0.20	−2.24	0.05	−0.03
Unaffected	A232V	1	0	URV	0.000008123	1.00	−1.08	−0.50	−0.20
Unaffected	L936I	1	0	URV	0.000008196	−3.04	−0.98	0.00	−0.04
Schizophrenia	A929T	0	1	URV	0.000008212	1.59	1.78	0.07	0.14
Mix	D302G	15	10	Common	0.001385000	−0.57	−2.61	−0.13	−0.13
Mix	Q1193H	38	30	Common	0.003676000	−0.53	−2.33	−0.17	0.03
Mix	P1026L	90	50	Common	0.004077000	0.57	2.26	0.10	−0.01
Mix	H505N	41	34	Common	0.004177000	−0.39	1.11	−0.02	−0.14
Mix	R307H	160	127	Common	0.004883000	1.07	−1.64	0.02	−0.22
Unaffected	P7L	2	0	Common	0.010360000	−2.22	−0.52	0.01	0.04
Mix	I1040V	2634	2053	Common	0.631500000	1.34	1.66	0.17	0.01

Values in bold are statistically significant.

ACT = ΔV1/2ACT, the shift in the voltage dependence of activation respect to wild-type (mV). INACT= ΔV1/2INACT, the shift in the voltage dependence of inactivation respect to wild-type (mV). PEAK corresponds to the normalized peak current change and REC refers to normalized recovery from inactivation time based on equations 7 and 8.

Clinical classification of the variants in the Swedish cohort study. CTL = unaffected individuals; SCZ = schizophrenia patients; Mix = found in both patients and unaffected individuals.

Number of alleles identified in controls (CTL) and patients (SCZ) carrying the mutation in the Swedish cohort.

Allele frequency calculated from the number of alleles in gnomAD (https://gnomad.broadinstitute.org/) carrying the mutation.

C1498R was the only variant that presented an abnormal inactivation profile (see the ‘Results’ section).

CACNA1I encodes the core α1 subunit of the hCaV3.3 channel, which are composed of four pseudosymmetric domains (DI–DIV), each composed of six transmembrane segments (S1–S6), as shown in Fig. 1A. Segments S1–S4 constitute the voltage-sensing domains (VSDs), including the S4 helix, which is essential for sensing changes in the transmembrane voltage, while the S5–S6 segments constitute the pore domain, responsible for forming the central ion conduction pore.[43,44] We found that all seven common variants of CACNA1I code for changes in the intracellular or extracellular segments of the hCaV3.3 protein, while many URVs of CACNA1I lead to changes in conserved transmembrane regions throughout the hCaV3.3 protein (Fig. 1A).

Figure 1

Ultrarare and common coding variants of (A) Schematic representation of the membrane topology of the hCaV3.3 channel encoded by CACNA1I. Four transmembrane domains are shown as DI–DIV. Each domain is composed of six transmembrane segments (S1–S6). The location of the ultrarare variant (URV) mutations is indicated in magenta (schizophrenia, SCZ) and green (unaffected), and that of the common variants is shown in blue circles. (B) The cartoons represent three major structural conformations of the hCaV3.3 channel—closed ‘C’: the VSD is ready to translocate, and the pore domain cannot conduct ions; open ‘O’: the VSD is translocated to the upward position and the pore domain can conduct ions; inactivated ‘I’: the VSD assumes a different upward configuration and the pore domain cannot conduct ions. The transitions between the three conformations depend on the membrane voltage. (C) Biophysical characterization of the URVs and common variants of CACNA1I derived from the Swedish cohort. Variants were grouped into common variants (blue), URVs derived from schizophrenia (magenta) and URVs derived from unaffected individuals (green), as displayed. x = 0 in each bar graph corresponds to the wild-type (WT) hCaV3.3 channel behaviour. Voltage-dependent channel activation, ‘Activation’, shows the change in the V1/2ACT (ΔV1/2ACT) for each variant hCaV3.3 channel with respect to the wild-type channel (x = 0). Each bar represents one hCaV3.3 variant, and the graph is organized from the largest rightward shift to the largest leftward shift (top to bottom) for the three groups. The presentation is similar for ‘Inactivation’, where ΔV1/2INACT for each variant hCaV3.3 channel is shown. ‘Peak current’ shows the peak current densities for each variant hCaV3.3 channel, normalized to those of the wild-type channel (x = 0). Reductions in the peak current densities are shown as negative values. For ‘Recovery time’, the time constant (τ) of recovery for each hCaV3.3 variant channel is normalized to that of the wild-type channel (x = 0); negative values represent smaller τ values and faster recovery times. Asterisk denotes significant differences from the wild-type channel; see Table 1. Individual curves for each of the variants are shown in Supplementary Fig. 3.

CaV channels exist in conformations that roughly correspond to closed states, open states and inactive states (Fig. 1B).[45-47] The relative occupancy in open, closed and inactive hCaV3.3 channel states defines how the channel responds to physiological stimuli in neurons.[48-50] We established electrophysiological recording protocols (‘Materials and methods’ section) to probe the transition between the open (O), closed (C) and inactivated (I) states of hCaV3.3 channels by characterizing four key biophysical properties: voltage dependence of activation, voltage dependence of inactivation, peak current density and time course of recovery from inactivation (Supplementary Fig. 2). These biophysical descriptions produce a comprehensive functional landscape for each hCaV3.3 channel that harbours a missense change relative to the reference hCaV3.3 channel.

URVs mainly alter the voltage-dependent activation and the current density of hCaV3.3 channels

Voltage-dependent activation

Six of 10 URVs from the schizophrenia group (E39K, V180M, N643D, A832G, A1352T and V1705L) and 4 of 13 URVs from the unaffected individuals (L214M, L936I, V1418I and S1482Y) produced hCaV3.3 channels that opened at membrane voltages that were significantly different from those of wild-type hCaV3.3 channels (Fig. 1C, ‘Activation’, magenta and green bars, ACT in Table 1). The voltage dependence of activation of hCaV3.3 channels coded from any of the seven common variants was indistinguishable from that of wild-type hCaV3.3 channels (Fig. 1C, blue bars and Table 1). None of the URVs or common variations induced changes in the steepness of the hCaV3.3 channel activation curves (Supplementary Table 2). Notably, hCaV3.3 channels coded from the schizophrenia variant V180M opened at the most negative membrane voltages relative to those of the wild-type hCaV3.3 channels (ΔV½ = −5.6 ± 0.4 mV, P < 0.0001), while hCaV3.3-N643D channels opened at the most positive voltage (ΔV½ = +3.9 ± 0.4 mV, P < 0.0001). In total, 10 of 23 (43%) URVs significantly altered the hCaV3.3 channel response to changes in membrane depolarization, as defined by the V1/2ACT values relative to the wild-type control, while none of the common variants did. These results indicate that ultrarare hCaV3.3 variants, but not common hCaV3.3 variants, alter voltage-dependent activation of the hCaV3.3 channel with respect to the wild-type control channel.

Voltage-dependent inactivation

Three of 13 URVs (R196Q, S990R and M1605I) from unaffected individuals produced hCaV3.3 channels with altered voltage dependence of inactivation (V1/2INACT) relative to the wild-type channels, and 2 of 10 URVs from the schizophrenia group significantly influenced the voltage dependence of inactivation (Fig. 1C, ‘Inactivation’, INACT in Table 1). Notably, hCaV3.3 channels encoded by URVs carrying M1605I were inactivated most negatively (ΔV1/2INACT = −4.2 ± 0.7 mV, P < 0.0001) relative to wild-type hCaV3.3 channels, while channels encoded by URVs carrying C1498R were inactivated most positively (ΔV1/2INACT = +3.0 ± 0.5 mV, P < 0.0001). In addition to altering V1/2INACT, C1498R produced channels with altered closed-state inactivation profiles at more negative potentials (Supplementary Fig. 3). None of the hCaV3.3 channels coded by a common variant inactivated at voltages differed significantly from those of wild-type channels. Altogether, 5 of 23 (22%) URVs and none of the common variants produced channels that inactivated at voltages that were different from wild-type hCaV3.3.

Peak current density

We found that nearly half of the hCaV3.3 channels produced by URVs (10 of 23 total URVs) had peak current densities that were significantly lower than those of wild-type channels, whereas the whole-cell current densities of hCaV3.3 channels encoded by each of the seven common variants were similar to those of the wild-type channels (Fig. 1C, ‘Peak current’, PEAK in Table 1). The majority of the URVs that influenced hCaV3.3 current density were identified in unaffected individuals (8 of 13; R196Q, L214M, A232V, C1335S, N1404K, V1418I, S1482Y and M1605I; Fig. 1C). Notably, we were not able to detect any currents in hCaV3.3 channels from variants that harboured C1335S (Δfold change = −0.96 ± 0.06, P < 0.0001). Among the 10 URVs from the schizophrenia group, only those harbouring C1498R and A1352T reduced the peak current density of the channel (Δfold change = −0.54 ± 0.04 and −0.41 ± 0.06, respectively, P < 0.0001). Interestingly, A1352T is located next to C1335S, close to a glycosylation residue previously characterized by our group as being crucial for protein stability.[26] None of the common variants altered the current density significantly, and none of the URVs analysed significantly increased the current densities relative to wild-type hCaV3.3 channels. In summary, nearly half of hCaV3.3 channels carrying CACNA1I URVs are associated with lower current densities than wild-type channels, while common variants do not appear to influence hCaV3.3 current density.

Time course of recovery from inactivation

We found that 5 of 23 (22%) hCaV3.3 URV channels (A232V, V1418I, S1482Y, C1498R and M1605I) had an altered time course of recovery from inactivation relative to wild-type channels, while no common hCaV3.3 variant possessed such an altered property. One URV from the schizophrenia group (C1498R) and three URVs from the unaffected group (V1418I, S1482Y and M1605I) slowed the rate of recovery from inactivation, while one variant (A232V) sped up the recovery from inactivation by ∼20% with respect to wild-type (Fig. 1C, ‘Recovery Time’, REC in Table 1). All five URVs that significantly changed the rate of recovery from inactivation of the channel led to either inactivation at voltages different from those of wild-type or peak current deficits, suggesting that the structural regions affected by these variants influence multiple aspects of the channel properties.

In total, 8 of 13 (62%) URVs of CACNA1I from the unaffected individuals altered the peak current density, whereas 6 of 10 (60%) URVs derived from schizophrenia patients changed the voltage dependence of hCaV3.3 channel activation. These results strongly suggest that voltage dependence of activation and peak current density are the most relevant properties in the context of the human genetics of CACNA1I coding variants.

URVs of CACNA1I derived from patients and unaffected individuals alter distinct channel properties

To visualize and compare the functional impact of coding variants across multiple biophysical properties, we calculated the functional differences among hCaV3.3 variants relative to wild-type channels with dimensionless z-scores. These z-scores were calculated for each of the four biophysical properties of each variant channel (see the ‘Materials and methods’ section and Supplementary Table 3). Using the z-scores, the overall functional impact of each missense variant was visualized in a rhombus plot, where each corner represents a distinct biophysical parameter (Fig. 2, inset). The z-score rhombuses were grouped for the URVs derived from schizophrenia patients (magenta) and unaffected individuals (green) and for common variants (blue) (Fig. 2A). The C1335S mutation of the control group did not generate any detectable calcium currents (z-score = −26 in peak currents; Supplementary Table 3). As it was not possible to calculate the remaining biophysical properties, this mutation was excluded from the z-score analyses. The resultant vector of all z-scores (thick arrow) represents the collective magnitude and direction of all variants in each group (Fig. 2A and Supplementary Fig. 4). The magnitudes of the resultant vectors for the URVs from schizophrenia (z-score = 25.4) and those from unaffected individuals (26.1) were similar, but the vector directions were distinct (north versus west), reflecting different overall functional phenotypes. In comparison, the resultant vector from common variants (blue arrow, Fig. 2A) exhibited a much smaller magnitude of 8.9 and pointed eastward (voltage dependence of inactivation). A statistical comparison across the three groups in each functional property is shown in Fig. 2B. CaV3.3 channels generated from URVs derived from the schizophrenia group open at voltages that are significantly different from those of channels coded by common variants (schizophrenia: average z-score = 4.6 ± 1.0; common: 1.6 ± 0.4; P = 0.031, ANOVA). In contrast, URVs derived from the unaffected group produced hCaV3.3 channels with significantly smaller peak current densities than those from common coding variants (unaffected: average z-score = 4.9 ± 1.3; common: 1.4 ± 0.4; P = 0.036, ANOVA). In summary, our analyses show that hCaV3.3 URVs from the patient group mainly altered the activation of the channel in response to voltages, whereas URVs from the unaffected group preferentially reduced peak current densities.

Figure 2

Functional comparison among URVs from schizophrenia patients, URVs from unaffected individuals and common variants using (A) The four biophysical properties for each variant hCaV3.3 channel with detectable currents are visualized using dimensionless z-scores in a rhombus. The C1335S variant did not lead to detectable currents and was excluded from this analysis. The data shown here correspond to 29 of the 30 variants presented in Fig. 1. Grouped rhombus plots from common variants (blue), schizophrenia (SCZ) URVs (magenta) and unaffected URVs (green) are shown. Each rhombus represents a variant, where each corner reflects the z-score of one biophysical parameter (inset). The solid black line in each rhombus represents the wild-type behaviour, where z-score = 0, and dashed lines indicate ± 2.5 z-score away from the wild-type. The average resultant vector for each group is shown as a thick arrow. (B) Group-level comparison of absolute z-scores for each of the four properties of the hCaV3.3 channels. For V1/2ACT, absolute z-score = 4.6 ± 1.1 (SCZ), 2.7 ± 0.7 (unaffected) and 1.6 ± 0.4 (common) with Welch’s ANOVA (F = 4.019, P = 0.038). For the peak current density, z-score = 3.0 ± 1.4 (SCZ), 4.9 ± 1.3 (unaffected) and 1.4 ± 0.4 (common) with Welch’s ANOVA (F = 3.697, P = 0.049). Asterisk denotes P < 0.05 for post hoc statistical significance in pairwise comparison.

URVs affecting activation and current densities are located in distinct structural regions

We were able to map 8 (of 10) URVs (V180M, L214M, N643D, A832G, A1352T, V1418I, S1482Y and V1705L) that impaired voltage dependence of activation in a structural homology model of hCaV3.3 based on the published hCaV3.1 structure[51] (Fig. 3A). These eight URVs encode for proteins in three different regions (Fig. 3B). ‘Region I’ is the voltage-sensing domain and is where V180M, N643D and S1482Y reside. V180M is located within the S4 segment of domain I (DI) and is flanked by two positively charged arginine residues, R178 and R181 (Fig. 3C, Region I). The V-to-M mutation alters the hydrophobicity of this residue and may modify how the S4 segment translocates through the hydrophobic membrane core in response to membrane depolarization,[43] thus inducing a change in V1/2ACT (ΔV1/2ACT = −5.6 ± 0.4 mV). In comparison, N643D resides at the bottom of the S1 segment of domain II (DII), ∼9 Å away from R743 in the S4 helix of DII, and induces a rightward shift in V1/2ACT (ΔV1/2ACT = 3.9 ± 0.4 mV). These results are consistent with the well-established role of the positively charged residues of S4 in voltage-sensing mechanisms.[47,52] S1482Y is located at the beginning of the S1 segment in DIV, and it induced a leftward shift in V1/2ACT (ΔV1/2ACT = −2.5 ± 0.4 mV). The electrostatic equilibrium of the VSD can be modified not only by residues located or interacting with S4 segment, but also those located in adjacent segments (S1–S3) within the VSD,[53,54] which may underlie how this variant changes the voltage dependence.

Figure 3

3D modelling revealed structural hotspots for voltage dependence of activation and channel expression. (A) Side view (left), top view (middle) and bottom view (right) of the homology structural model of hCaV3.3 generated from the cryo-EM structure of hCaV3.1 (PDB ID, 6KZP). (B) The four transmembrane domains (DI–DIV) are structurally aligned to DI, and the residues corresponding to the variants that induced significant changes in V1/2ACT are indicated in magenta (schizophrenia, SCZ) and green (unaffected). (C) Three regions (I, II and III) harbour the mutations from variants that changed the voltage dependence of activation. The positively charged arginine residues in S4 are shown in pale cyan (R1-K5) in Region I. The selectivity filter residues for Ca2+ selectivity are highlighted in pale cyan for Regions II and III. (D) Ten variants that changed the peak current density are mapped onto the transmembrane domain alignment of hCaV3.3. The selectivity filter residues are highlighted in pale cyan. (E) Top: Analyses of the expression levels of hCaV3.3 channel proteins in HEK293 cells using western blotting with an antibody against the FLAG tag. Vinculin is the loading control. Dotted line indicated where the blots were cut with the raw images presented in the Supplementary material. Bottom left: Densitometry of the total amount of hCaV3.3 channel protein in cell lysate for all variants (n = 3–6; *P < 0.05), normalized to the wild-type (WT) hCaV3.3. Bottom right: The relationship between the expression levels of hCaV3.3 and the peak current densities for URVs in Fig. 1. Fitted linear regression (magenta line) is shown with 95% confidence interval (dotted line), demonstrating a slope of 1.02, r = 0.41, and P = 0.033.

The second region (Fig. 3C, Region II) is located in the S5–S6 linkers near the membrane re-entering loop (P-loop helix) that contains the selectivity filter motif.[55] A832G, A1352T and V1705L are located on the top of the S6 segments of DII, DIII and DIV, respectively. Interestingly, A832G and V1705L both shifted the response of hCaV3.3 channels to more negative potentials (ΔV1/2ACT = −4.3 ± 0.7 mV and −2.4 ± 0.5 mV, respectively), while A1352T shifted the activation curve to more depolarized potentials (ΔV1/2ACT = 3.2 ± 0.6 mV). Unexpectedly, this region modulates hCaV3.3 voltage activation (V1/2ACT); however, analogous regions have been reported to influence the coupling between the voltage sensing domains and the pore opening in KV and TRP channels.[56,57] ‘Region III’ resides within the S5–S6 transmembrane segments that form the ion pore (Fig. 3C, Region III). Specifically, L214M is located in the middle of the S5 segment in DI and V1418I in the S6 segment of DIII; both variants right shifted the voltage dependence of the activation curve (ΔV1/2ACT = +3.6 ± 0.4 mV for L214M, ΔV1/2ACT = +2.8 ± 0.7 mV for V1418I). In our homology model, V1418I and L214M are within 4 Å of the P-loop helix of the selectivity filter that gates the calcium ion flow. Both residues altered V1/2ACT and reduced Ca2+ currents by 30% and 60%, respectively, suggesting that this region is important for both voltage dependence and ionic current flow.

We next examined the structural locations of the mutations of URVs that altered hCaV3.3 peak current densities and mapped all 10 coding URVs associated with smaller hCaV3.3 peak current densities (eight from unaffected and two from schizophrenia individuals) in our homology model (Fig. 3D). Eight of these 10 URV mutations were located in the S5–S6 pore regions in the structural model, while the remaining two resided in DIII-S1 (C1498R) and adjacent to DIV-S1 (S1482Y; Fig. 3D). To understand the biochemical mechanism underlying the reduced hCaV3.3 channel currents, we analysed the total hCaV3.3 channel protein levels in cell lines that expressed the variant channels by Western blots. Compared with wild-type hCaV3.3, the channels produced by 9 of the 10 coding URVs had reduced protein levels in whole cell lysates (R196Q, L214M, A232V, C1335S, A1352T, N1404K, V1418I, C1498R, M1605I, P < 0.05, Fig. 3E). The linear correlation between total protein level and peak current deficit suggests that these URVs reduced hCaV3.3 currents mainly by disrupting overall channel protein stability or protein expression (Fig. 3E, scatter plot, slope 1.018, R = 0.41, P = 0.034). These results support that S5–S6 pore domains are the major pore structural motifs where mutations are most likely to disrupt hCaV3.3 stability and/or expression. We cannot rule out the contributions of additional mechanisms to the deficits of these variants, such as ionic conduction, for example, based on our current analyses. For channels harbouring S1482Y, deficits in single channel conductance, trafficking and permeability may explain the reduced current density, but additional analyses are needed to dissect the mechanisms.

In summary, the URVs that impaired the responses of hCaV3.3 channels to voltage changes reside in three distinct structural elements, and the majority of URV mutations that reduced hCaV3.3 current density are located within the S5–S6 ion pore domains where the missense changes seem to disrupt total protein levels.

NEURON simulation reveals the impact of the biophysical properties of the hCaV3.3 channel on TRN excitability

CaV3.3 channels are necessary for hyperpolarization-induced rebound bursting in the TRN, which is critical for sleep spindle generation.[23,27] To understand the physiological impact of the four functional properties of hCaV3.3 channels, we simulated rebound bursting using an 80-compartment TRN neuron model in NEURON (Supplementary Fig. 5).[26,30,58] Using this model, we discovered that changes in the voltage dependence of activation of hCaV3.3 channels (V1/2ACT) drastically alter the number of action potentials produced in rebound firing, the latency to rebound firing and the hyperpolarization threshold necessary to generate rebound firing (Fig. 4A and B). Specifically, the shift in the V1/2ACT of hCaV3.3 channels to more negative voltages increased the number of action potentials and lowered the threshold of hyperpolarization required and the latency to fire an action potential relative to wild-type hCaV3.3 channels (Fig. 4B). In contrast, a small shift in the V1/2ACT of hCaV3.3 channels to more depolarized voltages (1 mV) abolished rebound bursting. The overall impact of V1/2ACT on the rebound output with different current injection is summarized in a heat map (Fig. 4C, ‘Voltage-dependent activation’). We next investigated the impact of the voltage dependence of inactivation (V1/2INACT) of CaV3.3 on TRN rebound bursting. When V1/2INACT was right shifted +5 mV, the number of action potentials in the rebound increased 1.6-fold (from 12 to 20) in response to a modest hyperpolarization (−0.25 nA current injection), whereas a leftward shift of −1 mV in V1/2INACT abolished rebound bursting (Fig. 4C, ‘Voltage-dependent inactivation’). Therefore, in contrast to the effects of ΔV1/2ACT, the hyperpolarized shift in the V1/2INACT of hCaV3.3 channels reduced the TRN firing output, while depolarized shifts increased the TRN output.

In addition, we replicated our previous findings in which a reduced hCaV3.3 current density led to reduced rebound output[26] (Fig. 4C, ‘Peak current density’). Furthermore, a faster hCaV3.3 channel recovery time (smaller τ) from inactivation increased the number of action potentials in the rebound bursting, whereas a slower recovery time reduced the rebound burst output (Fig. 4C, ‘Recovery from inactivation’). The observed changes in recovery from inactivation kinetics only mildly altered TRN rebound burst output, suggesting that this property is not as influential as the other three parameters in modifying rebound bursting. Last, observed differences in hCaV3.3 current density[26] or V1/2ACT among the variants had little or no effect on TRN tonic firing (Supplementary Fig. 6). In summary, we used NEURON simulations to model the effects of each of the four biophysical properties on TRN neuron rebound bursting and predicted the impact of hCaV3.3 variants on TRN output.

Segregation of variant channel properties across the disease status

To examine variants beyond the URVs, we expanded our analyses to all 50 rare coding variants, including 27 additional 2rare variants in this cohort (Table 2), and confirmed that rare hCaV3.3 variant channels differed from wild-type channels mostly in two properties: V1/2ACT and peak current density (15 of 27 variants, Table 2), similar to the findings for the URVs (Fig. 1 and Table 1). Furthermore, we found that rare alleles with altered V1/2ACT and peak current density were significantly segregated across disease status with an OR of 0.22 (Fisher’s exact test, P = 0.034, Table 3 and Fig. 5A), suggesting that individuals with alleles carrying a mutation that leads to peak current deficits have an ∼80% reduced disease burden relative to individuals with alleles carrying mutations that modify the voltage sensitivity of the CaV3.3 channel. As most of the peak current deficits induced by coding variants are associated with total channel protein levels (Fig. 3), we then analysed the European subset of alleles carrying protein-truncating variants (PTVs) of CACNA1I in a large cohort of schizophrenia patients recently published by the SCHEMA consortium.[10] Consistent with our findings, analysis of the SCHEMA findings identified 17 CACNA1I PTVs in 196 644 control individuals and 0 PTVs in 48 496 schizophrenia patients of European origin (Fisher’s exact test, P = 0.034, OR = 0), suggesting that there are fewer CACNA1I PTVs among schizophrenia patients (Table 3 and Fig. 5A). Taken together, the results of our functional analyses of CACNA1I missense rare variants in the Swedish cohort and the analyses of the PTVs of CACNA1I in the SCHEMA cohort suggest that rare cases of loss of or reduction in function of CaV3.3 channels are significantly enriched in unaffected individuals with respect to the patients.

Table 2

Twenty-seven additional rare variants of CACNA1I identified in this study

Classification[a]	Variant	CTL[b]	SCZ[c]	Type	Allele freqency[d]	ACT	INACT	PEAK	REC
Unaffected	D805N	1	0	Rare	0.0000089	2.7	−2.7	−0.48	−0.06
Unaffected	R1689W	1	0	Rare	0.0000094	3.4	−1.4	−0.59	0.05
Schizophrenia	R327H	0	1	Rare	0.0000122	1.9	1.3	−0.66	0.34
Mix	R1272G	1	1	Rare	0.0000161	−3.1	4.4	−0.01	0.08
Unaffected	A615V	1	0	Rare	0.0000171	1.2	1.9	−0.12	0.01
Unaffected	R111G	3	0	Rare	0.0000179	2.2	1.2	−0.29	0.05
Mix	E436G	1	2	Rare	0.0000185	−7.5	−3.2	0.44	−0.14
Mix	R410Q	1	1	Rare	0.0000282	−0.5	2.6	0.01	0.10
Schizophrenia	R749Q	0	1	Rare	0.0000323	2.6	1.4	−0.15	0.07
Unaffected	A551T	1	0	Rare	0.0000323	−1.7	−1.3	−0.12	0.06
Unaffected	G635S	1	0	Rare	0.0000358	−2.1	−2.3	−0.07	0.15
Unaffected	R873H	1	0	Rare	0.0000399	1.3	1.8	0.05	−0.14
Mix	S1225A	1	2	Rare	0.0000401	1.7	0.9	−0.03	0.09
Schizophrenia	A589V	0	1	Rare	0.0000405	0.7	1.5	−0.11	0.07
Schizophrenia	R1545H	0	1	Rare	0.0000577	−1.3	0.8	−0.11	0.03
Mix	V825I	1	1	Rare	0.0000577	3.5	1.8	−0.11	0.02
Mix	R965H	5	2	Rare	0.0000634	4.7	4.1	−0.04	0.01
Unaffected	R627Q	1	0	Rare	0.0000668	−0.5	−0.6	−0.07	0.02
Mix	I1166V	12	13	Rare	0.0001260	−1.1	−0.9	−0.19	0.08
Mix	T534M	1	2	Rare	0.0001264	−0.3	1.4	0.03	0.01
Mix	R29W	2	1	Rare	0.0001342	1.0	3.0	0.02	0.04
Unaffected	P35L	2	0	Rare	0.0001667	−0.4	0.0	−0.11	−0.08
Schizophrenia	R452H	0	2	Rare	0.0001914	−3.6	−3.0	−0.12	0.08
Mix	G1231S	2	2	Rare	0.0002319	−1.4	−2.0	−0.05	0.09
Mix	M128L	3	3	Rare	0.0003158	−3.0	−0.8	0.00	0.05
Mix	I438T	21	18	Rare	0.0003511	−2.4	−1.9	0.12	0.08
Mix	V1535L	6	5	Rare	0.0004149	0.6	−0.8	−0.11	−0.12

Values in bold are statistically significant.

ACT corresponds to ΔV1/2ACT, the shift in the voltage dependence of activation and INACT to ΔV1/2INACT, the shift in the voltage dependence of inactivation respect to wild-type (mV). PEAK indicates to the normalized peak current change and REC to normalized recovery from inactivation time based on equations 7 and 8.

Clinical classification of the variants in the Swedish cohort study. CTL = unaffected individuals; SCZ = schizophrenia patients; Mix = found in both patients and unaffected individuals.

Number of alleles identified in controls (CTL) carrying the mutation in the Swedish cohort.

Number of alleles identified in schizophrenia patients (SCZ) carrying the mutation in the Swedish cohort.

Allele frequency calculated from the number of alleles carrying the mutation found in gnomAD (https://gnomad.broadinstitute.org/).

Table 3

Categorical analyses of rare coding alleles of CACNA1I in the Swedish schizophrenia cohort functionalized in this study, and alleles carrying PTVs of CACNA1I in the European population according to SCHEMA

	SCZ	Unaffected	Fisher’s test	Odds ratios	95% CI
Functional characterization of Swedish cohort
Alleles with peak current deficits	3	14	P = 0.036	0.22	0.03–0.98
Alleles with altered V1/2ACT	18	18
European SCHEMA data[a]
PTV alleles	0	17	P = 0.034	0	0–0.98
Total alleles	48 496	196 644

Values in bold are statistically significant. SCZ = schizophrenia.

SCHEMA consortium aggregates data from different genetic background, the numbers displayed here correspond only to the data derived from European and Northern European population according to SCHEMA classification (https://schema.broadinstitute.org).

Figure 5

Impact of heterozygous expression of (A) Forest plots of Table 3. The horizontal bars show the 95% confidence intervals of the odds ratio (OR, dot), and nominal significance using Fisher’s exact test is shown. The magenta bar indicates significant associations. (B) The relationship between the number of action potentials in TRN rebound bursting and the current injection for URVs and common variants using NEURON simulation. Each line represents the TRN response in the presence of a URV or a common variant of CACNA1I in heterozygosity. Common variants, blue; URVs from schizophrenia (SCZ), magenta; URVs from unaffected, green. Rebound responses may overlap. (C) Comparison of the rebound output of TRN shown in B. Each dot represents the TRN output in the presence of a variant channel with −0.25 nA current injection. Independent t-tests showed a significant reduction in the number of spikes for the common vs. control channels (P = 0.011). (D) The relationship between the number of action potentials in TRN rebound bursting and the injected current for all 57 variants of CACNA1I in heterozygosity. Gain-of-function and loss-of-function: increased and decreased number of action potentials in the rebound firing, respectively, with respect to the wild-type firing pattern (magenta line).(E) Forest plots of categorical analyses between disease status (schizophrenia versus unaffected) and number of functional alleles associated with predicted TRN phenotypes. The horizontal bars show the 95% confidence intervals of the odds ratio, and nominal significance using Fisher’s exact test is shown for each comparison. The relevant contingency tables are shown in Supplementary Table 4.

Using NEURON modelling, we next estimated the TRN rebound firing output in the presence of variant hCaV3.3 channels in heterozygosity (see Method). First, we compared the predicted TRN output (number of action potentials in rebound firing) in the presence of heterozygous URVs and common variants in response to various hyperpolarizing inputs (Fig. 5B). At a moderate hyperpolarizing input of −0.25 nA, TRN neurons with one allele expressing a common variant for hCaV3.3 channels produced an average output of 12 ± 0 spikes in the rebound versus 8.2 ± 2.3 spikes for those expressing URV channels derived from schizophrenia patients and 4.4 ± 1.7 spikes for those expressing URV channels derived from unaffected individuals (Fig. 5C). There was a significant difference between the output of the common variant group and that of the unaffected individual group (Brown–Forsythe ANOVA, F = 4.05 (2.0, 19.8), P = 0.034) and a significant difference in the standard deviation across the three groups (Bartlett test, P = 0.0037). These results suggest that URVs derived from unaffected subjects reduced TRN output with respect to common variants, consistent with the results in Fig. 5A.

Next, we analysed the physiological impact of all rare variants and categorized the rare coding variants of CACNA1I based on their impact on TRN output. The variants that did not change the rebound spike count with respect to wild-type were categorized as neutral, and those that changed the rebound bursting with respect to wild-type were categorized as active. Active variants can be further categorized as gain-of-function or loss-of-function depending on whether they increased or decreased TRN output, respectively (Fig. 5D). There is no genotype association in the number of alleles carrying neutral variants versus active variants (P = 0.50) or gain-of-function versus neutral variants (P = 0.83). There was a nonsignificant trend between loss-of-function and neutral alleles (OR = 0.5, P = 0.14) across the disease status and a similar trend between loss-of-function and gain-of-function alleles (OR = 0.41, P = 0.17; Fig. 5E). This suggests that the loss-of-function of TRN output might be a candidate phenotype that segregates between schizophrenia patients and unaffected individuals, but the result is not significant. Additional cohort information, such as from the SCHEMA, would increase the power for such analyses. Most surprisingly, these analyses indicated a phenotype in control populations that may reduce schizophrenia risk burden, highlighting the urgent need for functional characterization of coding missense variants in both control and affected populations. Our analyses provide a functional landscape in variant-to-function of CACNA1I to understand the role of risk genes in the pathophysiology of the disease.

Discussion

Monogenic channelopathies, including Dravet, Brugada and Timothy syndromes, defined by pathogenic coding variants with high penetrance, are relatively well characterized.[59-67] Schizophrenia, by contrast, is a common polygenic disorder that is thought to arise from the combined actions of multiple genetic factors, of which none, individually, are highly penetrant or causal.[68] This confounds our ability to address if and how a given allele contributes to disease pathophysiology. However, by analysing the function of a relatively large number of common, rare and ultrarare variants of CACNA1I and using a highly sensitive assay of channel protein function, we identified clear molecular phenotypes among rare variants of CACNA1I between the schizophrenia group and control group. Typically, the impact of a given variant is measured relative to the properties of a specific wild-type protein. However, common coding variants exist across unaffected or healthy individuals.[66,69-72] In our study, we analysed the properties of common hCaV3.3 channel variants in the general population that better reflect the range of functional properties of wild-type channels (Fig. 1 and Table 1). Notably, none of the common CACNA1I variants we analysed significantly altered any of the four biophysical parameters (Fig. 1), consistent with reports that no association between CACNA1I common coding variants and schizophrenia risk.[5] By analysing a comprehensive allelic series of variants in this disease cohort, we discovered that CACNA1I variants from unaffected individuals produced a range of altered properties, but these were distinct from variants derived from the patients (Figs 1, 2 and 5). These findings emphasize the importance of analysing rare variants found in the common and ‘unaffected’ population to parse out disease-associated phenotypes.

We found that ultrarare amino acid variants with mutations located in the pore domain along transmembrane S5 and S6 segments were associated with reduced hCaV3.3 channel protein levels, suggesting that missense changes in the pore domain influence the overall stability and/or expression of the channel (Fig. 3D and E). We also identified three regions in hCaV3.3 that contained amino acid changes that influenced the voltage dependence of channel activation (Fig. 3A and B). While our findings are overall consistent with the known impact of mutations in the VSD and pore domain of ion channels,[73-75] such structural–functional relationship analyses in hCaV3.3 channels provided specific mappings that will ultimately help in the development of structure-assisted prediction algorithms[40] to predict CACNA1I variants that are likely to modify V1/2ACT from those that are most likely to alter current density to facilitate large-scale genotype–phenotype analyses of human channelopathies.

Ultrarare variants from schizophrenia patients were more likely to modify the voltage dependence of activation than common variants in the population, whereas those in the unaffected group mostly influence current density (Fig. 2), suggesting that the range of voltages over which hCaV3.3 channels activate and the channel current densities are the most important biophysical phenotypes associated with the genetics of CACNA1I. Subsequently, we found that the V1/2ACT and channel current density of all rare variants were segregated across disease and control status, and alleles that encode CaV3.3 channels with reduced current density were enriched in control subjects (Fig. 5A and Table 3). Corroborating our findings, the data from the SCHEMA consortium[10] showed that rare PTV of CACNA1I is depleted in patients (Fig. 5A and Table 3). These analyses together support the notion that rare cases of haploinsufficiency or reduced channel current density from CACNA1I variants may protect against schizophrenia risk. We previously found that R1346H-CaV3.3, identified in one schizophrenia patient, reduced channel function and sleep spindle density in mice, suggesting that a reduction in CaV3.3 function may underlie the reduced spindle densities observed in schizophrenia patients.[27] These seemingly contradictory roles of CaV3.3 channels in schizophrenia pathophysiology may partly arise from the lack of characterization of sleep spindles as a marker for CaV3.3 function during and prior to the early course of the disease, as most studies focused on chronic patients. While sleep spindle density increases during adolescence in the general population,[76] it is unknown how sleep spindles develop over the course of schizophrenia. A similar paradoxical role exists in SCN2A, where both loss of function and gain of function of SCN2A mutations underlie neurodevelopmental encephalopathies[77] with different clinical symptoms and developmental trajectories. It is also plausible that other genetic or environmental risk factors confound the interpretation of the analyses of rare coding mutations and that distinct subpopulation of schizophrenia patients may have different sleep spindle abnormalities. For example, a subpopulation of patients with increased sleep spindles in the early course of the disease may exist, with such hyperfunctionality ultimately resulting in a reduction in TRN function as a compensatory mechanism. While the exact role of CaV3.3 in schizophrenia pathogenesis needs to be further examined, both functional analyses of rare genetic variants of CaV3.3 channels and translational studies on loss of/reduction in CaV3.3 function models support a role of CaV3.3 in schizophrenia pathophysiology.

We excluded the variants in the C-terminus in our analyses, perhaps limiting the scope of our analyses of 57 variants of CACNA1I in this cohort (Supplementary Fig. 1A and B). The C-terminal region of CaV3.3, beyond the proximal functional motif, has little impact on the channel biophysical properties[39] in HEK293 cells, supporting not analysing these variants in the heterologous expression system. Unstructured, this region likely mediates protein–protein interactions and its functional output is best examined in a neuronal system, such as Neuro2A cells or induced human neurons, where its interaction partners are present. More interestingly, genes such as CACHD1[78] modulate the surface expression of CaV3.3 channels, potentially through interacting with the C-terminal tail of CaV3.3. Identifying the functional interactome of CaV3.3 channels and characterizing the genetic variations of channel modulators and/or interactors will further help examine the disease burden mediated through CaV3.3 channels.

We modelled TRN excitability as a readout for the impact of genetic variants on the physiology of Cav3.3, given the well-described role of CaV3.3 in TRN physiology while its role in other regions of the brain is poorly understood.[19] Variants associated with the loss-of-function phenotype of TRN output (Fig. 5D and E) showed a non-significant trend when compared with neutral or gain-of-function alleles across disease statuses, suggesting that TRN rebound output may be a candidate physiological phenotype associated with CACNA1I in schizophrenia risk. The TRN modulates information transfer between the thalamus and cortex, and TRN excitability influences sleep spindles, sensory gating, focused attention and cognitive flexibility.[23,27,30,79-84] Previous studies showed that reduced CaV3.3 function led to reduced TRN function, with impaired sleep spindles and attention.[27,81] A gain-of-function phenotype of TRN rebound bursting may exert stronger inhibition on thalamocortical cells, disrupting thalamocortical processing, as reported in Dravet syndrome.[85] Recently, gain-of-function mutations in CACNA1I have been associated with neurodevelopmental disorders such as intellectual disability and epilepsy,[86] supporting our analyses that loss of function in CACNA1I may reduce the disease burden in schizophrenia and potentially other neurological disorders. These analyses provide evidence supporting the development of CaV3.3-specific inhibitors that may benefit patients with gain-of-function mutations in CACNA1I and other phenotypes associated with such mutations.

In conclusion, we defined the biophysical properties associated with an allelic series of rare missense CACNA1I variations. Our results showed that a reduction in the current density of hCaV3.3 is enriched in control individuals, while a shift in the voltage dependence of activation may be a candidate biophysical phenotype associated with schizophrenia patients in a Swedish cohort of >10 000 individuals. Our study is the first to reveal a clear functional correlation between rare CACNA1I variants and disease burden in a schizophrenia cohort. The power of our study depends on analysing a series of 50+ CACNA1I rare variants and defining the phenotypes of variants from both affected and unaffected individuals. Additional analyses in a larger cohort may further validate the relationship between CACNA1I and disease risk. Our approach may serve as a template strategy for channelopathies in polygenic disorders, highlighting the importance of analysing coding allelic series of risk genes in common complex disorders.

Click here for additional data file.