The genetic and phenotypic correlates of neonatal Complement Component 3 and 4 protein concentrations with a focus on psychiatric and autoimmune disorders

The iPSYCH2012 study

Key elements of this study were based on the Lundbeck Foundation Initiative for Integrative Psychiatric Research (iPSYCH) sample²⁶, a population-based case-cohort designed to study the genetic and environmental factors of schizophrenia (SCZ), bipolar disorder (BIP), depression (DEP), autism spectrum disorder (ASD) and attention deficit hyperactivity disorder (ADHD). The original iPSYCH sample (known as iPSYCH2012) included information on case status complete through 31 December 2012. We also included 4,791 anorexia nervosa cases (AN; ANGI-DK) from the Anorexia Nervosa Genetics Initiative (ANGI)²⁹, which had the same design as iPSYCH2012. Henceforth, we refer to iPSYCH2012 as the combined dataset with the ANGI samples. The iPSYCH2012 sample is nested within the entire Danish population born between 1981 and 2005 (n=1,472,762). Diagnoses were identified in the Danish Central Psychiatric Research Register^30,31, which includes all inpatient contacts in Danish psychiatric hospitals since 1969 and all outpatient and emergency contacts since 1995. The ICD-10 codes used to classify the psychiatric disorder cases can be found in Supplementary Table 1. The phenotype information for the iPSYCH2012 participants was updated for the target mental disorders until December 2016. The case-cohort sample includes a population-based random sub-cohort³² (N = 30,000) with an inclusion probability of 2.04% of the study base (30,000 / 1,472,762). This sub-cohort also includes some participants with the target mental disorders of interest. The genotypes and C3 and C4 protein concentrations were measured in neonatal dried bloodspots (DBSs) taken as part of routine screening at birth from all babies born in Denmark since 1981 and stored in the Danish Neonatal Screening Biobank³³. Dried blood spot samples have been collected from practically all neonates born in Denmark since 1st May 1981 and stored at −20 °C. Samples are collected 4–7 days after birth. After the dried blood spots were retrieved from the biobank, samples were extracted in a PBS buffer and stored for further use at - 80°C. Subsequently, DNA was extracted according to previously published methods³⁴. After storage the protein extracts were assayed for C3 and C4 concentrations. Thus, all genotypes and C3/C4 protein concentration data originated from a single DBS extraction. Additional details related to blood spot extraction and storage are provided in Supplementary Methods 1.

Ethical framework

Material from the Danish Neonatal Screening Biobank has been used primarily for screening for congenital disorders, but are also stored for follow-up diagnostics, screening, quality control and research. According to Danish legislation, material from The Danish Neonatal Screening Biobank can be used for research after approval from the Biobank, and the relevant Scientific Ethical Committee. There is also a mechanism in place ensuring that one can opt out of having the stored material used for research. The Danish Data Protection Agency and the Danish Health Data Authority approved this study. According to Danish law, informed consent is not required for register-based studies. All data accessed were deidentified.

C3 and C4 protein concentrations

