Causal Relationships Between Screen Use, Reading, and Brain Development in Early Adolescents

Abstract The rise of new media has greatly changed the lifestyles, leading to increased time on these platforms and less time spent reading. This shift has particularly profound impacts on early adolescents, who are in a critical stage of brain development. Previous studies have found associations between screen use and mental health, but it remains unclear whether screen use is the direct cause of the outcomes. Here, the Adolescent Brain Cognitive Development (ABCD) dataset is utlized to examine the causal relationships between screen use and brain development. The results revealed adverse causal effects of screen use on language ability and specific behaviors in early adolescents, while reading has positive causal effects on their language ability and brain volume in the frontal and temporal regions. Interestingly, increased screen use is identified as a result, rather than a cause, of certain behaviors such as rule‐breaking and aggressive behaviors. Furthermore, the analysis uncovered an indirect influence of screen use, mediated by changes in reading habits, on brain development. These findings provide new evidence for the causal influences of screen use on brain development and highlight the importance of monitoring media use and related habit change in children.

brain development in children and adolescents.The dataset utilized in this study included brain imaging data, questionnaires, and genetic data collected during both the baseline and two-year follow-up assessments.The ABCD study recruited 11,875 children aged 9 and 13 from 21 research sites across the United States between October 2016 and October 2018.The same cohort was reassessed using similar measurements two years later.All study procedures were approved by institutional review boards at individual sites, and written consent was obtained from parents with verbal assent from children.Comprehensive details about the ABCD study such as data types, preprocessing, and ethical considerations can be found elsewhere (Developmental Cognitive Neuroscience Special Issue 2018, vol.32, pp.1-164).We excluded preterm-born infants with gestational age less than 37 weeks (n = 1629) from our analyses.For siblings, we randomly excluded one participant from each pair with a closed gene relationship (pi-hat > 0.4) estimated with individual genotype data using PLINK (v1.9).Finally, 7107 subjects remained at baseline and 4505 at two-year follow-up (see Table 1 for detailed sociodemographic information).Note that the reduction in subjects at the 2-year follow-up is mainly because ABCD Release 4 did not fully open the second year's data.Meanwhile, the removed siblings were identical at both time points.

Exposure variables
The study focused on analyzing six different exposure variables related to screen use and reading time.Screen use time was obtained from the ABCD Youth Screen Time Survey [4,5] .This survey involved self-reporting by the children regarding their daily time spent on various electronic activities including watching TV shows or movies (TVM), watching videos (Video), playing video games (Game), texting on a device (Text), visiting social networking sites (SNS), and video chatting (VC) on both weekdays and weekends.A weighted mean time was calculated using the formula [(5*Tweekday + 2*Tweekend)/7] to represent the daily time for each screen type.At the 2-year follow-up, the survey was slightly modified, and the time spent playing video games was split into single-player and multiplayer [6] .Therefore, we summed these two times together to get the daily time for video games.
Reading time was assessed using the ABCD Parent Sports and Activities Involvement Questionnaire (SAIQ), which was a parent-reporting questionnaire for information about the involvement in sports, music, reading, and other hobbies for the children [5] .The weekly reading time reported by the parents was divided by 7 to calculate the average reading time per day.For the 2-year follow-up, the longitudinal version of the SAIQ was used.Because children's self-reported reading time was not collected in the first two years of the ABCD study, we used the parent-reporting questionnaire.Nevertheless, all the assessments in the ABCD study were well-designed and proven to be reliable and validated by other studies [7,8] .To further ensure the reliability of our results, we used the data from two time points for crossvalidation, and only the consistently significant results would enter the subsequent analysis, which may reduce the bias due to the inaccuracies in the questionnaires.

Non-neuroimage outcome variables
Non-neuroimaging variables in this study included children's cognitive ability and behavior.
Cognitive ability was assessed using five age-corrected standard scores obtained from the NIH Toolbox [5,9] , including picture vocabulary, flanker, pattern comparison processing speed, picture sequence memory, and oral reading recognition tests.Two other tests (list sorting working memory test and dimensional change card sort) were not available during the second assessment and were therefore excluded from the analysis.Behavior was measured according to the ABCD Parent-Child Behavior Checklist Scores Aseba (CBCL), which assesses 8 behavioral dimensions (t-scores), including anxious/depressed symptoms, withdrawn/depressed symptoms, somatic complaints, social problems, thought problems, attention problems, rulebreaking behavior, and aggressive behavior [5,10] .

