Pupil size and pupillary light reflex in early infancy: heritability and link to genetic liability to schizophrenia

Introduction

The pupillary light reflex (PLR) is a physiological response to retinal or attended stimulus luminance serving to regulate the amount of light reaching the retina. Individual variation in the PLR (the relative amplitude the pupil constricts or the latency to initiate the constriction response) is thought to reflect altered sensory reactivity, and the measure has a wide range of clinical applications, including detection of parasympathetic dysfunction (Wang et al., 2016) and cholinergic deficits (Hall & Chilcott, 2018). Although an involuntary reaction, it has been shown to be modulated by top-down cognitive processes (Ebitz & Moore, 2017; Mathôt, Linden, Grainger, & Vitu, 2015). While the sudden constriction to light seen in the PLR is linked primarily to the parasympathetic nervous system and draws heavily on cholinergic synaptic transmission (Hall & Chilcott, 2018), fluctuations in the pupil more generally reflect also other brain systems, including the sympathetic activity and the locus coeruleus–norepinephrine (LC-NE) system more specifically (Aston-Jones & Cohen, 2005; Gilzenrat, Nieuwenhuis, Jepma, & Cohen, 2010; Murphy, Robertson, Balsters, & O'connell, 2011), which index ongoing psychological processes related to arousal and attention (Laeng, Sirois, & Gredebäck, 2012). Tonic levels (baseline levels) of pupil dilation thus both reflect the amount of light reaching the eye, and the level of attention and arousal of the individual in that context.

In recent years, there has been an increasing interest in the link between the PLR early in life and autism spectrum disorder (ASD). In children, adolescents, and adults, a hypo-sensitive response (slower latency and reduced relative amplitude) has been shown at a group level (Daluwatte et al., 2013; Fan, Miles, Takahashi, & Yao, 2009) and associated with ASD traits (DiCriscio & Troiani, 2017). In contrast, early in life, a hypersensitive PLR has been seen in infants with a later ASD diagnosis and associated with symptom severity (Nyström et al., 2018; Nyström, Gredebäck, Bölte, & Falck-Ytter, 2015). Atypical PLR responses have also been described in psychiatric conditions such as schizophrenia, where a hypo-sensitive response has been found at a group level (Hakarem, Lidsky, & Sutton, 1966) and in association with more negative symptoms (i.e. symptoms related to withdrawal, lack of motivation and interest were associated with a slower latency to constrict; Fattal et al., 2020).

Twin studies help us understand the nature of individual differences in terms of genetic and environmental contributions to variability in traits. So far, no study that we are aware of investigated pupil measures such as the PLR in infants in a genetically informed design. This is striking given the link between the PLR and strongly heritable neurodevelopmental conditions like ASD and schizophrenia (Daluwatte et al., 2013; Fan et al., 2009; Fattal et al., 2020). The current infant twin study aimed to clarify the etiological background and genetic links of the PLR (the relative constriction amplitude and the latency) and the pupil size (baseline pupil diameter) at 5 months of age. The primary aim was to quantify the contribution of genetic and environmental factors to the PLR and baseline pupil size, and to the covariance between the different pupil measures. In light of previous research indicating significant genetic effects on early brain development (e.g. morphology measures; Maggioni, Squarcina, Dusi, Diwadkar, & Brambilla, 2020), we expected all pupil measures to be heritable. Our second aim was to utilize genome-wide polygenic scores (GPSs, a genotype-based individual score that summarizes the estimated effect of many trait-associated alleles) to assess whether the common genetic architecture known to be associated with psychiatric and neurodevelopmental conditions, particularly ASD and schizophrenia, associates with infant pupil parameters.

Methods and materials

Participants were part of the BabyTwins Study Sweden (BATSS), a longitudinal infant twin study (Falck-Ytter et al., 2021). BATSS has a classic twin design, in which similarity is compared between monozygotic twin pairs (MZ; who share 100% of their segregating alleles) and dizygotic same-sex pairs (DZ; who on average share 50% of their segregating alleles). Twins were identified via the Swedish population registry and contacted by letter. Only the greater Stockholm area was targeted due to the need to come in person to the lab. Inclusion and exclusion criteria, checked via telephone interview, were opposite-sex twin pairs, diagnosis of epilepsy, known presence of a genetic syndrome, vision or hearing impairments, very premature birth (prior to week 34), presence of developmental or medical conditions likely to affect brain development (e.g. Cerebral Palsy), and infants where none of the biological parents were involved in the infant’s care. Parents also needed to be willing to disclose information about family medical and psychiatric history, demographics, and birth delivery. From 2016 to 2020, 311 families (29% of the entire population of same-sex twins born in the area; 57% MZ pairs, as estimated based on DNA sampled from all infants) took part in an in-person multi-methods assessment at 5 months (testing site at Karolinska Institutet). The study sample has a high socioeconomic status, and it comprises of mostly Swedish-born families (90% of twin pairs had at least one parent born in Sweden). Recruitment summary and sample demographics are fully reported in the BATSS description paper (Falck-Ytter et al., 2021). Parents gave informed consent to take part. BATSS was approved by the regional ethics board in Stockholm and was conducted in accordance with the Declaration of Helsinki.

For the current analyses, 28 twins were excluded due to parent-reported twin-to-twin transfusion syndrome (12 pairs), report of seizures at the time of birth (2 individuals from different pairs), very low birth weight (<1.5 kg, 1 individual), and report of spina bifida (1 individual). In addition, 23 infants did not take part in the eye-tracking assessment due to technical reasons (3 pairs), lack of time (1 pair plus 1 individual), bad calibration (1 individual), and tiredness (4 pairs plus 5 individuals). From the infants that took part in the assessment, 61 infants (9 pairs plus 43 individuals) did not have enough valid trials to yield pupil measures (see below criteria for measurements), so these were excluded from the analysis. The final sample was of 510 individuals (280 incomplete pairs, i.e. pairs with at least one twin with data). We found no statistically significant differences between the excluded and included infants in terms of either parental education level, family income, sex, age or any of the GPSs used in this study.

Measures

The PLR was measured during a lab visit at 5 months of age which included other paradigms not related to the current study (Falck-Ytter et al., 2021). During this visit, both twins performed the PLR task; one twin was assessed at a time.

Eye tracking protocol

To record gaze and pupil data a Tobii TX300 eye-tracker was used (sampling rate of 120 Hz) with scripts written for the Eurosibs study (Jones et al., 2019) in MATLAB (version R2013b, MathWorks, Natick, MA, USA) and Psychtoolbox (version 3.0.12). An initial 5-point calibration was called before the start of the task battery, which lasted for about 10 min and involved rotations of free-viewing of dynamic and static scenes (mixture of social and nonsocial content), gaze-contingent gap-overlap trials, PLR measurements, and postcalibration sequences. Here, we analyzed only pupil data collected during the presentation of 12 trials that elicited a PLR. These 12 PLR measurements were not collected in sequence, rather they were always interspersed with other stimuli with varying luminosity.

The pupillary light reflex task

For each trial, lasting approximately 6 s, the stimulus consisted of a small central fixation point on a black background (0.9 lux) that flashed white (190 lux) for 75 ms with a random onset between 1,600 and 2,400 ms. In line with previous research, the PLR was evaluated both in terms of its relative constriction and its latency. The relative pupil constriction was calculated as in Fan et al. (2009) by the formula (A02–Am2)/A02, where A0 is the average pupil diameter before the onset of the PLR (during an interval starting 100 ms before and ending at the PLR onset, as determined by the PLR latency), that is, the mean baseline, and Am is the minimum pupil diameter in the interval 500–1,500 ms relative to the stimulus onset (Fan et al., 2009). Following previous work (Nyström et al., 2015, 2018), the PLR latency was defined at the acceleration minimum between flash onset and minimum pupil diameter. Computational processing was done in MATLAB using the TimeStudio framework (Nyström, Falck-Ytter, & Gredebäck, 2016) and is described elsewhere (Nyström et al., 2015, 2018). All trials were inspected by a researcher (blind to zygosity). Invalid trials (with apparent artifacts and implausible responses) for each eye were first automatically marked, and then visually rejected. Only one eye per trial was kept and infants with less than three valid trials were excluded. Mean values from valid trials were computed to get the outcome variables (relative constriction; latency of the response; and baseline pupil diameter). Short-term internal stability of the pupil measures was estimated by intra-class correlations (two-way random effects) using mean values of the first six trials (test block) and of the last six trials (retest block) – relative constriction had moderate reliability (ICCtwin 1 = .60/ICCtwin 2 = .58), baseline pupil diameter had good reliability (ICC = .83/.87), whereas latency of the response had poor reliability (ICC = .12/.18).

Genome-wide polygenic scores

The DNA samples were genotyped using Infinium Global Screening Array (Illumina, San Diego, CA, USA). Genotype quality control and processing was done using standard procedures and described elsewhere (Falck-Ytter et al., 2021). GPSs were calculated based on the most recent and largest genome-wide association studies for a range of neurodevelopmental and psychiatric conditions [ADHD (Demontis et al., 2019), ASD (Grove et al., 2019), bipolar disorder (Stahl et al., 2019), major depressive disorder (Howard et al., 2019), and schizophrenia (Consortium, Ripke, Walters, & O’Donovan, 2020)]. The GPSs were calculated using the PRS-CS (polygenic prediction via Bayesian regression and continuous shrinkage priors) method (Ge, Chen, Ni, Feng, & Smoller, 2019). For the GPS analyses, the first 10 principal components of ancestry were included as covariates in the models. In the sample, only one pair was of completely non-European genetic ancestry. This pair was retained in the analysis since their GPSs for all five conditions abovementioned were below 3 standard deviations of the sample.

Data analysis

The analysis plan for this study was preregistered in OSF (https://osf.io/5kzh2), after data collection, cleaning, and extraction, but prior to any analysis of the data. PLR latency, PLR relative constriction, and pupil baseline had skewness values between −1 and 1 (.04 for Latency, .23 for Relative Constriction and .59 for Baseline) so they were not transformed. All measures were standardized. R software (version 4.0.0) was used for all analyses.

For twin models, the OpenMx package (version 2.18.1; Neale et al., 2016) with full-information maximum likelihood estimation was used. Even though these methods allow for partially complete pairs (one twin missing data) to be included, they were not included in the current analysis because the twins with missing data were also missing certain covariates (which have to be complete for all cases).

The effects of covariates on the PLR parameters were estimated within the twin model. These included the number of valid trials (a proxy of data quality), the year of data collection (a proxy of possible changes in luminance levels of the testing room, even though no systematic differences were detected), age (in days, same for all twins, except for six pairs where the average age was calculated), and sex. In a deviation from the preregistered plan (covariates would be added if found to be associated with the relevant outcome PLR measures), these variables were added as covariates in all twin analyses; the findings do not change between approaches.

Univariate twin models estimate the relative contribution of genetic and environmental factors to the variation in a phenotype. This variation can be partitioned into additive genetic influence (A, heritability, which increases twin similarity), nonshared environment (E, environmental influences that differ between twins and decrease twin similarity, including measurement error), and either shared environment (C, environmental influences that increase twin similarity regardless of zygosity, e.g. family socioeconomic status) or nonadditive genetic (D) effects (D and C variance cannot be estimated simultaneously with twin data alone, because they confound one another). These components can be estimated by comparing the resemblance in MZ and DZ pairs. When MZ pairs are more similar than DZ pairs we can infer that variance in a phenotype is influenced by genetic factors; when DZ pairs are more than half as similar as MZ pairs, we can infer that shared environment influences the variance in that phenotype.

Univariate twin models were fit separately for each pupil measure. First, to test for the assumptions of equality of mean and variances and to derive twin correlations, a constrained saturated model, in which the means and variances were constrained to be equal across twin order and zygosity, was fit. Next, using either a model with A, C, and E (ACE model, where A stands for additive genetic effects, C for shared environmental effects, and E for nonshared environmental effects) or a model with A, D, and E (ADE model, where A stands for additive genetic effects, D for nonadditive or dominance genetic effects, and E for nonshared environmental effects), we estimated and tested the genetic and environmental contributions to the variance in the measures. The best-fitting model was the most parsimonious (fewer number of parameters estimated) and nonsignificant model (meaning that there was no decrement in fit compared to the genetic model, indexed by the χ2 distribution).

Bivariate twin modeling allows the relative contributions of genetic and environmental factors to the covariation between two measures to be estimated. It assumes that the observed correlation between measures can be decomposed into genetic, shared environmental, and unique environmental effects. In the bivariate case, cross-trait cross-twin correlations (CTCTs) are compared, that is, the correlation between one measure for one twin and another measure for their co-twin, between MZ and DZ pairs.

Bivariate twin models were fitted between phenotypes that were phenotypically correlated and shown to be heritable in univariate analyses. To test for the assumption of the equality of phenotypic and CTCTs in the data and to derive these correlations, a constrained saturated model, in which the phenotypic and CTCTs were constrained to be equal across twin order and zygosity (for the phenotypic correlation), was fitted. Next, using a bivariate ACE or ADE model and a correlated factors solution, we examined genetic and environmental correlations and bivariate heritability (i.e. the proportion of the phenotypic correlation that is explained by genetic factors).

To investigate the associations between the pupil measures and GPSs, a generalized estimating equation (GEE) framework (using the drgee package; Zetterqvist & Sjölander, 2015), with cluster-robust SE to account for family relatedness, was used to derive linear regression Beta estimates and p-values. Pupil outcomes were regressed on age and sex before analysis. P-thresholds were adjusted for the number of outcomes. To estimate effect sizes (ΔR2) for model comparison, the R2 of corresponding models were derived and compared to the R2 of the null model (i.e. the model with only covariates included).

Associations with the Infant and Toddler Sensory Profile, which were planned in the preregistration, are reported in Appendix S1 (no statistically significant associations were found, Tables S1 and S2).

Results Twin modeling analyses

Sample descriptive statistics are presented in Table 1. In terms of the effects of the covariates, there were no significant covariate effects for latency to constrict (see Table S3 in Appendix S2). Age was positively related to relative constriction (β = .02, p < .001, see Table S4) and to baseline (β = .01, p = .025, see Table S5), that is, older infants had increased constriction amplitudes and bigger baseline pupils. Sex was related to baseline diameter (β = .37, p = .001, infant boys had, on average, .16 mm larger pupils, see Table S5). These patterns of age and sex effects are in line with previous literature with infants and toddlers (Kercher et al., 2020). Number of valid trials and year of testing were statistically associated with baseline diameter (more valid trials yielded a smaller pupil size, β = −.04, and pupil size was smaller in the last year of data collection, β = −.14; both p < .01, see Table S5).

Table 1. Descriptive statistics Overall MZ females MZ males DZ females DZ males N (twins) 510 129 152 116 113 Age (days)

167.70/8.59

[145–203]

168.01/9.16

[153–194]

167.41/8.09

[150–187]

168.03/8.24

[153–189]

167.41/8.96

[145–203]

Pupillary light reflex Number of valid trials

8.48/2.85

[3–12]

8.33/2.67

[3–12]

8.47/2.92

[3–12]

8.68/2.79

[3–12]

8.44/3.04

[3–12]

Mean latency (s)

0.35/0.04

[0.18–0.52]

0.35/0.05

[0.24–0.49]

0.36/0.04

[0.24–0.52]

0.35/0.05

[0.22–0.49]

0.35/0.04

[0.18–0.45]

Mean relative constriction (%)

28.78/7.75

[2.42–51.90]

28.40/7.63

[11.61–51.90]

29.03/7.67

[9.90–48.57]

28.54/7.98

[11.90–50.46]

29.11/7.85

[2.42–47.53]

Mean baseline (mm)

3.23/0.46

[2.27–4.96]

3.17/0.48

[2.31–4.96]

3.33/0.48

[2.27–4.85]

3.13/0.38

[2.27–4.19]

3.28/0.45

[2.39–4.69]

Statistics presented as mean/SD [min max]. DZ, dizygotic; MZ, monozygotic.

Table 2 shows the twin correlations for the pupil parameters. Higher MZ than DZ correlations suggests heritable influences on baseline and relative constriction. Low shared environmental influences were suggested because the DZ correlations were only marginally greater than half the MZ correlations. For latency, there was very low similarity within both MZ and DZ pairs, suggesting no familial effects on the phenotype. For this reason, further results for latency were not reported (but see Appendix S2).

Table 2. Twin correlations for the PLR measures, separate for MZ and DZ pairs. 95% confidence intervals are shown in parentheses N (Twin pairs) Latency Relative constriction Baseline MZ 127 .15 [−.03–.32] .61 [.49–.70] .64 [.54–.72] DZ 103 .28 [.10–.44] .36 [.19–.51] .34 [.12–.51]

Univariate twin modeling (see model assumptions tests in Tables S4–S6) confirmed the different patterns of genetic and environmental contributions to the pupil measures. Table 3 shows the twin models fit statistics and the standardized estimates.

Table 3. Univariate twin model fit statistics and parameter estimates. Best-fitting model in bold Model # Parameters −2LL df AIC Comparison model Δχ2 Δdf p Value A C E Relative constriction Fully Sat 14 1198.18 446 306.18 NA – – – – – – ACE 8 1200.17 452 296.17 Fully Sat. 1.99 6 .920 .49 .12 .39 AE 7 1200.64 453 294.64 ACE 0.47 1 .491 .62 0 .39 CE 7 1207.66 453 301.66 ACE 7.49 1 .006 0 .49 .51 E 6 1271.98 454 363.98 ACE 71.81 2 <.001 0 0 1 Baseline Fully Sat 14 1155.17 446 263.17 NA – – – – – – ACE 8 1163.45 452 259.45 Fully Sat. 8.28 6 .228 .61 .03 .36 AE 7 1163.47 453 257.47 ACE 0.02 1 .884 .64 0 .36 CE 7 1174.48 453 268.48 ACE 11.03 1 .001 0 .54 .46 E 6 1249.67 454 341.67 ACE 86.22 2 <.001 0 0 1 Model definitions. The Fully Sat. model is the fully saturated model of the observed data, which models the means and variances separately for each twin in a pair and across zygosity. −2LL = fit statistic, which is minus two times the log-likelihood of the data; df = degrees of freedom; AIC, an alternative fit index, lower values denote better model fits; Δχ2 = difference in −2LL statistic between two models, distributed χ2; Δdf = difference in degrees of freedom between two models.

For relative constriction, there were significant high genetic effects (best-fitting AE model A = .62, CI: .51–.70); while the contribution of the shared environment was nonsignificant (ACE model C = .12, CI: 0–.42). For baseline diameter, high genetic effects were significant (best-fitting AE model A = .64, CI: .54–.72), and the shared environment had a minimal and nonsignificant contribution (ACE model C = .03, CI: 0–.37).

The phenotypic correlation between baseline and relative constriction was moderate (rPh = .38, CI: .28–.46), suggesting that infants with a larger pupil at baseline had higher constriction amplitudes. A saturated model (see the bivariate model assumptions tests in Table S7) was used to derive the CTCT correlations, which suggested genetic influences on the phenotypic correlation because the MZ CTCT (rCTCT MZ = .24, CI: .13–.34) was higher than the DZ CTCT (rCTCT DZ = .06, CI: −.09–.21). Because the correlation in DZ was less than half of the correlation in MZ twins, some nonadditive genetic effects could explain the phenotypic correlation between the measures. However, given the sample size and the small (although nonsignificant) contribution of shared environment seen in the univariate models, we did not fit an ADE model. Table 4 shows the model fit statistics for the bivariate model. An AE model fitted the data best, so this model was used to derive separately the genetic correlations between pupil baseline and relative constriction (Figure 1).

Table 4. Bivariate twin model fit statistics. Best-fitting model in bold Model # Parameters −2LL df AIC Comparison model Δχ2 Δdf p Value rA rC rE Fully Sat 36 2270.68 884 502.68 NA NA NA NA – – – ACE 19 2289.04 901 487.04 Fully Sat. 18.36 17 .367 .64 −1 .36 AE 16 2292.15 904 484.15 ACE 3.11 3 .375 .38 .38 CE 16 2305.04 904 497.04 ACE 16.00 3 .001 – .34 .42 E 13 2450.21 907 636.21 ACE 161.17 6 <.001 – – .38 Model definitions. The Fully Sat. model is the fully saturated model of the observed data, which models the means and variances for both variables, and the phenotypic and cross-twin-cross-trait correlations between the two variables, separately for each twin in a pair and across zygosity. The best-fitting model (in bold) was the nonsignificant and most parsimonious model. −2LL = fit statistic, which is minus two times the log-likelihood of the data. df = degrees of freedom; AIC, an alternative fit index, lower values denote better model fits; Δχ2 = difference in −2LL statistic between two models, distributed χ2; Δdf = difference. image

The AE bivariate model was the best fit for baseline and relative constriction. This was used to derive the genetic correlations between phenotypes using a correlated factors solution

The proportion of the phe

留言 (0)

沒有登入
gif