These methods have been described in a related study³⁵. Two 3.2 mm discs of DBS were punched into 96 well polymerase chain reaction plates (72.1981.202, Sarstedt). The extracts were analyzed with a multiplex immunoassay (also measuring vitamin D binding protein³⁵) using U-plex plates (Meso-Scale Diagnostics (MSD), Maryland, US) employing antibodies specific for complement C3 (HYB030-07 and HYB030-06) and complement C4 (MA1-72520 (ThermoFisher Scientific) and HYB162-04). The antibodies were purchased from SSI Antibodies (Copenhagen, Denmark) except if otherwise stated. Extracts were analyzed diluted 1:70 in diluent 101 (#R51AD, MSD). Capture antibodies (used at 10 μg/mL as input concentration) were biotinylated in-house using EZ-Link Sulfo-NHS-LC-Biotin (#21327, Thermo Fisher Scientific) and detection antibodies were SULFO-tagged (R91AO, MSD), both at a challenge ratio of 20:1. As calibrators, we used complement components purified from human: C3: #PSP-109 (Nordic Biosite, Copenhagen, DK), C4: abx060108 (Abbexa, Cambridge, UK). Calibrators were diluted in diluent 101, detection antibodies (used at 1 μg/mL) were diluted in diluent 3 (#R50AP, MSD). Controls were made in-house from part of the calibrator solution in one batch, aliquoted in portions for each plate, and stored at -20°C until use. The samples were prepared on the plates as recommended by the manufacturer and were read on the QuickPlex SQ 120 (MSD) 4 min after adding 2x Read buffer T (#R92TC, MSD). Analyte concentrations were calculated from the calibrator curves on each plate using 4PL logistic regression using the MSD Workbench software.

Intra-assay variations were calculated from 38 measurements analyzed on the same plate of a pool of extract made from 304 samples. Inter-assay variations were calculated from controls analyzed in duplicate on each plate during the sample analysis, 1022 plates in total. Lower limits of detections were calculated as 2.5 standard deviations from 40 replicate measurements of the zero calibrator. The higher detection limit was defined as the highest calibrator concentration. The lower and upper detection limits for: (a) C3 were 95.4 and 7.98×10⁴ μg/L respectively, and (b) C4 were 55.2 and 7.98×10⁴ μg/L respectively. The intra- and inter-assay coefficient of variation (CV) for (a) C3 were 5.2% and 18.1% respectively; and for (b) C4 were 3.9% and 8.5% respectively. To validate the stability of the samples during storage, we randomly selected 15-16 samples from five years (1984, 1992, 2000, 2008, and 2016; a total of 76 samples). After extracting the samples and adding them to an MSD plate, the rest of the extracts were frozen for 2 months, thawed and measured as described above to imitate the freeze-thaw cycle of the samples in the study. The oldest samples (from 1984) recorded higher concentrations, most probably due to a change in the type of filter paper after 1989. In light of this artefact, we adjusted all DBP values by plate (the sequence of testing followed the date of birth of the sample). This is described below. Additional details related to pre-analytic variation are provided in Supplementary Methods 1.

Imputation of genotypes

DNA genotyping was conducted at the Broad Institute (Boston, MA, USA) using the Infinium PsychChip v1.0 array (Illumina, San Diego, CA, USA)³⁶. We restricted the genotyped SNPs to 252,339 high quality, common SNPs based on build hg19 (the same human genome reference build was used throughout this study). Details of the filtering can be found elsewhere (Schork et al., 2019). Briefly, we excluded SNPs with minor allele frequency (MAF) < 0.01, Hardy Weinberg Equilibrium (HWE) p-value < 1.0×10⁻⁶ or non-SNP alleles (i.e., insertions and deletions, INDELs). 245,328 autosomal and 7,011 X-chromosome (chrX) SNPs were retained and used to impute SNPs using the Ricopili pipeline³⁷ with the Haplotype Reference Consortium (HRC)³⁸ as the imputation reference panel (accession number: EGAD00001002729). 6,743,499 autosomal SNPs, 227,371 chrX SNPs for males and 184,517 chrX SNPs for females were retained with missing rate < 0.02 and genotype call probability > 0.8. We further excluded the imputed SNPs with imputation info score < 0.8, MAF < 0.01 or HWE p-value < 1.0×10⁻⁶. 5,201,724 SNPs were retained in autosomes and 126,109 SNPs were retained on chrX. We then used the common SNPs to infer the genetic ancestries of 80,873 participants in the iPSYCH2012 study, 75,764 individuals of European ancestry and 5,109 individuals of non-European ancestry. Details are provided in Supplementary Method 2.

Imputation of C4 haplotypes

C4 haplotypes were imputed from reference data^9,14 using the genotyped SNPs in the iPSYCH2012 sample. The human C4 haplotypes have various copy numbers, including two isotypic polymorphisms, C4A (A) and C4B (B). Each isotype has two length-polymorphisms due to a human endogenous retroviral (HERV) insertion, long form (L, with HERV insertion) and short form (S, without HERV insertion). The isotypic and length polymorphisms lead to four alleles in a C4 copy, AL, AS, BL and BS. Using the genotyped SNPs, the C4 haplotype reference was used to impute the C4 alleles and the number of C4 copies (with a maximum copy number of 4). The C4 haplotype imputation panel comprised whole genome sequencing data from 1,234 individuals of multiple ancestries, which enabled us to identify C4 alleles with high accuracy. We used Beagle software³⁹ for the imputation with the C4 haplotype reference. The imputation results provided the counts of alleles, but were unable to confidently distinguish all combinations of variants, for example, between the haplotypes AS-BL and AL-BS. We counted the two C4 alleles (C4A and C4B) with combination of HERV using a subset of the imputed result, where combinations can be confidently distinguished (details are provided in Supplementary Method 3). Both counts of C4 allele combinations and reported studies⁴⁰ indicated that the C4A gene is more likely to carry HERV insertion than the C4B gene. Therefore, the C4 haplotype is assumed to be AL-BS rather than AS-BL, consistent with methods described by Sekar et al⁹. The imputed counts were converted to the C4 haplotypes. Eight common C4 haplotypes (allele frequencies ≥ 0.01) were imputed in the iPSYCH2012 study (Supplementary Table 3). The allele frequencies of the 8 haplotypes were consistent with other studies^14,41. Given the common C4 haplotypes, we counted the copy numbers of the C4 alleles (Supplementary Figure 2) for each participant. 28 individuals (0.04%) carried 4 copies of C4B and 35 individuals (0.05%) carried 6 copies of HERV insertion. Therefore, we excluded these individuals with very rare copy numbers. The C4A copy number is strongly correlated with C4B and HERV copy numbers (Pearson correlation between C4A and C4B = -0.52; between C4A and HERV = 0.73).

Quality control of the C3 and C4 protein concentrations

The C3 and C4 protein concentrations were measured in 78,268 iPSYCH2012 participants of multiple ancestries. We focused on 68,768 individuals of European ancestry with measures of C3 and C4 protein concentrations. The protein assay plates captured a substantial amount of variance (C3 = 49.4%, C4 = 45.3%). Therefore, we used a linear mixed model (LMM)⁴² approach to adjust protein concentrations, y = ∑z_plateu_plate + e, where y represents the C3/4 protein concentration; z_plate represents protein assay plate, a random variable; u_plate represents the random effect of protein assay plate; and e represents residual. The mixed model regression was conducted by the R package of lme4 (Bates et al., 2015). The rank-based inverse normal transformation (RINT)⁴³ was applied to the residuals to have mean 0 and variance 1. The standard deviations (SDs) adjusted for variance captured by protein assay plate were used for the interpretation of results of C3 and C4 protein concentrations, for C3 protein concentration, 1 SD unit = 2.56 μg/L (3.60 μg/L × √(1 - 0.49)), and for C4 protein concentration, 1 SD unit = 2.46 μg/L (3.33 μg/L × √(1 - 0.45)).

Heritability and SNP-based heritability of the C3 and C4 protein concentrations

The iPSYCH2012 cohort had 75,764 participants of European ancestry. 19,113 participants who shared a genetic relatedness (entry of genetic relationship matrix, r_GRM) ≥ 0.05 with at least one other individual were considered as relatives; 3,253 first degree (r_GRM ≥ 0.4), 2,077 second degree relatives (0.2 ≤ r_GRM < 0.4) and 13,783 third degree relatives (0.05 ≤ r_GRM < 0.2). We jointly estimated both the heritability (h²) and the SNP-based h² (h²_SNP) of the C3 and C4 protein concentrations by using the method proposed by Zaitlen et al⁴⁴. This method assumes a normal distribution of SNP effect sizes. The GWAS studies of protein concentrations^45-47 observed that SNPs in or near the coding genes would capture more phenotypic variance than the remaining SNPs. Therefore, we used two approaches to further explore the h²_SNP, 1) estimating it using all common SNPs by BayesR⁴⁸, which assumes a mixture distribution of SNP effect sizes and can be used to test number of SNPs with nil, small, median and large genetic variance (small R² < 0.01%, median 0.01% ≤ R² < 0.1%, and large 0.1% ≤ R² < 1%) in addition to estimation of h²_SNP, and 2) partitioning h²_SNP into (a) h²_cis-chr, explained by SNPs on the chromosome where the coding gene (cis-chr SNPs) was positioned, and (b) h²_trans-chr, explained by the remaining SNPs (trans-chr SNPs). To partition h²_SNP, we divided SNPs into two subsets, 1) cis-chr and 2) trans-chr SNPs. The h²_cis-chr and h²_trans-chr were estimated by using GCTA-GREML⁴⁹ and BayesR. Only genetically unrelated participants were included in the two analyses. The genetic relationship matrix used in the Zaitlen and GREML analyses were estimated from 5,201,724 common SNPs. Only the subset of HapMap phase 3 (HM3) SNPs were included in the BayesR analyses because of the computation complexity (853,129 HM3 SNP in total; for C3 protein concentration, 132,239 cis-chr SNPs and 839,891 trans-chr SNPs; for C4 protein concentration, 60,631 cis-SNPs and 792,498 trans-SNPs). The Zaitlen method and GREML were implemented in Genome-wide Complex Trait Analysis (GCTA)⁵⁰. BayesR was implemented in Genome-wide Complex Trait Bayesian analysis (GCTB). The URLs for these programs are provided below.