Neuroimaging variables
The acquisition and preprocessing of image data have been described in detail in other articles [11,12] .Briefly, MRI scans were acquired using 3T scanners at all participating ABCD sites, with optimized and standardized acquisition protocols.A unified pipeline [11] was then employed for data preprocessing.The preprocessing of T1-weighted and T2-weighted images involved the correction of gradient nonlinearity and intensity inhomogeneity, rigid registration to a reference brain in standard space, and segmentation of brain tissues.Quality control measures were implemented both before and after the image processing to ensure the exclusion of poor-quality images.This study focused specifically on brain volume estimated using the Destrieux atlas (148 ROIs in the cortex [13] ) in combination with the subcortical regions in the ASEG atlas [14] (12 subcortical regions including thalamus, caudate, putamen, pallidum, hippocampus, and amygdala in both hemispheres).

Confounding variables
To reduce unbiased results, we adjusted for common confounding variables in our study, including age, sex, household income, parental education, race, ethnicity, prematurity, sibling, and site effect.Race was categorized into six groups: white, black, Native American, Pacific Islander, Asian, and other.For the sake of brevity, Table 1 shows only the household income of the five levels (< $11999, $12000 -$24999, $25000 -$49999, $50000 -$99999, > $100000 annually) and the proportion of parents with a college education.In addition, because there were multiple MRI machines at the same site, we used the serial number of the MRI scanner (to take into account differences MRI scanners at the same site) as a covariate in the analysis of brain volume.For simplicity, we also referred to this as the "site effect".

