Development of thalamus mediates paternal age effect on offspring reading: A preliminary investigation

Abstract The importance of (inherited) genetic impact in reading development is well established. De novo mutation is another important contributor that is recently gathering interest as a major liability of neurodevelopmental disorders, but has been neglected in reading research to date. Paternal age at childbirth (PatAGE) is known as the most prominent risk factor for de novo mutation, which has been repeatedly shown by molecular genetic studies. As one of the first efforts, we performed a preliminary investigation of the relationship between PatAGE, offspring's reading, and brain structure in a longitudinal neuroimaging study following 51 children from kindergarten through third grade. The results showed that greater PatAGE was significantly associated with worse reading, explaining an additional 9.5% of the variance after controlling for a number of confounds—including familial factors and cognitive‐linguistic reading precursors. Moreover, this effect was mediated by volumetric maturation of the left posterior thalamus from ages 5 to 8. Complementary analyses indicated the PatAGE‐related thalamic region was most likely located in the pulvinar nuclei and related to the dorsal attention network by using brain atlases, public datasets, and offspring's diffusion imaging data. Altogether, these findings provide novel insights into neurocognitive mechanisms underlying the PatAGE effect on reading acquisition during its earliest phase and suggest promising areas of future research.


| INTRODUCTION
There has been a global trend of postponed childbearing, especially in developed countries (Kohler, Billari, & Ortega, 2002). This so-called "postponement transition" is primarily owing to changing patterns of education, employment, and marriage (Khandwala, Zhang, Lu, & Eisenberg, 2017;Sobotka, 2010). Although the research field is still in its infancy, increasing evidence reveals that advanced paternal age at childbirth (PatAGE) increases risks for a wide range of neuropsychiatric conditions, such as schizophrenia and autism spectrum disorders (D'Onofrio et al., 2014).
In contrast to mental health, few studies investigated the effects of PatAGE on offspring's cognition, such as reading, which is essential to success in modern society. The pioneering study in 1978 reported a negatively skewed distribution of PatAGE in 48 boys with reading disorders (RD; a.k.a. and referred here to dyslexia) (Jayasekara & Street, 1978). Four decades later, this topic remains underinvestigated and the existing findings are controversial. In a broader sample of 7-year-old children, Saha and colleagues revealed a significantly negative effect of PatAGE on offspring's reading after controlling maternal age at childbirth (MatAGE), gestational age, sex, and race (Saha et al., 2009). However, when parental education and number of siblings were added to the statistical model, the effect of PatAGE on reading was no longer significant (Edwards & Roff, 2010). Such inconsistency underlines the importance of more research that controls for possible confounds in examining the PatAGE effect on reading.
In addition to the controversy over the PatAGE effect on reading, no studies have yet examined the underlying mechanisms. Nascent research in molecular genetics, however, show that PatAGE explains nearly all the variance in the amount of de novo mutation, which is an alteration in a gene as the result of a mutation in a germ cell (egg or sperm) that increases by cell divisions of the gametes (approximately 38-fold in males at the age of 50 compared to females ;Breuss et al., 2020;J onsson et al., 2017;Kong et al., 2012). Hence, de novo mutation is the most likely molecular mechanism underlying the PatAGE effect. In a separate line of research, de novo mutation is known to increase risk by up to 20-fold in neurodevelopmental disorders (De Rubeis et al., 2014;Deciphering Developmental Disorders Study, 2017). Taken together, it is conceivable that de novo mutations are at least partially responsible for the negative effect of PatAGE on offspring's mental health, offering a plausible explanation of the PatAGE effect on children's reading abilities.
At the neurocognitive level, whether and how cognitive-linguistic factors mediate the PatAGE effect on reading development is unknown. Studies to date have focused on genetic and environmental factors that contribute to the multifactorial liability of dyslexia (Pennington, 2006;Petrill, Deater-Deckard, Thompson, DeThorne, & Schatschneider, 2006). One such example, phonological processing, is thought to exert its effect more dominantly through inherited genetic impact, often estimated by family reading history (van Bergen, Bishop, van Zuijen, & de Jong, 2015). Under the same framework, whether PatAGE serves as a contributor to the multifactorial liability, and if so, what the neural and cognitive mediators (that would likely be heritable but not inherited traits if caused by de novo mutation) are, have not been examined. Brain measures derived from neuroimaging techniques, including magnetic resonance imaging (MRI), are particularly informative as they have been suggested as mediators between genetic etiology and behavioral outcome, acting as endophenotypes (Grasby et al., 2020;Shaw et al., 2012). Further, longitudinal investigation, combined with cross-sectional analysis, can provide comprehensive insights into the neural basis underlying typical reading acquisition and its impairment (Clark et al., 2014;Yeatman, Dougherty, Ben-Shachar, & Wandell, 2012). Therefore, we conducted a preliminary study examining behavioral and multimodal neuroimaging data cross-sectionally and longitudinally in a cohort of 51 children followed from kindergarten through third grade in conjunction with analyses of publicly available datasets.
The objective of the study was threefold: (a) to examine the relationship between PatAGE and offspring's reading while systematically controlling for potential contributing/confounding factors; (b) to examine the role of previously known cognitive-linguistic precursors, neuroanatomy, and its maturational process in relation to PatAGE and reading; and (c) to understand the functional role of the neuroanatomical findings in this study by identifying convergent evidence through the use of brain atlases, public datasets, and offspring's diffusion imaging data.

| Participants
Participants in this study were drawn from a longitudinal NIH-funded project (K23HD054720) focusing on children's reading development and followed from kindergarten (time-point 1 [t1], mean age = 5.58 years, SD = 0.43) to third grade (time-point 2 [t2], mean age = 8.30 years, SD = 0.46). All children were healthy native English speakers without neurological/psychiatric disorders (e.g., attention deficit/hyperactivity disorder) or contraindications to MRI based on parental reports. Among the participating children, 76% were White, 6% were Asian, and 18% were of multiradical heritage. About 8% of the children identified as Hispanic or Latino. Based on the annual household income, parental educational levels, and occupation, the participants in this study were of relatively high socioeconomic status (also see Black et al., 2012). The initial sample consisted of 51 children recruited from local newspapers, school mailings, flyers, and mothers' clubs. In the behavior analyses, eight children were excluded because of attrition (n = 5), no record of PatAGE (n = 1), or being siblings (n = 2). In the latter case, the child with poorer T1 image quality was excluded. The final sample for behavioral analyses included 43 unrelated children (17 girls). In the neuroanatomical analysis, another seven children (two girls) were excluded because of incomplete T1 data collection or poor image quality at either time-point. In the diffusion imaging analysis, a sub-group of 23 children (eight girls) with the same acquisition sequence was included. The differences in either familial or behavioral measures between the entire and subgroups were nonsignificant (all p's > .1). The Institutional Review Boards of Stanford University where data were collected and the principal investigator was at the time of the study, and the University of California San Francisco where data were analyzed due to transition of the principal investigator, approved the present study. Both informed assent and consent were obtained from children and their guardians.

| Behavioral measurements
Demographics, family information, and performance on behavioral tests of the participants are summarized in Table 1. Family information collected at t1 includes: PatAGE; MatAGE; Adult Reading History Questionnaire from both parents (PatARHQ, MatARHQ) that were used to estimate familial history of reading difficulty (Lefly & Pennington, 2000); numbers of older and younger siblings; parental education (PatEDU, MatEDU); socioeconomic status (SES), a composite index computed from annual family income, parental education, and occupation with principal component analysis (Noble, Wolmetz, Ochs, Farah, & McCandliss, 2006); Home Observation Measurement of the Environment (HOME), an index for home environment including home literacy environment (Segers, Damhuis, van de Sande, & Verhoeven, 2016   were used to measure phonological skills. Rapid Automatized Naming (RAN; Objects and Colors sub-tests) (Wolf & Denkla, 2005) and Letter Identification sub-test of Woodcock Reading Mastery Test R/NU (WRMT-R/NU; Mather, 1998) were also administered.
The same set of tests was used at t2 (

| Image acquisition
High-resolution T1-weighted images were collected at both timepoints with the following parameters: 128 slices; thickness = 1.2 mm; NEX = 1; repetition time = 8.5 ms; echo time = 3.4 ms; inversion time = 400 ms; in-plane resolution = 256 Â 256; voxel size = 0.9 Â 0.9 Â 1.2 mm 3 ; flip angle = 15 ; field of view = 22 cm. Highangular resolution diffusion-imaging data were collected at t2 with the following parameters: 46 axial slices; slice thickness = 3 mm; repetition time = 5,000 ms; echo time = 81.7 ms; in-plane resolution = 128 Â 128; voxel size = 2.0 Â 2.0 Â 3.0 mm 3 ; 150 directions with b = 2,500 s/mm 2 ; 6 volumes with b = 0 s/mm 2 . All images were acquired using a GE Healthcare 3.0 T 750 scanner with eight-channel phased-array head coil at Richard M. Lucas Center for Imaging at Stanford University. The quality of images was qualitatively evaluated by an investigator who was blinded to the behavioral and demographic information prior to any analyses.

| Behavior analyses
To reduce the dimensionality of behavioral metrics, factor analyses were conducted on reading-related tests for each time-point; t1: Spelling sub-tests of WJ-III Tests of Achievement. The Maximum Likelihood, Varimax, and Bartlett methods were used for extraction, rotation, and factor score calculation. The criterion of eigenvalue greater than 1 was used to identify factors. From t1 behavioral metrics, we obtained two factors, explaining 53.8% of the total variance. Since phonological awareness (PA) and RAN loaded heavily on one factor and the other, respectively, we named these two factors t1PA and t1RAN ( Figure 1a). Given that PA and RAN are the most reliable predictors for reading development in alphabetic languages (Caravolas et al., 2012), we used these two composite scores as cognitivelinguistic precursors of reading in subsequent analyses. Using the same approach, we extracted three factors from t2 behavioral metrics, T A B L E 1 Demographic profiles, familial variables, and performance on reading-related tests (n = 43)  (Figure 1b).
Since a consensus on the definition of advanced PatAGE remains lacking (Couture, Delisle, Mercier, & Pennings, 2021), in this study we treated PatAGE as a continuous variable rather than separating children into different groups. To examine the relationship between PatAGE and reading, we first calculated the zero-order correlation.
Once the correlation was significant, hierarchical linear regressions were conducted to answer four questions in the following order: (a) whether the PatAGE effect on reading remains significant after controlling for demographic variables and general intelligence; (b) whether the PatAGE effect on reading exists after additionally regressing out MatAGE, which is known to correlate highly with PatAGE and is a possible confound; (c) whether the PatAGE effect on reading is present above and beyond familial risk (representing inherited risk; Swagerman et al., 2017), and environmental factors to address the issue of multifactorial liability; (d) whether the PatAGE effect on reading is explained by t1 cognitive-linguistic precursors of reading to examine whether the most common predictors were the mediating factor.
Specifically, in the first model we entered t2 age, sex, handedness, and average performance IQ (pIQ) across t1 and t2 in the first step and PatAGE in the second step (  (Edwards & Roff, 2010), number of older and younger siblings (Price, 2008), SES (Pan et al., 2017), HOME (Segers et al., 2016), which are known to be associated with reading were additionally controlled. In the final model, t1PA and t1RAN were entered in the fourth step, just before PatAGE, to examine whether the PatAGE effect was present beyond t1 cognitive-linguistic skills. Given that t1RAN and t1PA did not correlate with PatAGE (both r's < .01; Table S1), these two cognitive-linguistic precursors were not examined further for the mediating relationships. All statistics were done with SPSS 24.0 (IBM, Inc., Armonk, NY), and p-values were twotailed while statistical significance was set at .05.
The modulated images containing regional tissue volume of gray matter for each voxel were smoothed with an 8-mm full-width halfmaximum isotropic Gaussian kernel. Voxels with gray matter values < .1 were excluded to avoid edge effects.
In the longitudinal data preprocessing, the "Preprocessing of Longitudinal Data" module that contains specific steps for processing longitudinal structural MRI data was used. Intra-subject realignment, bias correction, segmentation, and normalization were conducted sequentially as described elsewhere (Ridgway et al., 2007). After applying spatial smoothing with an 8-mm full-width half-maximum Gaussian kernel, we obtained maps of gray matter volume for both time-points. We then generated a brain map reflecting gray matter volume (GMV) change from t1 to t2 for each child (such that a positive value indicates enlarging from t1 to t2).

| Whole-brain analyses
Prior to voxel-wise analyses, we examined relationships between PatAGE and the global measurement (i.e., total intracranial volume; defined as the sum of total gray matter, white matter, and cerebrospinal fluid) at each time-point (t1TIV and t2TIV). Then, we examined whether PatAGE correlated with the change of TIV from t1 to t2 (ΔTIV) while controlling for the baseline measure (t1TIV). To explore relationships between regional GMV at each time-point (i.e., crosssectional analyses), as well as the change of regional GMV (ΔGMV) across time-points (t2GMV-t1GMV) with PatAGE (i.e., longitudinal analysis), voxel-wise whole-brain regressions were conducted while controlling for global measurements. Specifically, t1TIV or t2TIV was controlled in cross-sectional analyses for t1 and t2, respectively. In the longitudinal analysis, t1TIV and ΔTIV were controlled to exclude the effects from initial gross volume and its development. Since correlations between t1TIV, ΔTIV, and PatAGE were not significant (all p's > .1), the models were free from multicollinearity. Topological Family Wise Error (FWE) correction was used to determine the corrected thresholds of statistical significance. All clusters significant at a threshold of p-cluster < .05 corrected for the whole brain (pvoxel < .001 for height) were reported in MNI space. Region-ofinterest (ROI) analyses were conducted in the significant clusters to examine the robustness and specificity of the effect. For this, values of each voxel in the cluster were extracted and averaged, then included in hierarchical multiple regression analyses as the dependent variable. Demographic variables (t1 or t2 age, sex, handedness, average pIQ across t1 and t2 for cross-sectional data; t1 age, time interval between t1 and t2, sex, handedness, average pIQ across t1 and t2 for longitudinal data) and global measurements (t1TIV or t2TIV for cross-sectional data; t1TIV and ΔTIV for longitudinal data) were entered in the first step, while PatAGE was entered in the second step. Then, we further controlled for MatAGE and MatARHQ since they showed significant correlations with PatAGE as in previous analyses of this study.

| Mediation analyses
One of the main objectives of this study was to investigate possible WRMT WID, Woodcock Reading Mastery Test, Word Identification sub-test brain level, two analytical strategies were used. For the primary approach, we conducted whole-brain analyses on t2READ crosssectionally and longitudinally in the same way as we did for PatAGE.
Next, we administered conjunction analysis to identify overlapping areas that showed significant associations with both PatAGE and t2READ, following which the mediation relationship was examined.
Alternatively, we took an ROI approach if no significant cluster survived multiple corrections at the whole-brain level on t2READ. Specifically, we examined the relationship between ΔGMV and children's t2READ in the cluster significantly associated with PatAGE (hereafter PatAGE-cluster). The partial correlation coefficient was calculated while controlling for demographic variables (t2 age, sex, handedness, average pIQ across t1 and t2), global measurements (t1TIV and ΔTIV), and cognitive-linguistic precursors (t1PA and t1RAN) that were significantly associated with t2READ in previous regression analysis of this study. Once a region was identified in either the whole-brain analysis or the ROI analysis, we examined whether the PatAGE effect on reading was mediated by brain measures. The model was adjusted for demographic variables (t2 age, sex, handedness, average pIQ across t1 and t2), global measurements (t1TIV, ΔTIV), and t1 cognitive-linguistic precursors (t1RAN and t1PA). Bootstrapping (10,000 samples) was used to obtain 95% confidence interval of the indirect effect. If the confidence interval does not contain zero, a significant indirect effect is indicated.
Existing evidence suggests that dyslexia is largely genetically transmitted from parents (often assayed by parent-report of reading difficulty; Soden et al., 2015;Swagerman et al., 2017). Further, twin studies find a dissociation between sources of variance in phonological and orthographic processes, with variance in phonological skills being primarily genetic compared to orthographic skills (Olson et al., 2011;Olson, Wise, Conners, Rack, & Fulker, 1989). These findings are consistent with the idea that PA mediates the negative effect of parental reading difficulty (a proxy for inherited genetic transmission) on reading in offspring (van Bergen et al., 2015). In line with the previous literature, we observed significant correlations between MatARHQ and t1PA, MatARHQ and t2READ, t1PA and t2READ. We therefore examined the role of PA on the relationship between the history of maternal reading difficulty and lower reading performance in offspring. Demographic variables (age at t2, sex, handedness, average pIQ across t1 and t2) and the other cognitivelinguistic precursor (t1RAN) were controlled statistically. If both a PatAGE effect (a proxy for non-inherited genetic risk) and MatARHQ effect via PA processing (a proxy for inherited genetic risk) are to be observed in the current samples, this study then supports the multifactorial liability model. PROCESS procedure (release 2.16.1) implemented in SPSS was used to conduct mediation analyses (Hayes, 2013).

| Complementary analyses
We adopted multiple complementary analytical approaches to depict fine-grained spatial localization and connectivity patterns of the T A B L E 2 Results of multiple linear regression analyses examining the unique contribution of paternal age on offspring's reading at timepoint 2

Model
Step Predictor ΔR 2 β PatAGE-cluster, capitalizing on the fact that these have been shown to inform possible functional roles of a particular brain region (in this case, the left posterior thalamus) in the absence of a comprehensive set of cognitive measures. First, we spatially localized the PatAGEcluster with two brain atlases. (a) Given that the thalamus consists of multiple nuclei with different functions, we calculated the percentage of overlapped voxels between the PatAGE-cluster and each thalamic nucleus from the Morel Atlas, a histological atlas that is optimal for thalamic targets in MNI space (Jakab, Blanc, Berényi, & Székely, 2012;Krauth et al., 2010); (b) Given that the connectivity pattern provides information about the function of a given brain region (Barron, Eickhoff, Clos, & Fox, 2015), we used the Oxford Thalamic Connectivity Probability Atlas (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Atlases) with the atlasquery tool implemented in FSL to obtain the probability that the PatAGE-cluster is structurally connected to different cortical areas.
To further identify the possible functional roles of the PatAGEcluster and complement the results from analyses using the histological and connectivity probability atlases, we examined PatAGE-clusterassociated cortical patterns by using online databases provided in Neurosynth v0.5 (Yarkoni, Poldrack, Nichols, Van Essen, & Wager, 2011). In particular, we generated a co-activation map by including all fMRI studies in the database (N > 10,900), with the PatAGE-cluster as the ROI. False Discovery Rate (FDR) corrected q < .01 was used as the threshold to obtain significant regions reported in fMRI studies when the PatAGE-cluster is also reported (i.e., forward inference). In addition, we generated a map of wholebrain resting-state functional connectivity (RSFC) by using the 1000 Functional Connectome dataset (Biswal et al., 2010). The center of gravity of the PatAGE-cluster (MNI: x = À19.6, y = À28.1, z = 6.9) was used as the seed, and its connectivity to the rest of the brain was calculated. The resultant brain map was thresholded with a liberal cutoff value of r ≥ .01, the same as in the previous literature (Yang, Rosenblau, Keifer, & Pelphrey, 2015). To be conservative, we only considered the overlapping regions between the co-activation map and the RSFC map as the cortical pattern of the PatAGE-cluster.
Sørensen-Dice coefficients (hereafter Dice coefficient) between the conjunction map and the seven large-scale intrinsic connectivity networks (i.e., visual, somatomotor, dorsal attention, ventral attention, limbic, frontoparietal, and default networks) from Yeo et al. (2011) were calculated to examine which functional network most overlapped with the PatAGE-cluster-associated cortical pattern. Here we used an adult network template because studies in Neurosynth that were used to produce the co-activation map were conducted in participants with a wide range of ages and the 1000 Functional Connectome dataset mainly consists of adult data.
Given that the functional community structure in children is to some extent different from that in adults and the highest uncertainty was found in attention networks (Tooley, Bassett, & Mackey, 2021;Vijayakumar et al., 2021), in the final step, we examined the structural connectivity pattern in a sub-group of our own samples where diffusion imaging data were available. Moreover, we adopted a pediatric intrinsic functional network template that was extracted from 670 children aged 9-11 years (Tooley et al., 2021) with the same approach as in Yeo et al. (2011). Specifically, we analyzed white matter connectivity, where fibers passing the PatAGE-cluster were reconstructed using deterministic tractography. Diffusion-weighted imaging preprocessing was performed by using ExploreDTI (Leemans, Jeurissen, Sijbers, & Jones, 2009). Next, to visually inspect for possible artifacts, rigorous motion correction with CATNAP and eddy current correction were conducted by using the required reorientation of the b-matrix ). The diffusion tensors were calculated using a nonlinear regression procedure (Pierpaoli & Basser, 1996). The individual datasets were nonrigidly normalized to MNI space. Next, whole-brain tractography was performed for each individual dataset using a deterministic approach ( Table S1). As to the main objective of this study, greater PatAGE was significantly correlated with lower reading composite scores in offspring (t2READ; r = À.39, p = .011).
To examine whether the PatAGE effect on reading existed above and beyond commonly identified confounds and additional variables known to influence reading development, hierarchical linear regressions were conducted with t2READ as the dependent variable in a systematic and hypothesis-driven fashion. In the first model, before PatAGE was entered, demographic variables (t2 age, sex, and handedness) and general intelligence (average pIQ across two time-points) were entered as predictors in the first step. The negative effect of advanced PatAGE remained significant, explaining an additional 14.9% of the variance (t = À3.12, p = .004; Model 1 in Table 2). In the second model, MatAGE was further added in the second step, since it was significantly correlated with PatAGE and t2READ. As shown in Table 2 Table 2). Thus far, we demonstrated that the PatAGE effect on reading was not accounted for by factors that predict children's reading outcomes and known to be either inherited or environmental. In the final model, we investigated its relationship with early cognitive-linguistic predictors of reading outcomes by entering t1PA and t1RAN in the fourth step. The PatAGE effect on offspring's reading was above and beyond that of cognitive-linguistic variables, explaining an additional 9.5% of the variance (t = À2.71, p = .014; Model 4 in Table 2; Figure S1B). In accord with the prior literature, t1PA and t1RAN also significantly predicted t2READ in the final model and accounted for 13.8% of the variation (t1PA: t = 2.87, p = .010; t1RAN: t = 2.19, p = .042). That is, contributions from PatAGE and cognitive-linguistic precursors were relatively independent, and jointly predicted children's reading outcomes.

| The PatAGE effect on offspring's reading is mediated by ΔGMV in the left posterior thalamus
First, no significant correlations were observed between t2READ with TIV at each time-point or ΔTIV from t1 to t2 (all p's > .1). In either the whole-brain cross-sectional or longitudinal analyses, we did not find clusters showing GMV-reading correlations survived the FWE correction. Therefore, the subsequent analyses were conducted based on the PatAGE-cluster revealed in the previous step. After we found ΔGMV in the PatAGE-cluster was correlated with t2READ while nuisance variables (t2 age, sex, handedness, average pIQ across t1 and t2), global measurements (t1TIV, ΔTIV), and cognitivelinguistic precursors of reading (t1PA, t1RAN) were statistically controlled (partial r = À.48, p = .011; Figure 2c), we further examined the mediation relationship. As shown in Figure 2d, ΔGMV significantly mediated the negative effect of advanced PatAGE on offspring's reading; 95% confidence interval was [À0.406, À0.004] when nuisance variables (age at t2, sex, handedness, average pIQ across t1 and t2), global measurements (TIV at t1, ΔTIV), and cognitive-linguistic precursors (t1PA, t1RAN) were statistically controlled. These results are in contrast to the commonly found results in the literature that we also find in the present study, that is, t1PA mediates the negative effect of family history on offspring's reading (95% confidence interval was [À0.249, À0.001] when nuisance variables (age at t2, sex, handedness, average pIQ across t1 and t2) and the other cognitive-linguistic precursor (t1RAN) were controlled; Figure S2).

| PatAGE-cluster is localized in the pulvinar nuclei and linked to the dorsal attention network
To understand the neurostructural profile of the PatAGE-cluster in the left thalamus, we compared it with a histological atlas and a con- Next, we examined the cortical pattern of the PatAGE-cluster by utilizing two approaches available in Neurosynth v0.5 (Yarkoni et al., 2011). These included the generation of a meta-analytic map of regions that co-activate with the PatAGE-cluster across more than 10,900 fMRI studies and an RSFC map from the PatAGE-cluster to the rest of the brain by using the 1000 Functional Connectome dataset (Biswal et al., 2010). The co-activated areas included subcortical structures and cortical regions such as bilateral intraparietal sulci, inferior temporal gyrus, and frontal eye fields in the frontal cortex ( Figure S3A). The RSFC map showed a similar but more widespread pattern than the co-activation map ( Figure S3B). A conjunction analysis revealed that the bilateral frontal eye fields, intraparietal sulci, middle temporal visual area (V5/MT), and cerebellum were among the overlapped regions across the two approaches, in addition to subcortical structures ( Figure S3C).
Dice coefficients between the overlapping areas and the previ- Together with the aforementioned findings utilizing structural atlases, these results using large-scale fMRI databases from functional neuroimaging studies point to the attention network, particularly the DAN, to be the candidate brain functional system associated with the PatAGE-cluster in the left thalamus. Finally, we analyzed diffusion imaging data available in a subgroup of 23 participants to confirm DAN was more likely the candidate system associated with the PatAGE effect on reading than VAN.
Because the diffusion imaging data were collected at time-point 2, different from the previous analyses with Neurosynth, here we adopted Tooley's functional network template that was extracted from data of children aged 9-11 years (Tooley et al., 2021). Using deterministic tractography, we reconstructed white matter fibers through the PatAGE-cluster, covering inferior fronto-occipital fasciculus, corticospinal tract, forceps major, superior corona radiata, as well as anterior and posterior limbs of the internal capsule. Figure 3d shows reconstructed fibers in a representative child and Figure 3e shows intersection across participants, for demonstrative purposes. In line with previous findings of this study, the PatAGE-cluster showed significantly stronger connectivity (defined as the total number of streamlines divided by global density of the target network [percentage of total voxels]) with DAN than with VAN (t = 5.24, p < .001; Figure 3f). No significant correlations were found between PatAGE-cluster-Network streamlines and PatAGE or t2READ (all p's > .1).

| DISCUSSION
In this study, we observed a significantly negative effect of PatAGE on offspring's reading at the earliest stages of formal schooling from ages 5 to 8, independent of confounds (e.g., maternal age) and factors that play key roles in learning to read (i.e., family reading history, environmental factors, and cognitive-linguistic precursors of reading), explaining an additional 9.5% of the variance. Furthermore, we revealed volumetric maturation of the left thalamus as a potential neural endophenotype mediating this effect. We identified that this area is most relevant to the dorsal attention network with brain atlases, public datasets, and offspring's diffusion imaging data. These findings contrast and complement the current literature linking phonological and orthographic processing in reading to the left temporo-parietal and occipito-temporal regions. The mediation revealed here was distinct from the mediating role of phonological processing on the relationship between reading and familial risk, which has been attributed to hereditary effects. Taken together, this study provides novel and converging evidence suggesting PatAGE as a significant factor that is associated with offspring's reading, independent of phonological processing, possibly through the maturational process of the left posterior thalamus.
4.1 | The PatAGE effect on offspring's reading Jayasekara and Street (1978), for the first time, reported that advanced PatAGE was associated with a greater incidence of dyslexia, independent of SES and birth order. While the analysis was conducted in dyslexic boys, Saha and colleagues extended the finding to a broader population of 7-year-old children with varying reading abilities measured using Wide Range Achievement Test (Saha et al., 2009). Negative PatAGE effects on several cognitive measures, including reading, were observed after controlling for confounds that included MatAGE, SES, and parental psychiatric illness. A follow-up study re-analyzed the same dataset and found that the PatAGE effect was no longer significant after further adjusting parental education and the number of siblings (Edwards & Roff, 2010). Therefore, the PatAGE effect on reading was equivocal, and the inconsistency was related to covariates controlled in the model, especially parental characteristics such as educational level.
In the present study, with the range of PatAGE restricted to 25-47 years, we found PatAGE was negatively associated with reading performance measured using a variety of tests, even after additionally controlling for strong predictors of reading that were not included in previous studies (Edwards & Roff, 2010;Saha et al., 2009). These predictors included family reading history and cognitive-linguistic skills (e.g., phonological processing) that shown to be more genetically than environmentally mediated, as well as home literacy environment (Hulme, Snowling, Caravolas, & Carroll, 2005;Swagerman et al., 2017;van Bergen et al., 2015). These findings support an adverse PatAGE effect on reading and suggest such effect may occur through a mechanism different to factors such as inherited genetic and environmental risks.
While the number of studies examining the PatAGE effect on reading is too few to infer the potential mechanisms, studies on PatAGE-linked neuropsychiatry disorders offer insights. One predominant explanation is that advanced PatAGE exerts its effect on the risk of disorders such as autism and schizophrenia through accumulated de novo genetic mutations and epigenetic modifications (e.g., DNA methylation and repressive histone modification) in paternal gametes (Deciphering Developmental Disorders Study, 2017;Girard et al., 2016;Saha et al., 2009).
At a more macroscopic scale, understanding of the mechanisms can be deepened by identifying intermediate (endo)phenotypes through behavioral and neuroimaging measures such as we did in the current study. That is, advanced PatAGE may impact precursors of neurodevelopmental disorders, which in turn leads to a higher occurrence of such disorders (Cannon, 2009). For example, the likelihood of having impaired social functioning in offspring, a core symptom of psychiatric disorders, increases with PatAGE (Weiser et al., 2008).
While the underlying mechanisms are yet to be fully understood, multifactorial liability confers risk for neurodevelopmental disorders and may involve liability such as de novo mutations in addition to inherited and environmental risks. Adding to prior research, the current study offers insights into potential mechanisms at the macroscopic level.

| The intermediary role of the left posterior thalamus
The thalamus is an important relay center in the human brain, receiving information from sensory cortices and relaying it to higher-level association cortices. Previous studies paint a mixed picture on thalamic development: gross volume relative to its brain size is smaller (Sussman, Leung, Chakravarty, Lerch, & Taylor, 2016) or larger (Brain Development Cooperative, 2012) in older compared to younger children of ages 4-18, and the pulvinar compared to other thalamic nuclei show no apparent change with age (Raznahan et al., 2014). Despite controversial evidence on typical thalamic maturation, its anomalous development undoubtedly affects the growth of other cortical and subcortical brain regions (Ball et al., 2012), which in turn could impact higher-level cognitive processes that underlie typical reading acquisition. In support of this, anomalies in thalamic structure (Giraldo-Chica,  (Yeatman et al., 2012).

| Limitations and future directions
In the present study, we observed a negative PatAGE effect on offspring's reading and the left posterior thalamus as a possible brain intermediary. Given the preliminary nature of this investigation and the small sample size, the findings should be interpreted with caution until they are replicated in large independent samples. Second, PatAGE in this study was restricted to 25-47 years, with which we observed a negative linear relationship. The findings hence do not necessarily extend to children with extremely young and old fathers.
For example, Saha et al. (2009) observed a nonlinear relationship with a range of PatAGE between 14 and 66 years, while the relationship appears to be linear in the age range as ours. Of relevance, young fatherhood is also associated with adverse outcomes in offspring but possibly due to other factors, including immature sperm and economic disadvantages (Chen et al., 2008). Third, the complementary analyses implied DAN as the candidate functional system associated with PatAGE. It should be noted that the brain atlases and public datasets implemented in Neurosynth are primarily from research on adults, while children have specific characteristics regarding brain organization (Tooley et al., 2021;Vijayakumar et al., 2021). We attempted to address this issue by adopting the functional network template for children in analyzing diffusion imaging data available in a subgroup of the current samples and found that the PatAGE-cluster more strongly connected with DAN than VAN. Nevertheless, the neurostructural profile and function of the PatAGE-cluster need to be re-visited as more pediatric-specific atlases and tools are available. Fourth, while we found that the left posterior thalamus mediated the PatAGE effect on reading, it remains unknown why this subcortical structure is susceptible to advanced PatAGE (and related to de novo mutations).
Given that typical thalamic maturation is also affected by prenatal and postnatal factors such as preterm birth (Ball et al., 2012), questions including how PatAGE influences the growth of thalamus and relevant functional systems, together with other factors, require further elaboration.
To advance understanding of the PatAGE effect on reading, future studies are warranted in which a more comprehensive battery of behavioral tests (e.g., measuring visuo-spatial attention, executive function, etc.), neural measures (e.g., task-driven activation), and molecular approaches measuring the number and origins of de novo mutations (e.g., trio-based whole-genome/exome sequencing; Jin et al., 2018) are included. Fusing neural, cognitive, and molecular genetic approaches at multiple levels will provide the much-needed vertical and multi-level explanatory models that will further our understanding of risk factors associated with poor reading. In particular, future research aiming at disentangling different sources of genetic variations related to reading development and their interplays will greatly further our understanding.
In addition, advanced research designs such as the intergenerational neuroimaging approach can be adopted to gain in-depth knowledge on how multiple factors, including PatAGE, affect the development of offspring's reading and the corresponding networks interactively from preliteracy to mature stages of reading (Ho, Sanders, Gotlib, & Hoeft, 2016;Hoeft & Hancock, 2018).

| CONCLUSION
The current study examined the PatAGE effect on offspring's reading at both behavioral and neurobiological levels. The results provide initial evidence that advanced PatAGE is a relatively independent factor associated with poor reading outcomes in beginning readers, above and beyond previously identified familial and cognitive-linguistic precursors.
This effect was mediated by the maturation of the posterior thalamus, suggesting a neurobiological pathway to intergenerational influence on reading acquisition, complementing prior findings that offspring's reading is influenced by parental reading via offspring's phonological skills (van Bergen et al., 2015). Based on these findings, we argue that PatAGE should be regarded as an important factor influencing literacy development, and included in a cumulative risk (and protection) model

ACKNOWLEDGMENTS
The authors thank all the families for their participation in this longitudinal study. They also thank Albert Galaburda and Tuong-Vi Nguyen for their thoughtful suggestions during manuscript preparation.

CONFLICT OF INTEREST
The authors declare no competing financial interests.