We estimated the genetic correlation between C3 and C4 concentrations by BOLT-REML⁵¹. To further examine if the genetic correlation was primarily driven by the two protein-coding genes (i.e. C3 and C4), we conducted the BOLT-REML and Haseman-Elston regression⁵² analyses using the trans-chr SNPs. The Haseman-Elston regression was implemented in GCTA.

GWAS of C3 and C4 protein concentrations

We performed the GWAS analysis of the C3 and C4 protein concentrations by fastGWA⁵³. The fastGWA is a LMM method which can include all individuals of European ancestry regardless of relatedness. 5,201,724 imputed SNPs were analysed in the GWAS. In addition, fastGWA can include candidate markers which optimizes power if particular SNPs capture a large proportion of the total variance⁵⁴. Therefore, we excluded the SNPs in and near the coding gene for the required GRM in the GWAS, C3: chr19, 4.67Mb – 8.74Mb, C4: chr6, 24.8Mb – 33.9Mb. For each GWAS, we fitted birthyear, sex, wave (i.e., genotyping batch) and the first 20 PCs as covariates in the model. The PCs were estimated by FastPCA⁵⁵, excluding the same SNPs as we did for the required GRM. We conducted the GWASs using all SNPs on autosomal and sex chromosomes. SNPs on X chromosome for males (coded as 0/2) were tested as diploid, assuming X chromosome of males has half dosage compensation⁵⁶. We used GCTA-COJO⁵⁰ to identify the SNPs which were independently associated with the two concentrations. We randomly sampled 10,000 participants from the population-based sub-cohort of iPSYCH2012 as the LD reference cohort. The GWAS significance threshold was 5.0×10⁻⁸.

To explore if the enrichment of mental disorder cases in the iPSYCH2012 case-cohort could induce bias within the GWASs, we conducted simulations with ascertained individuals and performed GWASs in the population-based sub-cohort (Supplementary Method 4).

Associations between C4 haplotypes and protein concentrations

We examined the associations between the imputed C4 haplotypes and the two observed C3 and C4 protein concentrations. We first examined the associations of C4 copy numbers using a LMM approach, y_protein = x_copyb_copy + ∑x_cb_c + ∑z_-MHCu_-MHC + e, where y_protein was C3/4 protein concentration; x_copy was copy number of C4 allele, either C4A, C4B or HERV; b_copy represents effect of copy number; x_c was covariate with b_c being its effect; Both b_copy and b_c were fixed effects. The covariates in the model were the same as those fitted in the GWAS of C4 protein concentration. Because of the strong linkage disequilibrium in the MHC region, fitting SNPs in the region is likely to underestimate the effects of C4 allele count and reduce power. Therefore, we fitted the SNPs outside the MHC region (z_-MHC) in the model with u_-MHC being their random effects. In practice, the effect of copy number from the linear mixed model could be estimated by generalized least squares (GLS) method, b = (X^TV^-1X)^-1X^TV^-1y_protein, where X = {x_copy, x_c} and V was phenotypic covariance matrix of C3/4 protein concentration. It was implemented by GCTA-GREML. All the individuals of European ancestry were included in the analysis. The three C4 allele counts were correlated. Therefore, we estimated the joint effects using the same LMM approach as above, y_protein = x_C4Ab_C4A + x_C4Bb_C4B + x_HERVb_HERV + ∑x_cb_c + ∑z_-MHCu_-MHC + e. In the model, x_C4A, x_C4B and x_HERV represent respectively counts of C4A, C4B and HERV. The remaining variables were defined as above. Secondly, we further examined the associations of the imputed C4 haplotypes using the LMM approach. Previous studies have reported a strong effect for the C4A gene⁹, while the effect of C4B remains unclear⁵⁷. Therefore, we used ‘BS’ as the reference haplotype to estimate joint effects of the remaining haplotypes. The regression model can be expressed as, y_protein = ∑x_alleleb_allele + ∑x_cb_c + ∑z_-MHCu_-MHC + e, where x_allele represents C4 haplotype. Seven C4 haplotypes were included in the model, except for BS. The remaining parameters were defined as above. The estimated effect can be interpreted as the effect of the C4 haplotype compared to BS. All the European participants were included in the analysis. The significance threshold for these analyses was the same as the main GWAS significance threshold (i.e., 5.0×10⁻⁸).

FUMA and SMR

We conducted gene-based analysis by Functional Mapping and Annotation of Genome-Wide Association Studies (FUMA)⁵⁸. There were 18,305 genes available for the gene-based analysis, thus the Bonferroni corrected threshold was 1.4×10⁻⁶ (= 0.05 / (18,305 × 2)). We conducted Summary-data-based Mendelian Randomisation (SMR)⁵⁹ to identify pleiotropic genes for C3 and C4 concentrations. For the SMR analysis, the eQTL data (i.e., summary statistics from associations of gene expressions), was Genotype-Tissue Expression version 8 (GTEx v8)⁶⁰. The LD reference sample with 10,000 participants was the same as the above GCTA-COJO analysis. 22,338 gene-tagged probes within 49 tissues (200,144 probes in total) which had significant SNPs were included in the SMR analysis. The Bonferroni significance threshold was 1.2×10⁻⁶ (= 0.05 / (200,144 × 2)).

Associations between C4 haplotypes and mental disorders observed within the iPSYCH2012 case-cohort study

Based on the associations with protein concentrations, we conducted the associations between C4 haplotypes and 6 iPSYCH2012 disorders (SCZ, BIP, DEP, ASD, ADHD and AN). We used three approaches to examine the relationships, 1) associations with C4 allele counts, 2) associations with imputed C4 haplotypes, 3) associations with predicted C4 gene expression in the brain. Because the iPSYCH2012 case-cohort study has person-level data on the age-at-first contact with psychiatric services, we were able to assess the risk of mental disorders within the more informative time-to-event framework, using Cox proportional hazards regression (Cox PH) to analyse the hazards of C4 allele counts and haplotypes with respect to the mental disorder of interest. For C4 allele count, we examined the joint effects due to their correlations, h(t) = h₀(t)exp(x_C4Ab_C4A + x_C4Bb_C4B + x_HERVb_HERV + ∑x_cb_c). In the model, h₀(t) represents the baseline hazard while h(t) represents the hazard at time t between baseline and December 2016. The remaining variables were defined as above. For C4 haplotypes, we examined the joint effects using the Cox PH model, h(t) = h₀(t)exp(∑x_alleleb_allele + ∑x_cb_c). All the variables were defined as above. In addition to C4 haplotypes, we used the predicted C4A and C4B gene expressions as outlined in the post-mortem brain study of Sekar et al.⁹ The association was conducted with a Cox PH model, h(t) = h₀(t)exp(x_{C4A_predicted}b_{C4A_predicted} + x_{C4B_predicted} b_{C4B_predicted} + ∑x_cb_c), where x_{C4A_predicted} and x_{C4B_predicted} represent the predicted C4A and C4B gene expressions, respectively. We included only unrelated individuals of European ancestry in the three analyses.