Genotype data and processing
In our study, we used the public release 3.0 genotype data in the ABCD (v4) datasets.The dataset included both pre-and post-imputed genotype data.We referred to the previous methods in the related paper [15] to estimate population stratification, and combined the preimputed genotype data with the 1000 Genome Phase 3 data (https://www.internationalgenome.org/).We used the first five principal components to identify population clusters using UMAP, resulting in 7 broad populations including Africans, Americans, Bengali, Finnish Europeans, Non-Finnish Europeans, East Asians, and South Asians.
Only non-Finnish Europeans (n = 6161) were included in the subsequent gene-based analyses in the ABCD data.
We performed the GWAS for phenotypic variables (e.g., exposure and outcome variables) using the fastGWA method [16] implemented in Genome-wide Complex Trait Analysis (GCTA [17] ).We included age, sex, and the top 20 principal components as covariates and accounted for relatedness using a sparse genetic relatedness matrix (GRM) with a cutoff of 0.05.The GWAS results of screen use and reading will be used for subsequent Mendelian randomization (MR) analysis.
In this study, we used genetic data in a public dataset to obtain potential causal relationships between variables.However, the use of genetic data also carries certain risks, including invasion of privacy, genomic discrimination, and other ethical and legal issues [18] .Therefore, before using genetic data, it is best to fully understand the relevant policies and regulations and comply with the data use norms.

Statistical analyses
To examine the relationship between exposure variables (e.g., reading and screen use) and cognitive skills, behavior, and brain volume.We first used a linear mixed-effects (LME) model to explore the association between these variables.To infer a causal relationship between the variables, we performed a one-sample MR analysis based on the LME analysis results.
Considering the small sample size of ABCD for GWAS analysis, the results of one-sample MR analysis might be affected by sample selection bias [19] .We attempted to use related GWAS summary data from larger sample cohorts to perform two-sample MR analysis, to validate the results of one-sample MR.It is important to note that for some of the variables we used, GWAS summary data were not available, and the existing data were based on adults rather than adolescents.To examine whether screen use indirectly influences the outcomes by changing reading habits, we first conducted a mediation analysis with reading time as the mediator and then performed a two-step MR analysis.In all statistical analyses, we only included participants with complete data to ensure data quality.

Association analysis
LME models were used to examine the association between each exposure variable and each outcome variable (e.g. 5 cognitive scores, 8 behavioral scores, and brain volume in 160 brain regions).Each exposure and outcome variable pair was analyzed separately, while fixed effects for age, sex, race, ethnicity, household income, parental education, and random effects for the site were included in all models.The fitlme function in MATLAB (v2018a) was used for model fitting.
~ (1|site) +  +  + ℎ +  +  +  + ɛ Association analysis was performed separately at two time points.Bonferroni correction for multiple comparisons was applied at each time point for all non-neuroimage outcome variables, and FDR correction was applied for neuroimage variables, because conservative Bonferroni correction for hundreds of brain regions may be too strict considering that we only included results that were significant at both time points.

Genetic correlation
We estimated the genetic correlation among seven exposure variables and their genetic correlation to the outcomes.The genetic correlations were estimated by the LDSC [1] method implemented in the Genomic SEM [20] package in the R software.

One Sample Mendelian randomization analysis
To examine the causal relationship between exposure and outcome variables, we used the one-sample MR method [21] in this study.This method involves two-stage least-squares (TSLS) regression with SNPs as instrumental variables (IVs) [22] , we used the OneSampleMR R package to perform the analysis and used a liberal threshold of p < 1e -5 and MAF > 0.005 in selecting the instruments for the exposure variables to maximize the power of the genetic instruments [22] .Linkage disequilibrium (LD) pruning with an R 2 threshold of 0.001 and a window size of 1 Mb to remove SNPs in high LD.SNPs that showed significant correlations with the outcome variable (using the same p-value threshold as for the exposure variable) were removed.Each SNP was coded as 0, 1, or 2, depending on the number of effect alleles carried by a subject.
We included covariates such as age, sex, household income, parental education, and site in both regressions.The ABCD baseline data were used for the analysis because the genetic data is only available at baseline.
The MR analysis was conducted among the significant results in the above association analysis, and the FDR correction was used for multiple comparisons in all non-neuroimage (n = 33) and neuroimage (n = 119) outcome variables, respectively.Furthermore, the reverse MR analysis was also conducted to investigate the causal influence of those outcomes on the children's screen use and reading habits.The reverse MR analysis was similar to the above MR method except that we exchanged the exposure and outcome variables [23] .
The sensitivity analysis includes weak instruments and the Sargan test.In the one-sample MR analysis, if the SNPs-exposure associations were weak, chance variation means that genetic association with the exposure and outcome are correlated in the direction of the confounded association between the two, resulting in increased false positive rates [24] .Weak instrument test is performed using the F statistic for one-sample MR and the F should be greater than 10 [25] .In addition, the Sargan test is applicable if the number of instrumental variables exceeds the number of endogenous variables, and a positive test result means that at least one instrumental variable is not valid.Therefore, significant results were needed to meet three criteria: original p-value survived after FDR correction, F > 10 in the weak instruments test, and not significant in the Sargan test.

Two-sample Mendelian randomization analysis
For screen use, we did not find GWAS summary data for the above six types of screen use, so we used GWAS data on leisure screen time using UK Biobank data (n up to 526725 [26] ).For cognitive ability, we used GWAS data on reading ability (n = 27180 [27] ) and general cognitive performance (n = 257841, [28] ).We selected the instrumental SNPs using the PLINK software and took the genomic data of the European superpopulation in the 1000 Genomes Project as the LD reference, which was restricted to bi-allelic SNPs with minor allele frequency > 0.1.We set the LD pruning r 2 = 0.001, the window size = 10000 kilobase pairs, and the P value associated with the exposures of 5e -8 .We then removed the SNPs that showed significant association with outcomes (p < 5e -8 ) and harmonized the data using the harmonise_data function in the TwoSmpleMR package to ensure that the same allele was used to estimate the genetic variant association.Inverse variance-weighted (IVW) regression [29] with multiplicative random effects was conducted as the primary analysis method.For the significant results with the IVW method, we would conduct four other methods to assess the robustness of the results including weighted median method, weighted mode, MR-Egger, and MR-PRESSO.The Wald ratio method would be used when only one instrument was available.All of these methods were conducted in the TwoSmpleMR package.In addition, we performed the MR-Egger regression to examine the potential bias of directional pleiotropy [30] , the MR-PRESSO method to the horizontal pleiotropy [31] , and a leave-one-out analysis to check whether the causal association was driven by a single SNP.There was no overlap between the data used in the one-sample and two-sample MR analyses.
Both one-and two-sample MR use SNPs as instrumental variables to estimate the causal effects of the exposure on the outcome, and both methods need to satisfy three assumptions, namely relevance assumption, independence assumption, and exclusion restriction.However, one-sample MR uses individual-level data and is usually estimated by the TSLS method, while two-sample MR generally uses GWAS summary data from two samples and is usually estimated by the IVW method [23,25] .Given the difference in these basic settings, there are both advantages and disadvantages for the two MR methods.On the one hand, the same population in the one-sample MR could guarantee the homogeneity in the association of the genetic variants with the exposure and outcomes, while the different samples in the twosample MR may affect the validity of the instrumental variable assumptions because a genetic variant may not be as strongly associated with the exposure in the outcome dataset.
Meantime, the difference in population characteristics such as age, sex, and socio-economic background might affect the interpretation of causal estimates in the two-sample MR analysis.
One-sample MR does not suffer from these problems, and it is even possible to compare the results of traditional statistical and MR analyses in the same study [24] .On the other hand, the overlap of participants will introduce sample selection biases and inflated false positive rate (type 1 error) in the one-sample MR analysis [19,32] , while two-sample MR without sample overlap could avoid such bias [24] .Therefore, we combined the analyses with the two methods to enhance the reliability of the conclusion in the present study.

Mediation analysis
To estimate the mediating effect of reading habits on the relationship between screen use and outcome variables, we conducted a mediation analysis using the mediation toolbox (https://github.com/canlab/MediationToolbox) in MATLAB v2018a.Specifically, the six types of screen use were the exposure variables, reading time was the mediator variable, those that showed a significant causal relationship with reading in the above analyses were taken as the outcome variables in the mediation analysis, and the covariates including age, gender, ethnicity, race, household income, parental education, and site were also considered in the analysis.Mediation analysis was conducted at baseline and two-year follow-up, respectively.

Two-step Mendelian randomization analysis
Two methods could be used to estimate the mediator effect in an MR framework, namely two-step MR analysis and multivariable MR analysis.The two-step MR analysis estimates the effect of the exposure on the outcome and compares this to the product of the effect of the exposure on the mediator multiplied by the effect of the mediator on the outcome.The multivariable MR analysis compares the standard model containing only the exposure variable and the model including the mediator and the exposure, to reflect the direct and indirect effect of the exposure on the outcome [24] .
To validate the above mediation analysis and generate a potential causal reference about the displacement relationship between screen use and reading habits in this study, we conducted a two-step MR analysis [23] to estimate the mediated effect of reading habits on the relationship between screen use and outcome variables.In the first step, we used MR analysis to reveal the causal effect of screen use on reading habits.In the second step, we used MR analysis to estimate the causal effect of the mediator (i.e., reading) on the outcome variables, which are significantly influenced by the mediator from the above analysis.We then used the "product of coefficients" method [33] to estimate the indirect effect and the "delta" method [34] to estimate the confidence interval [22] .Particularly, the indirect coefficient was equal to the product of two estimations (a * b, where a and b were the effect estimations in the above two steps, respectively), and the standard error in calculating the confidence interval was estimated as follows:

2 Figure S1 .
Figure S1.The Pearson correlation among all exposure variables at baseline (left) and 2-year follow-up (right).Abbreviation: TVM = watching television shows/ movies; Video: watching videos; Game = playing video games; Text = texting on a cell phone, tablet, and so on; SNS = visiting social networking sites; VC = video chatting.

Figure S2 .
Figure S2.MR estimation of the association between screen use and reading ability.Estimates showed the effect size (95% Cis) for the change in reading ability due to leisure screen use.The right panel shows the leave-one-out plot of the instrumental SNPs.

Figure S3 .
Figure S3.MR estimation of the association between screen use and cognitive performance.Estimates showed the effect size (95% Cis) for the change in cognitive performance due to leisure screen use.The right panel shows the leave-one-out plot of the instrumental SNPs.

Figure S4 .
Figure S4.Reverser MR estimation of the association between screen use and cognitive performance.Estimates showed the effect size (95% Cis) for the change of screen use time due to cognitive performance.The right panel shows the leave-one-out plot of the instrumental SNPs.