Family genetic designs in MoBa provide insights into health and functioning

The Norwegian MoBa cohort

MoBa is a population-based pregnancy cohort study conducted by the Norwegian Institute of Public Health. Participants were recruited from all over Norway from 1999 to 2008. The women consented to participation in 41% of the pregnancies (n = 112,908 recruited pregnancies)^43,54. The cohort includes approximately 114,500 children, 95,200 mothers and 75,200 fathers. The establishment of MoBa and initial data collection were based on a licence granted by the Norwegian Data Protection Agency and an approval from the Regional Committees for Medical and Health Research Ethics (REK). The MoBa cohort is currently regulated by the Norwegian Health Registry Act. The current study was approved by the Regional Committees for Medical and Health Research Ethics (14140 and 2016/1226), and all research was performed in accordance with relevant guidelines and regulations. All mothers and fathers provided written informed consent at recruitment. Mothers consented to participation on behalf of themselves and their children. At 18 years of age, the children become independent participants and receive an information letter about their rights, including how to withdraw. Participation is voluntary, and participants can withdraw their consent at any time. In accordance with REK regulations, individuals who withdraw consent are excluded. As shown in Fig. 1a, MoBa includes questionnaire data collections at many time points, including during the pregnancy, infancy, preschool-age, school-age, adolescence and beyond. Further details about the cohort representativeness and participation in specific waves of data collection are described elsewhere^43,54,55,56. Detailed instrument documentation is available on the MoBa website (https://www.fhi.no/en/ch/studies/moba/). A wide range of national health and administrative registries can be linked to MoBa for further phenotyping without reliance on participant engagement in completing questionnaires (Fig. 1b).

Phenotypic measures

For our exemplar analyses, we included four phenotypes assessed in middle childhood (age 7–10 years), spanning physical health (height), mental health (depression symptoms) and aspects of functioning (sleep duration and educational achievement). Child height at 7 years of age was reported by mothers in the age-7 questionnaire for the question, ‘What is the child’s height and weight now at 7 years old?’ Mothers were asked to report their child’s current height in centimetres. Sleep duration at 7 years of age was maternally reported in the 7-year-questionnaire on the item ‘Approximately how many hours does the child usually sleep on a weeknight?’ with the following response categories: (1) 8 h or less, (2) 9 h, (3) 10 h, (4) 11 h and (5) 12 h or more. Depressive symptoms were measured using the 13-item Short Mood and Feelings Questionnaire⁵⁷ reported by mothers in the 8-year-questionnaire. We prepared the questionnaire-assessed phenotypes using the phenotools R package (0.2.8) (ref. ⁵⁸) in R 4.1.0⁵⁹. Educational achievement at age 10 was assessed as a composite score based on performance on three national standardized tests of skills in literacy, numeracy and English, derived from the Statistics Norway Educational Registry. The national tests are administered in the autumn of grade 5 (age 10), grade 8 (age 13) and grade 9 (age 14), and are mandatory with exemptions only on application on the grounds that the results will not be useful for assessing the child’s learning due to disability or lack of knowledge of the Norwegian language. Because the raw score ranges differed across subjects and test years, the raw scores were standardized within each test year and subject test. Before standardization, outliers defined as values more than 4 standard deviations (s.d.) from the mean were dropped. The four phenotype scores were standardized for use in the trio analyses to place them on a common scale and to facilitate interpretation and comparison of effect sizes.

Parental phenotypes were not used in the analyses presented in this study. However, similar outcomes in parents to the four exemplar offspring phenotypes were defined to provide supplementary context for interpreting the offspring genetic and phenotypic results. Parental height was self-reported in centimetres in response to the question, ‘How tall are you?’ Educational achievement was derived from Statistics Norway data corresponding to the International Standard Classification of Education (ISCED)⁶⁰ definitions. The ISCED levels from 0 to 8 corresponding to ‘early childhood education (“less than primary”)’, ‘Primary education’, ‘Lower secondary education’, ‘Upper secondary education’, ‘Post-secondary non-tertiary education’, ‘Short-cycle tertiary education’, ‘Bachelor’s or equivalent level’, ‘Master’s or equivalent level’ and ‘Doctoral or equivalent level’ were coded as 1, 7, 10, 11, 13, 14, 16, 18 and 21 years of education, respectively. Sleep problems were self-reported by mothers when the child was 14 years old and by fathers in 2015–2016, using the total score of three items modified from the Karolinska Sleep Questionnaire⁶¹: ‘How often do you find it difficult to get to sleep at night?’, ‘How often have you woken up repeatedly during the night?’ and ‘How often do you feel tired or sleepy during the day?’ The response options ‘Never’, ‘Less than once a week’, ‘Once per week’, ‘Twice per week’, ‘Three times per week’ and ‘Four times or more per week’, were coded from 0 to 5, respectively. Parental depression symptoms were self-reported by mothers at about week 30 of pregnancy and fathers around week 15 of pregnancy, using the depression subscale of the short eight-item Hopkins Symptoms Checklist (SCL-8)⁶².

Information from the Medical Birth Registry of Norway (MBRN)⁶³, a national health registry containing information about all births in Norway from 1967 onwards, and the MoBa questionnaire data were used to identify registered sex, year of birth, multiple births (in the offspring generation) and reported parent–offspring relationships. Before genotyping quality control, all pedigrees were constructed based on the reported parent–offspring relationships for each pregnancy. Each family contained all reported relationships (these included parent–offspring, full-sibling relationships from pregnancies with single and multiple births, and half-sibling relationships). Wherever possible, registered sex was assigned using information from the MBRN. In instances where sex was not registered in the MBRN, the sex registered at birth reported in the MoBa questionnaires was used.

Biological data

Blood samples⁶⁴ were collected from participating mothers and fathers at approximately the 17th week of pregnancy during the ultrasound examination. A second blood sample was taken from the mother soon after birth. The blood sample for the child was taken from the umbilical cord after birth. Biological samples were sent to the Norwegian Institute of Public Health, where deoxyribonucleic acid (DNA) was extracted by standard methods and stored⁶⁵.

Genotyping and quality control

Genotyping of MoBa has been conducted through multiple research projects, spanning several years. Full details about MoBa genotyping are described in Supplementary Methods 1. The MoBaPsychGen pipeline for quality control and imputation was developed to ensure the complex relationship structure and varying selection criteria, genotyping batches and genotyping arrays were handled appropriately (pre-printed⁶⁶). Quality control was performed based on current best-practice protocols^{67,68,69,70,71}. Throughout the pipeline, both SNP and individual-level quality control were performed; with priority given to retaining individuals over SNPs, as SNPs can be imputed. The complete pipeline, with full details about quality control, phasing, imputation and post-imputation quality control, is described module by module in Supplementary Methods 2, and a pipeline overview is shown in Supplementary Fig. 1a.

The primary software used throughout the pipeline was PLINK v.1.90b7.2 and v.1.90b6.18 (ref. ⁷²) and the R Project for Statistical Computing v.4.0.5 (ref. ⁵⁹) was used to produce figures. Moreover, LiftOver tool^73,74, KING v.2.2.5 (ref. ⁷⁵), FlashPCA2.0 (ref. ⁷⁶), imputation preparation and checking script v.4.2.13 (ref. ⁷⁷), SHAPEIT2 release 904 (ref. ⁷⁸) with the duoHMM⁷⁹ algorithm, IMPUTE4.1.2_r300.3 (ref. ⁸⁰), QCTOOL v.2.0.8 and v.2.2.0 (ref. ⁸¹), ic tool v.1.0.8 (ref. ⁸²), PLINK2 v2.00a2.3LM (ref. ⁷²), IMPUTE2.3.2 (ref. ⁸³), yHaplo 2.1.12 (ref. ⁸⁴), the MitoImpute pipeline v.0.2 (refs. ^85,86), BCFtools v.1.8 and v.1.9 (ref. ⁸⁷), cat-bgen v.1.1.4, COSGAP v.1.0.0 containers (gwas.sif, python3.sif and r.sif from https://github.com/comorment/containers/releases/tag/v1.0.0; for details see the documentation at cosgap.readthedocs.io/en/latest/), Python 3.8.10 and perl v.5.32.1.

Trio genetic analyses

Exemplar trio analyses were conducted using the four offspring phenotypes. All analyses were carried out on complete case data and restricted to very well-imputed (imputation quality score ≥0.95) autosomal variants available after application of the quality control pipeline. All analyses adjust for offspring age and sex, and account for effects of technical covariates (genotyping batch and first 20 principal components; PC_1–20). Where analyses estimate or account for ‘indirect genetic effects’, these refer exclusively to variation in child outcomes associated with parents’ genotypes (that is, not those originating with siblings or other family members). The term ‘indirect genetic effects’ is used for consistency with the standards in the field, despite many methods being agnostic as to the specific mechanism, which can include indirect effects of parents, population structure (the presence of subgroups within a population that can bias genetic associations) and assortative mating (the tendency for individuals with similar traits to mate non-randomly, influencing genetic associations in the next generation) among other sources.

Trio-GWAS

We conducted within-family GWAS (Supplementary Note 4) using genotype data from parent–offspring trios, focusing on four offspring traits: educational achievement (n = 41,761), height (n = 23,786), sleep duration (n = 24,376) and depressive symptoms (n = 18,701). GWAS analyses involved fitting both population and within-family models to the same sample in linkage-disequilibrium-adjusted kinship v.6 (ref. ⁸⁸). To obtain population and within-family genetic association estimates, linear regression was first performed on offspring genotype without adjusting for parental genotype, and subsequently in a mutually adjusted analysis accounting for parental genotypes. The former is synonymous with a conventional principal-component-adjusted model. The latter extends this model by including the parental genotypes as additional covariates to estimate direct genetic effects, while accounting for (and estimating) indirect parental genetic effects. Standard errors were clustered at the family level to account for offspring relatedness. The models were specified as follows (omitting technical covariates for simplicity):

Population model:

Within-family model:

To quantify the attenuation of SNP associations in the within-family relative to population-based models, we compared estimates for genome-wide significant variants reported in previously published GWAS of height⁴⁴, educational attainment⁴⁶, sleep duration⁸⁹ and depression⁴⁸. For each phenotype, we extracted the set of independent genome-wide significant variants from the external GWAS and compared their association estimates in our MoBa population-based and within-family GWAS. Attenuation was calculated as the proportional decrease in SNP effect estimates between the population-based and within-family models, and 95% CIs were derived using a leave-one-out jackknife across the set of included variants. Furthermore, we illustrate the relationship between power and parental heterozygosity rate in Supplementary Fig. 2.

Linkage disequilibrium score regression

We used LDSC⁹⁰ (v.1.0.1 with default parameters) to estimate SNP-based heritability (\({h}_{\mathrm{SNP}}^{2}\)) and genetic correlations (r_G) between our population and within-family GWAS results and publicly available summary statistics from previous studies on height⁴⁴, BMI⁴⁵, educational attainment⁴⁶, smoking initiation⁴⁷ and depression⁴⁸. We input summary-level GWAS data into LDSC, using precomputed linkage disequilibrium scores from the 1000 Genomes Project European superpopulation. For the LDSC n parameter, we used the effective sample size (based on standard errors) for each phenotype and model. Effective sample size was computed per SNP using the following formula:

\(\text{Effective}\,n=\left(\frac{1}{{{\rm{s}}.{\rm{e}}.}^{2}}\right)\times \left(\frac{{{{\rm{s}}.{\rm{d}}.}_{{\rm{Resid}}}}^{2}}{(2\times {\rm{MAF}}\times (1-{\rm{MAF}}))}\right)\)

where s.e. is the standard error of the SNP effect estimates, MAF is the minor allele frequency, and s.d._Resid is the residual s.d. of the phenotype after covariate adjustment. SNP-based heritability and genetic correlations were calculated separately for the population-based and within-family GWAS models to assess the impact of adjusting for parental genotypes. We examined the attenuation of genetic correlations in within-family analyses compared with population-based models.

Polygenic-score-based approaches

Polygenic scores (PGSs) were calculated using LDpred2⁹¹, implemented in the bigsnpr package v.1.12.21 (ref. ⁹²), using the ‘LDPred2-auto’ option and, following established quality control procedures⁹³, based on quality control passing, well-imputed variants in MoBa that were present in an extended set of HAPMAP3+ (ref. ⁹⁴) variants. These restrictions resulted in a list of 1,310,867 SNPs. We used a precomputed linkage disequilibrium matrix from UK Biobank as the reference linkage disequilibrium panel⁹⁴. PGS were calculated based on the same external summary statistics used in the LDSC (height⁴⁴, BMI⁴⁵, educational attainment⁴⁶ and smoking initiation⁴⁷). All PGS were adjusted for the first 20 PCs and genotyping batch.

Trio-PGS

We performed trio-PGS analyses in a multiple regression framework, incorporating mothers’, fathers’ and children’s PGS for a given trait in a single model to estimate direct and indirect genetic effects on children’s height (n = 21,008), educational achievement (n = 40,101), sleep duration (n = 21,451) and depressive symptoms (n = 18,212). The models were specified as follows:

Estimates were calculated with robust standard errors with clustering on maternal IDs to account for the presence of siblings in the data. In sensitivity analyses to assess the sufficiency of this clustering strategy, we restricted to unrelated trios, removing at random all but one trio from families in which children were siblings, half-siblings, or first cousins – leaving the following number of trios per outcome: 14,541 for height, 27,926 for educational achievement, 15,026 for sleep duration, and 12,632 for depressive symptoms. To estimate effects without adjustment for parental indirect genetic effects, we re-ran models without parental PGS.

Polygenic transmission disequilibrium testing

PTDT analyses were performed in sub-samples of MoBa families identified by selecting genotyped children with extreme values for height (>2 s.d., above the sample mean; n = 527) educational achievement (>90th percentile; n = 3,227), sleep duration (below 10th percentile; n = 438), and depressive symptoms (>2 s.d. above the sample mean; n = 618), where parental average values on equivalent measures were not similarly extreme. In PTDT, the genetic relationship between different traits or conditions is assessed based on the extent of over- or under-inheritance of genetic propensities by individuals selected based on a trait or condition³⁹. This is calculated as the average deviation of selected children’s PGS from the mid-parental PGS average, scaled to the standard deviation of the mid-parental PGS. Siblings (n = 84 for height, 535 for educational achievement, 78 for sleep duration, and 95 for depressive symptoms) were used as negative controls in these analyses. In our PTDT analyses, we assessed the over- or under-inheritance of genetic variants associated with height⁴⁴, BMI⁴⁵, educational attainment⁴⁶ and smoking initiation⁴⁷.

Trio-genome-wide complex trait analysis

Trio-genome-wide complex trait analysis (trio-GCTA)³⁶ was run to estimate the additive genetic variance attributable to mothers, fathers and children—across all measured SNPs—on measures of height, educational achievement, sleep duration and depressive symptoms. Genomic-related matrices were estimated for complete trios of genetically inferred European ancestry in which children had phenotypic data for at least one of the four phenotypes. In the case of more than one sibling having phenotypic data available, only one sibling was included at random. To remove individuals related across family trios a bottom-up algorithm from the OpenMendel package⁹⁵ was used with a relatedness threshold of 0.1. The final number of trios, by outcome, was as follows: 12,332 for height, 23,221 for educational achievement, 12,799 for sleep duration and 10,794 for depressive symptoms. Four alternative models were run: (1) estimating the variance explained by the effects of all trio members and the covariance between them (full model); (2) dropping the covariance terms from the full model (no covariance model); (3) estimating the additive genetic effects of the child only (direct effects only model); and (4) estimating no genetic effects (null model). All models included child’s age (in years) at collection of the phenotype, sex, genotyping batch and parents’ first 20 PCs as fixed effect covariates. Akaike’s information criterion was primarily used to select the best model, whereas results of likelihood ratio tests (LRTs) was used to evaluate loss of fit when model parameters were removed. LRT and Bayesian information criterion were consulted for competing models with similar Akaike’s information criterion values. For the likelihood ratio test, we calculated P-values using the classical procedure with a single chi-square null distribution when covariance component parameters were removed. When variance component parameters were removed, we calculated P-values using the equation below, where \(\hat{\lambda }\) denotes the test statistic, and j is the number of parameters removed so that the null distribution is a mixture of chi-square distributions dependent on j (refs. ^96,97):

LDSC uses GWAS summary data, meaning its power is influenced by both the heritability of the trait and the precision of the underlying GWAS. By contrast, trio-GCTA does not use any summary data and estimates heritability as part of a variance decomposition, using individual-level genetic data. Therefore, although both methods can be used to estimate SNP heritability, their properties (and the underlying data they use) differ substantially.

Assortative mating

Whereas trio estimates of direct genetic effects tend to be robust to assortative mating and population structure, estimates of indirect genetic effects can reflect not only genetic nurture effects but also assortative mating and population structure. To aid interpretation of the potential impact of assortative mating on the indirect genetic effects estimates, we calculated spousal Pearson correlations with 95% CIs for the four PGS used in the trio-PGS and PTDT analyses.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.