In the time-to-event analysis, the cases were the diagnosed participants by December 2016, and the non-cases are defined as the entire cohort excluding those individuals with the disorder of interest. Therefore, we defined six psychiatric-disorder samples for the time-to-event analyses. The sample sizes for cases and non-cases are shown in Supplementary Table 1.

Associations between protein concentrations and mental disorders observed within the iPSYCH2012 case-cohort study

Based on the associations between C4 haplotypes and 1) the two protein concentrations (C3 and C4) and 2) six mental disorders, we explored the associations between C3 and C4 protein concentrations and mental disorders observed in the iPSYCH2012 case-cohort study, using Cox PH models. We first tested the marginal effects of two concentrations, using the model h(t) = h₀(t)exp(x_proteinb_protein + ∑x_cb_c). The variables were defined as above. Due to the high correlation, we then fitted both concentrations jointly, h(t) = h₀(t)exp(x_{C3_protein}b_{C3_protein} + x_{C4_protein}b_{C4_protein} + ∑x_cb_c), where x_{C3_protein} and x_{C4_protein} represent C3 and C4 concentration, respectively. b_{C3_protein} and b_{C4_protein}, effects of two protein concentrations, were fixed effects. In the three analyses, we included only unrelated individuals of European ancestry.

Mendelian Randomisation analysis based on summary statistics

We explored the relationships between protein concentrations and mental and autoimmune disorders using the generalised summary-data-based Mendelian Randomisation (GSMR) method⁶¹. Because of the possible link between C3 and C4 and brain function⁷, we also included two neurodegenerative disorders—Alzheimer’s disease, and amyotrophic lateral sclerosis in these analyses. Thus, there were 8 broadly-defined neuropsychiatric disorders (i.e., SCZ⁶², DEP⁶³, BIP⁶⁴, ASD⁶⁵, ADHD⁶⁶, AN⁶⁷, Alzheimer’s disease⁶⁸, amyotrophic lateral sclerosis⁶⁹), and 6 autoimmune disorders (i.e., multiple sclerosis⁷⁰, type-1 diabetes⁷¹, Crohn’s disease⁷², ulcerative colitis⁷², rheumatoid arthritis⁷³, systemic lupus erythematous⁷⁴). The GWAS summary statistics for these disorders were publicly available (additional details provided in Supplementary Table 19).

Unfortunately, detailed GWAS summary statistics for Sjögren’s syndrome were not available. The GSMR method was implemented in GCTA. The GSMR method includes options to exclude potentially pleotropic loci (HEIDI filtering). Similar to the methods for pleiotropy exclusion, it assumes fewer pleiotropic SNPs than valid variants and MR estimates may fluctuate with the inclusion of those SNP outliers⁶¹. Given that a substantial proportion of the variance in C4 concentration was associated with SNPs in the MHC region (including the C4 gene), genetic variants used in the GSMR analysis may be dominated by those positioned in or near the region. Therefore, the interpretation of results based on the large effect loci within the MHC can be misleading. As a planned analyses, we ran the GSMR method again based on summary statistics from the GWAS of C4 concentration adjusted for COJO SNPs within MHC region (as a covariate). In addition, in the presence of bidirectional GSMR findings (i.e., evidence of both protein concentration impacting on phenotype of interest, and vice versa; forward and reverse direction respectively), symmetrical (i.e., equivalent) effect sizes may reflect the presence of pleiotropy. It is of note that effects of C3/4 protein concentration on these mental and autoimmune disorders by GSMR were equivalent to log odds ratio (OR) from logistic regression if we have all the phenotypes in a cohort⁶¹. We, therefore, reported the logOR and 95% CI in the results. The effects of mental and autoimmune disorders on C3/4 protein concentration may be required to be transformed for better interpretation⁷⁵. However, we only used the reserve GSMR results to examine the presence of reverse causation and pleiotropy. Therefore, we reported the GSMR estimates directly. The LD reference sample which included 10,000 participants was the same as for the GCTA-COJO analysis. We conducted both forward (from C3/4 concentration to disorder) and reverse (from disorder to C3/4 concentration) analyses. We set HEIDI-outlier threshold at 0.01 to filter horizontal pleiotropy. The GSMR Bonferroni corrected significance threshold was 1.9×10⁻³ (= 0.05 / (2 × 13)).

PheWAS based on UK Biobank phenotypes

Based on the GSMR analysis results, we conducted phenome-wide association studies (PheWASs) to explore the relationships with disorder outcomes. The PheWASs were conducted in the UK Biobank (UKB) cohort⁷⁶, a large population cohort with 487,409 participants of multiple ancestries. The genotypes were imputed to the HRC³⁸ and UK10K⁷⁷ reference panels by the UKB group. The quality controls were described in detail elsewhere⁷⁸, including generic ancestry determination, quality controls of imputed SNPs, and estimation of principal components. In the study, we included 1,130,559 HM3 SNPs on autosomal chromosomes, with MAF ≥ 0.01, HWE P-value ≥ 1.0 ×10⁻⁶, because only effects of HM3 SNPs were predicted by BayesR. The genetic relationship matrix was estimated by GCTA. 347,769 unrelated participants of European ancestry were retained with genetic relationship < 0.05. In the PheWAS analysis, we included 1,148 UKB phenotypes, 1) 1,027 disorders which were classified by ICD-10 codes, 2) 51 anthropometric measurements and brain imaging traits, and 3) 70 infectious disease antigens. The quantitative traits were standardised by RINT to have mean 0 and variance 1. We then used the model to test the associations, for quantitative traits, y = x_{protein_prs}b_{protein_prs} + ∑x_cb_c + e, where y represents quantitative trait in UKB; x_{protein_prs} represents polygenic scores for neonatal C3 and C4 protein concentrations predicted by BayesR; x_c represent the covariate variables including birth year, sex and 20 PCs. For dichotomous traits, logit(y) = x_{protein_prs}b_{protein_prs} + ∑x_cb_c + e, where y represents the dichotomous trait and definitions of the remaining variables were the same as above. In addition, we conducted the PheWAS analyses for males and females separately using the same approach. Polygenic scores were predicted using GWASs in both sexes. The Bonferroni corrected significance threshold was 7.3×10⁻⁶ (= 0.05 / (1148 × 3 × 2)).

Heritability of C3 and C4 protein concentrations

The h² of C4 by Zaitlen’s method⁴⁴ was 40% (SE = 0.03, P-value = 2.7×10⁻⁴⁴, Supplementary Table 5) while the h²_SNP was 26% (SE = 0.006, P-value < 1.0×10⁻¹⁰⁰). For C3, h² was 21% (SE = 0.03, P-value = 1.1×10⁻¹¹) and the h²_SNP was 4% (SE = 0.005, P-value = 3.2×10⁻¹⁴). The high genetic variance of C4 concentration was confirmed by BayesR⁴⁸ — the h²_SNP for C4 was 24% (SE = 0.004, P-value 1.0×10⁻¹⁰⁰ and 6% (SE = 0.005, P-value = 4.4×10⁻³⁹) for C3. Moreover, ∼50 HM3 SNPs captured substantial genetic variance (R² > 0.1%), for C3 protein concentration 39 SNPs and for C4 protein concentration 62 SNPs. The observation was consistent with our expectation that SNPs in or near the salient encoding genes would capture a substantial proportion of phenotypic variance for these two circulating proteins. Thus, we partitioned SNPs into two subsets; (a) those on the chromosome where the coding gene is located (cis-chr SNPs), and (b) those on the remaining chromosomes (trans-chr SNPs), and jointly estimating SNP-based h² at the two subsets of SNPs, SNP-based h²_cis-chr and h²_trans-chr, using GREML⁴⁹. For C4, h²_cis-chr = 14% (SE = 0.005, P-value < 1.0×10⁻¹⁰⁰) and h²_trans-chr = 4% (SE = 0.006, P-value = 9.4×10⁻¹³). For C3, h²_cis-chr was 0.4% (SE = 0.001, P-value = 1.1×10⁻³) and h²_trans-chr = 3% (SE = 0.006, P-value = 1.4×10⁻⁷). The cis-chr SNPs of C4 concentration captured more genetic variance than the trans-chr SNPs. A higher genetic variance was found at cis-chr SNPs by BayesR. For C4, h²_cis-chr = 22% (SE = 0.003, P-value < 1.0×10⁻¹⁰⁰) and h²_trans-chr = 5% (SE = 0.004, P-value = 6.2×10⁻³⁰). For C3, h²_cis-chr = 2% (SE = 0.001, P-value = 3.9×10⁻⁸⁰) and h²_trans-chr = 5% (SE = 0.004, P-value = 5.9×10⁻²⁶). Both the Zaitlen method and GREML assume genetic variance at each SNP follows a normal distribution while BayesR assumes a mixture distribution. The BayesR estimates are expected to better model the true underlying genetic architecture. We only used individual-level-data-based methods (such as Zaitlen’s method, GREML and BayesR) in the analysis as genotypes and phenotypes were all available, since summary-level-data-based methods modelling LD from reference cohorts are less accurate. In summary, both C3 and C4 concentrations were heritable traits. SNPs positioned in the C4 gene accounted for a substantial proportion of the genetic variance of C4 concentration.

The genetic correlation (r_g) between the two concentrations based on BOLT-REML⁵¹ was 0.38 (Supplementary Table 6, SE = 0.03, P-value = 1.9×10⁻³⁵), smaller than the phenotypic correlation (0.65). Given the large genetic effects at the coding genes for C3 and C4 we then estimated r_g using SNPs other than chromosomes 6 or 19 (related to the location of C4 and C3 genes, respectively) to further investigate if the correlation was driven by cis-chr SNPs, r_g = 0.82 (SE = 0.05, P-value = 4.8×10⁻⁶⁵). The high correlation was confirmed by Haseman-Elston regression⁵² (r_g = 0.78, SE = 0.19, P-value = 4.1×10⁻⁵), using trans-chr SNPs. These results indicated that C3 and C4 were genetically correlated, and this genetic correlation was not driven by SNPs in or near their respective encoding genes.

GWAS of C3 and C4 protein concentrations

We used fastGWA⁵³ to conduct the GWAS analysis based on 75,764 participants of European ancestry and 5,327,833 common SNPs, 5,201,724 in autosomes and 126,109 on the X chromosome (Figure 2). We conducted a GCTA-COJO⁷⁹ analysis to help identify putative independent SNPs. For C4 protein concentration, 34 autosomal SNPs were identified as genome-wide significant (Supplementary Table 7), and all were autosomal. For C3 protein concentration, 14 significant SNPs were identified, again all were autosomal (Supplementary Table 8).

• Download figure
• Open in new tab

Figure 2

GWASs of neonatal C4 and C3 protein concentrations. a) unadjusted C4 protein concentration, b) C4 protein concentration adjusted for COJO SNPs in the MHC region (fitted as covariates in the regression model), c) unadjusted C3 protein concentration and d) C3 protein concentration adjusted for COJO SNPs on chromosome 19. The COJO SNPs fitted as covariates in GWAS of adjusted protein concentration (panels b and d) were identified from GCTA-COJO analysis of unadjusted protein concentration. The COJO SNPs were highlighted with red colour. The top-associated SNPs were annotated with their overlapped or nearest genes. The GWAS threshold was 5.0×10⁻⁸.

Of the 34 SNPs significantly associated with C4 protein concentration, 30 (88.2%, 30/34) were found on chromosome 6. Of these 29 were in the MHC region and 27 (79.4%, 27/34) SNPs were positioned within 2Mb of the C4 gene (chr6, 31.9 Mb). These 27 SNPs explained 16.7% of phenotypic variance in C4 concentration, which is consistent with the estimated h²_cis-chr. SNP rs113720465 (32,005,355bp, ∼1Kb away from C4B-AS1 [32,000-32,004Kb]) had the largest effect size (the A allele was associated with an increase of 0.76 standard deviation units of C4 protein concentration), however SNP rs3117579 had the smallest P-value (within an exon of GPANK1). Given this large effect size, it is possible that SNPs in LD at r² < 0.01 (the COJO threshold of independence) could also be reported as genome-wide significant through correlation. Thus, we conducted a GWAS fitting the COJO SNPs in and near MHC region as fixed effects (Figure 2). We identified 8 significant loci by COJO, 6 of which were significant from GWAS of unadjusted C4 protein concentration. The 2 additional loci were on chromosomes 9 (rs6477754) and X (rs12012736). Interestingly, nearly all the 8 COJO SNPs were annotated to the genes biologically related to complement-related pathways (Supplementary Figure 5). For example, C4BPA (rs12057769) encodes a binding protein of C4. The IL6 gene (rs2066992) encodes a cytokine stimulated in response to infections and injuries. C1S (7.1Mb on chr12) and C1R (7.2Mb on chr12), the nearest genes of rs11064501, are the protein-coding genes of two C1 sub-components.

With respect to C3 protein concentration, 7 COJO SNPs were positioned within 2Mb of the C3 gene (chr19, 6.7Mb) — these loci explained 3% of phenotypic variance in C3 concentration. After fitting these 7 COJO SNPs as covariates, 8 significant COJO SNPs were identified (Supplementary Figure 6). We found a SNP within the ABO gene, which has recently been identified as a ‘master regulator’ of plasma protein concentration^46,80. The gene annotations of the remaining SNPs encode proteins which involve immune- and/or C3-related pathways: (a) FCGR2B (rs844), which encodes an inhibitory receptor for the Fc region of immunoglobulin gamma (IgG), (b) CFH (rs558103 and rs11580821) which encodes Complement Factor H, a key factor that inhibits the alternative pathway and the amplification loop downstream from C3, (c) STK19 gene (rs114492815) which is close to the C4A gene, and (d) FAM117A (rs12949906), which has enhanced gene expressions in dendritic cells (i.e., antigen-presenting cells involved in the immune system⁸¹.

For both GWASs of C3 and C4 protein concentration, we found no evidence of potential ascertainment bias related to the enrichment of cases with mental disorders in the iPSYCH2012 case-cohort study (Supplementary Method 4, Supplementary Figures 8-10). Therefore, the following post-GWAS analyses were based on the results from the full iPSYCH2012 sample.

Associations between C4 haplotypes and C3/4 protein concentration

Both h²_SNP and GWAS results indicated strong effects of the SNPs in the MHC region for both C4 and C3 protein concentrations. Due to the complex LD structure in this region, we used the imputed C4 haplotypes to investigate phenotypic associations of these genetic variants. We first examined the associations between the imputed C4 haplotypes with the observed C4 protein concentration, using an LMM approach. As expected, more copies of C4 allele (either C4A, C4B, with or without HERV) were strongly associated with higher C4 protein concentration (Figure 3). The C4A copy number (b_C4A) had greater effect than C4B (b_C4B) and HERV (b_HERV), b_C4A = 0.3 (Supplementary Table 9, SE = 0.01, P-value < 1.0×10⁻¹⁰⁰), b_C4B = 0.2 (SE = 0.01, P-value < 1.0×10⁻¹⁰⁰), and b_HERV = 0.2 (SE = 0.004, P-value < 1.0×10⁻¹⁰⁰). The C4 copy numbers were correlated. Therefore, we fitted all 3 gene copy numbers in a regression model to estimate the joint effects. The C4A copy number had nearly identical effect to C4B copy number, b_C4A = 0.6 (SE = 0.01, P-value < 1.0×10⁻¹⁰⁰), b_C4B = 0.6 (SE = 0.01, P-value < 1.0×10⁻¹⁰⁰). The beta estimates associated with the HERV copy number was less than the comparable estimates for C4A and C4B, and was negatively associated with C4 protein concentration, b_HERV = -0.08 (SE = 0.005, P-value = 5.0×10⁻⁵¹). This may reflect the strong correlation with C4A (r = 0.73) and negative correlation with C4B (r = -0.17). The result suggested 1 more copy of C4A or C4B is likely to have 1.6 μg/L (∼0.6 × SD unit) higher C4 protein concentration given the same amount of HERV. We calculated the captured variance (R² = s²b²) that were comparable between the C4 copy numbers. In the formula, s² was variance of C4 copy number, analogous to variance of allele count. Of interest, the s² of C4A count was greater than C4B count (Supplementary Table 3, s² = 0.55 for C4A and 0.31 for C4B). Therefore, C4A count had a larger contribution to C4 protein concentration than the C4B count (R² = 23% for C4A and 11% for C4B). In total, both counts captured 17.3% of variance in C4 protein concentration, accounting for the negative correlation between the two allele counts (r = -0.52). The captured genetic variance was in line with h²_cis-chr, and the genetic variance at the MHC SNPs. In summary, the imputed counts of both C4A and C4B were associated with the observed C4 protein concentration and C4A count had a greater contribution than C4B count.

• Download figure
• Open in new tab

Figure 3

Plot of C4 copy number versus C4 protein concentration. There were three C4 alleles, (a) C4A, (b) C4B and (c) C4L/HERV. The colours represent C4 allele counts.

Based on the effects of C4 allele counts, we then examined the association between the commonly observed C4 haplotypes and C4 protein concentration. For this analysis, we used the BS haplotype as the reference category (because of [1] the positive association between C4 allele count and C4 protein concentration, and [2] the greater contribution of C4A count to C4 protein concentration).

All the remaining 7 common haplotypes were associated with increased C4 protein concentration. Due to their higher frequencies, AL-BS and AL-BL haplotypes captured greater variance of C4 protein concentration (Supplementary Table 10, R² = 7.6% for AL-BS and 7.1% for AL-BL).

We then examined C3 protein concentration. In keeping with expectations, we did not identify any significant associations between either C4 copy number or C4 haplotype, versus C3 protein concentration. However, we found the AL-BS haplotype was nominally significantly associated with C3 protein concentration (b_AL-BS = 0.23, SE = 0.03, P-value = 5.4×10⁻³). From the GWAS of C3 concentration, there was a COJO SNP positioned within the MHC region (rs114492815). While this SNP was in very weak association with each of the C4 allele counts (R² < 0.005 for C4A count and C4B count, 0.01 for C4L/HERV count), there was a moderate association with AL-BS (R² = 0.11). Therefore, we ran the analysis again, fitting rs114492815 as an additional covariate. After this adjustment, none of the haplotypes were associated with the C3 concentration (Supplementary Table 10). In general, C3 concentration was independent of C4 alleles.

Functional mapping of GWAS

Having found the significant SNPs from the GWASs, we explored the genes associated with both concentrations. For C4 concentration, we identified 263 significant genes by Functional Mapping and Annotation of Genome-Wide association Studies (FUMA), 257 (98%) on chromosome 6 and five (2%) on the remaining chromosomes (Supplementary Table 11). These findings are consistent with the high LD between loci in the MHC region and the high gene density in this region. Some differentiation between genes was achieved by using SMR which integrates the trait associations with significant eQTL associations. For SMR, using the eQTL summary data GTEx version 8, we identified 56 pleiotropic genes, 55 on chromosome 6 and one on chromosome 1 (Supplementary Table 12). We noted that the number of identified genes by SMR was smaller than by FUMA. Many genes on chromosome 6 were significant on the SMR test but failed to pass HEIDI test (i.e., there was evidence of pleiotropy). This was because of the complex LD and likely multiple causal alleles. Interestingly, SMR analysis found that C4A, C4B and C4BPA were all significantly associated with neonatal C4 protein concentration. These are three of the major genes involved in regulation of C4 protein concentration. The genetic correlates of neonatal C4 protein concentration were associated with higher C4A gene expressions in 8 brain tissues, amygdala (b_XY = 0.70, SE = 0.11, P-value = 3.7×10⁻¹⁰), anterior cingulate cortex (b_XY = 0.74, SE = 0.13, P-value = 6.8×10⁻⁹), caudate basal ganglia (b_XY = 0.70, SE = 0.09, P-value = 2.6×10⁻¹⁶), cerebella hemisphere (b_XY = 0.44, SE = 0.05, P-value = 8.2×10⁻¹⁸), brain cerebellum (b_XY = 0.43, SE = 0.04, P-value = 5.1×10⁻²²), hippocampus (b_XY = 0.68, SE = 0.10, P-value = 8.0×10⁻¹²), hypothalamus (b_XY = 0.79, SE = 0.11, P-value = 4.4×10⁻¹²), and putamen basal ganglia (b_XY = 0.76, SE = 0.12, P-value = 8.6×10⁻¹⁰). Except for cerebellar hemisphere and cerebellum, the effect sizes of these associations were comparable with a mean of b_XY = 0.73. Overall, the findings from FUMA and SMR indicate strong associations between genes in the MHC region and C4 protein concentration, and the C4A gene was likely to have causal effect on C4 protein concentration in brain tissues. Interestingly, these significant genes were enriched with the Kyoto Encyclopedia of Genes and Genomes (KEGG) gene-sets of systemic lupus erythematosus (SLE, P-value = 1.1×10⁻⁷⁰) and complement systems (P-value = 7.2×10⁻⁴) (Supplementary Figure 7). For C3 protein concentration, we identified 19 genes by FUMA (Supplementary Table 13). Within this set of genes, the DXO gene (chr6: 31.9Mb), which is positioned close to rs114492815, a significant SNP from the C3 GWAS, passed both SMR and HEIDI (a test of pleiotropy) tests (Supplementary Table 14). All these genes by FUMA and SMR were enriched with KEGG gene-sets of SLE (P-value = 3.5×10⁻⁶), leishmania infection (P-value = 5.0×10⁻⁶), complement systems (P-value = 9.1×10⁻⁶) and Fc gamma R-mediated phagocytosis (P-value = 9.5×10⁻⁴) (Supplementary Figure 7). In general, the identified genes (especially those within MHC region) suggest an association between SLE and the complement systems.

Associations with mental disorders within the iPSYCH2012 case-cohort study

We did not identify significant associations between C4A, C4B, or HERV copy numbers (Supplementary Table 15 or C4-related haplotypes (Supplementary Table 16) and any of the six mental disorders. Based on the formula between imputed C4 haplotypes and observed C4 gene expression (i.e., RNA concentration) in post-mortem brain tissue⁹, we found no significant associations between these estimates and any of the 6 mental disorders (Supplementary Table 17). Importantly, we did not find any significant associations between observed neonatal C4 protein concentration and any of the six mental disorders. In models that accounted for the strong correlation between C3 and C4 concentration, we found no significant association between C3 concentration and any of the six mental disorders. (Supplementary Table 18).

GSMR relationships with candidate neuropsychiatric and autoimmune disorders

We conducted Mendelian randomization analyses to examine relationships between the two protein concentrations (C3 and C4) and neuropsychiatric and autoimmune disorders (Supplementary Table 19). The results are shown in Figures 4 and 5. In the unadjusted analysis (i.e., all loci including the MHC region; with and without HEIDI filtering), higher C4 protein concentration was found to be associated with three mental disorders (Figure 4 and Supplementary Figure 11, SCZ, DEP, BIP). The odds ratios for these three findings were small (1.05 or less). We found that these GSMR results were strongly dependent on SNP instruments that were in and near the MHC region (e.g., for schizophrenia and bipolar disorder 126 out of 130 SNPs, and for major depression 103 out of 107 SNPs). These three findings did not persist in the analyses adjusted for the MHC region SNPs. Overall, these analyses do not lend weight to the hypothesis that C4 is causally related to the risk of the psychiatric disorders included in the analyses.

• Download figure
• Open in new tab

Figure 4

Association between C4-related measures and mental disorders, and GSMR analyses with mental disorders. There were 6 mental disorders in the analyses: SCZ, DEP, BIP, ASD, ADHD and AN. The results shown in the top row were from time-to-event analyses between mental disorders and C4-related genotypes and phenotypes, including a) C4 allele counts, b) imputed C4 haplotypes, c) predicted C4 gene expressions and d) C3 and C4 protein concentrations. The analyses were conducted in the iPSYCH2012 cohort. The results from Mendelian randomisation analyses (conducted by GSMR) were shown in the bottom row. The GSMR analysis using GWAS summary statistics predicted relationships between C4 protein concentration and mental disorders, e) using genetic variants from GWAS of C4 protein concentration, f) using genetic variants from GWAS of C4 protein concentration adjusted for COJO SNPs in the MHC region. The Bonferroni corrected thresholds were provided in Methods. Bars represent 95% confidence interval. Significant results were highlighted with “*”.

• Download figure
• Open in new tab

Figure 5

GSMR analyses of C4 protein concentration against autoimmune disorders. The top two panels showed the Mendelian randomisation results a) using genetic variants from GWAS of C4 protein concentration, b) using genetic variants from GWAS of C4 protein concentration adjusted for COJO SNPs in the MHC region. The Bonferroni corrected threshold was 1.8×10⁻³. Bars in the top 2 panels represented 95% confidence interval. Significant results were highlighted with “*”. The lower five panels show the effects of genetic variants for C4 protein concentration without adjustment against effects for five autoimmune disorders. Significant findings were identified for multiple sclerosis, Crohn’s disease, type-1 diabetes, rheumatoid arthritis and systemic lupus erythematosus. The estimates in the panels (OR, standard error and P-value) were estimated from GSMR. Slope of dash line represented logOR. Bars in the bottom 5 panels represented standard errors of SNPs from GWAS. Potential pleiotropic SNPs were excluded.

In contrast, we found strong, protective effects of C4 concentration for several autoimmune disorders. Higher C4 concentration was associated with lower risks of multiple sclerosis, type-1 diabetes, rheumatoid arthritis and systemic lupus erythematosus (Figure 5). The effects were very large, especially for type-1 diabetes (OR = 0.54, 95% confidence interval [CI] = 0.50 – 0.58, N_SNP = 47) and systemic lupus erythematosus (OR = 0.37, 95% CI = 0.34 – 0.42, N_SNP = 103) (Supplementary Table 20). We identified that higher C4 concentration increased the risk of Crohn’s disease (OR = 1.26, 95% CI = 1.19 – 1.34, N_SNP = 86). The strong effect of neonatal C4 protein concentration on these disorders were not caused by reverse causation (Supplementary Table 21 and Supplementary Figure 12). After removing the pleiotropic SNPs, genetic variants associated with all autoimmune disorders (except for type-1 diabetes) were not associated with C4 protein concentration in the reverse GSMR analysis (from autoimmune disorder to neonatal C4 protein concentration). When we examined the relationships adjusted for the MHC region SNPs, the significant association with SLE persisted. The effect size was comparable to that found using unadjusted C4 GWAS (with adjustment, OR = 0.24, 95% CI = 0.12 – 0.47, N_SNP = 7; without adjustment, OR = 0.37, 95% CI = 0.34 – 0.42, N_SNP = 103). Overall, these findings further support the hypothesis that higher C4 protein concentration is causally related to a reduced risk of systemic lupus erythematosus—it is predicted that an increase of 2.46 μg/L (1 SD unit) of C4 concentration would be associated with a 76% reduced risk (1 – 0.24) of systemic lupus erythematosus.

We then explored the relationships between C3 concentration and neuropsychiatric and autoimmune disorders by bidirectional GSMR. Mindful that analyses based on fewer instruments may be underpowered to detect small effects, no significant associations were identified with pleiotropic SNPs removed (Supplementary Table 20). Our findings provide no support for the hypothesis that C3 protein concentration is related to the risk of the neuropsychiatric nor autoimmune disorders examined in this study.

C3 and C4 phenome-wide association studies in the UK Biobank

With respect to C4 concentration, the PheWAS study in the UK Biobank identified significant associations with 35 phenotypes (Supplementary Table 22 and Supplementary Figure 13). Many of these were related to autoimmunity. Of the top 8 disease associations ranked by P-value, higher C4 concentration was associated with a reduced risk of six disorders, two were associated with an increased risk. The top 8 were intestinal malabsorption (which includes coeliac disease; ICD10 = K90, OR = 0.54, 95% confidence interval [CI] = 0.53 - 0.56); thyrotoxicosis [hyperthyroidism] (ICD10 = E05, OR = 0.77, 95% CI = 0.75 - 0.79); hypothyroidism (ICD10= E03, OR= 0.92, 95% CI = 0.90 – 0.93); insulin-dependent diabetes mellitus (ICD10 = E10, OR = 0.80, 95% CIs = 0.78 – 0.83); sarcoidosis (ICD10 = D86; OR = 0.79, 95% CIs = 0.75 – 0.83); psoriasis (ICD10 = L40; OR = 1.08, 95% CIs = 1.06 – 1.10); systemic lupus erythematosus (ICD10 = M32, OR = 0.74, 95% CI = 0.69 - 0.80) and ankylosing spondylitis (ICD10 = M45, OR = 1.22, 95% CIs = 1.16- 1.28). We also found a significant result for multiple sclerosis (ICD-10 = G35, OR = 0.88, 95% CI = 0.84 - 0.92). The attenuated results for the autoimmune disorders that were included in GMSR may be related to the low prevalence of these disorders (and the smaller number of cases in the UKB). In addition to the disorders which met the Bonferroni-corrected P-value threshold (P-value < 7.3×10⁻⁶), we note that Sjögren’s syndrome (ICD10 = M35; OR = 0.95, 95% CI = 0.92 - 0.97) was nominally significant. There were no significant findings between C4 and any neuropsychiatric disorders (ICD10 F codes). No significant difference was found between males and females. Overall, these findings lend weight to the hypotheses that neonatal C4 protein concentration is (a) not associated with the risk of neuropsychiatric disorders, but (b) is associated with reduced risks of several autoimmune disorders. There were no significant associations between C3 and any of the 1,148 phenotypes, which was in line with our GSMR findings (Supplementary Table 23 and Supplementary Figure 14).

Strengths and Limitations of the study

Our study has several strengths. Our sample was nearly 20 times larger than the only other published GWAS of C3 and C4¹⁸. With respect to the hypothesis linking complement to brain development, our complement assays were collected from neonatal samples (versus adult samples). Because the onset of mental disorders such as SCZ is often in the second and third decade of life, our samples are unlikely to be impacted by reverse causation (e.g., smoking may be linked to complement gene expression in the brain¹⁷) and medication effects may impact on post-mortem gene expression studies⁹⁷. With respect to limitations, because our samples were based on neonatal C3 and C4 concentrations, it remains to be seen if the genetic correlates we identified are stable across the lifespan. Also, we used an antibody that has been demonstrated to measure total C4 (i.e., both C4A and C4B), thus we are unable to isolate the concentrations of the two isoforms. The C3 and C4 concentrations in our study were derived from circulating plasma proteins, whereas the concentration of these proteins may vary between organs/tissues and also in response to local tissue activation pathways.

The study by Sekar et al.⁹, which was based on C4 haplotypes imputed from 28,799 schizophrenia cases and 35,986 controls, found that the AL-AL haplotype was associated with an odds ratio of 1.27 compared with the BS allele (i.e. a 27% increased odds). Our study, based on 2,517 cases and 51,799 non-cases, and using a more informative imputation training set⁹, found no significant association for this comparison (HR = 1.14, 95% CI = 0.84 - 1.54; Supplementary Table 16)). Because the Sekar et al. study included more schizophrenia cases in their analyses, it is feasible that our study was underpowered to confidently detect an association between C4 haplotypes and SCZ⁹. However, we had access to the concentrations of protein product of these haplotypes from a very large sample size (n = 68,768), and thus we could estimate the variance of neonatal C4 protein concentration. Based on this observed variance, our study had 80% power to confidently detect a 25% increased risk of SCZ (OR = 1.25) by 2.46 μg/L C4 protein concentration (1 standard deviation unit). Thus, our study had sufficient power to detect the effect size previously identified by Sekar et al.⁹.