Neural network correlates of high-altitude adaptive genetic variants in Tibetans: A pilot, exploratory study.

Although substantial progress has been made in the identification of genetic substrates underlying physiology, neuropsychology, and brain organization, the genotype-phenotype associations remain largely unknown in the context of high-altitude (HA) adaptation. Here, we related HA adaptive genetic variants in three gene loci (EGLN1, EPAS1, and PPARA) to interindividual variance in a set of physiological characteristics, neuropsychological tests, and topological attributes of large-scale structural and functional brain networks in 135 indigenous Tibetan highlanders. Analyses of individual HA adaptive single-nucleotide polymorphisms (SNPs) revealed that specific SNPs selectively modulated physiological characteristics (erythrocyte level, ratio between forced expiratory volume in the first second to forced vital capacity, arterial oxygen saturation, and heart rate) and structural network centrality (the left anterior orbital gyrus) with no effects on neuropsychology or functional brain networks. Further analyses of genetic adaptive scores, which summarized the overall degree of genetic adaptation to HA, revealed significant correlations only with structural brain networks with respect to local interconnectivity of the whole networks, intermodule communication between the right frontal and parietal module and the left occipital module, nodal centrality in several frontal regions, and connectivity strength of a subnetwork predominantly involving in intramodule edges in the right temporal and occipital module. Moreover, the associations were dependent on gene loci, weight types, or topological scales. Together, these findings shed new light on genotype-phenotype interactions under HA hypoxia and have important implications for developing new strategies to optimize organism and tissue responses to chronic hypoxia induced by extreme environments or diseases.

INTRODUCTION

The Tibetan Plateau, also known as the Qinghai‐Tibet Plateau in China, has an average elevation of 4,000 m with ~40% lower oxygen pressure than sea level and is thus very harsh for human habitation. Nonetheless, indigenous Tibetans have adapted to survive and reproduce for millennia in such hostile high‐altitude (HA) environments. The adaption is accompanied by multidimensional alterations ranging from genotype to a variety of phenotypes, presumably as a consequence of evolution by process of natural selection to beneficially adjust to surrounding environments.

Currently, one of the most studied and best documented HA adaptive changes is genetic variants. To date, EGLN1 and EPAS1 are the two best‐established genetic loci related to HA adaption (Simonson et al., 2010; Xu et al., 2011), both of these loci are crucial components of the hypoxia‐inducible transcription pathway. Specifically, EGLN1 encodes prolyl hydroxylase 2, a principal negative regulator of hypoxia‐inducible factors (HIFs), while EPAS1 encodes hypoxia‐inducible factor 2α (HIF2α), a transcription factor involved in the body's response to hypoxia. In addition to these two genes, PPARA was recently identified as another top candidate gene that is implicated in HA adaptation by regulating fatty acid metabolism (Murray, Montgomery, Feelisch, Grocott, & Martin, 2018). PPARA encodes peroxisome proliferator‐activated receptor α, a member of the PPAR family of ligand‐activated transcription factors that play key roles in regulating cellular energy metabolism.

Regarding phenotype, respiratory and circulatory systems are the two most salient physiological processes to manifest HA adaptive alterations, including larger chest circumference, greater lung diffusing capacity, and greater maximum heart rate in highlanders compared with those in lowlanders (Gilbert‐Kawai, Milledge, Grocott, & Martin, 2014). Following these physiological alterations, brain structure and function are inevitably affected due to the high oxygen consumption of the brain, such as modifications in local brain features (Wei, Wang, Gong, Fan, & Zhang, 2017; Zhang, Zhang, et al., 2013), intraregional neural synchronization (Chen, Fan, et al., 2016), and interregional functional and structural connectivity (Chen, Li, et al., 2016). Finally, given the numerous brain alterations, neuropsychological adaption to HA is unsurprising since the brain is the basis and carrier of behavior and cognition (Smith et al., 2015).

Despite the advances in the repertoire of HA adaptive alterations, the observed adaption is not completely understood due to limited knowledge on how HA adaptive genetic variants relate to phenotypic alterations. In this regard, existing research mainly focuses on the genetic modulation of physiological characteristics. For example, HA adaptive variants in both EGLN1 and EPAS1 are associated with lower hemoglobin concentrations in Tibetan highlanders (Beall et al., 2010; Tashi et al., 2017). In recent years, studies relating genotype to other dimensions of phenotype, brain structure and function in particular, have become increasingly popular. For example, based on a dataset of 8,428 healthy subjects, a recent study carried out genome‐wide association analyses of 3,144 structural and functional brain imaging phenotypes and identified 148 clusters of associations between single‐nucleotide polymorphisms (SNPs) and imaging phenotypes (Elliott et al., 2018). However, despite the prevalence of genotype–phenotype association studies, to the best of our knowledge, no study has systematically examined the associations of HA adaptive genetic variants with different levels of phenotype in Tibetan highlanders.

To bridge this gap, in the current study we conducted a comprehensive examination of how HA adaptive genetic variants relate to interindividual variance in different levels of phenotype, ranging from physiology to neuropsychology and from brain structure to brain function, in indigenous Tibetan highlanders. Specifically, we collected blood and multimodal MRI data from 135 Tibetan natives who also underwent a comprehensive physiological and neuropsychological test battery. We first identified HA adaptive SNPs among 60 candidate SNPs from three gene loci (EGLN1, EPAS1, and PPARA) that were previously reported to undergo strong positive selection for HA adaptation (Beall et al., 2010; Bigham et al., 2010; Ge, Simonson, Gordeuk, Prchal, & McClain, 2015; Hu et al., 2017; Lorenzo et al., 2014; Peng et al., 2011; Peng et al., 2017; Tashi et al., 2017; Xiang et al., 2013). Then, we constructed and topologically characterized large‐scale structural and functional brain networks from the multimodal MRI data given that the human brain operates essentially as an interconnected network wherein neural signals propagate and exchange to respond to various inputs from different sensory systems (Bullmore & Sporns, 2009). Moreover, a substantial body of evidence indicates that the human brain network is related to interindividual variance in cognition and behavior (Park & Friston, 2013), and disturbances in its normal organization are largely responsible for cognitive and behavioral dysfunction in brain disorders (Rubinov & Bullmore, 2013). Therefore, uncovering genetic bases of the human brain networks is of great significance for understanding gene–brain–behavior relationships, and for providing insights into the genetic architecture of the brain that are relevant to neurological and psychiatric disorders, brain development, and aging. Finally, we examined effects of HA adaptive genetic variants on different levels of phenotype for each HA adaptive SNP via dichotomous comparisons and for all HA adaptive SNPs as a whole with correlation analyses. In addition, we explored pairwise relationships between different levels of phenotype. We hypothesized that HA adaptive genetic variants would selectively affect specific phenotypes in a complicated manner, particularly for large‐scale brain network organization.

MATERIALS AND METHODS

Participants

A total of 135 native Tibetan students (62 males and 73 females; aged from 17 to 23 years with a mean = 19.9 ± 1.1) from the Tibet University at Lhasa (3,650 m, Tibet, China) were enrolled in this study. To increase sample homogeneity, all Tibetans were Tibetan Buddhists living in the Tibetan Plateau (2,600–4,700 m with a mean = 3,994.5 ± 326.7) with native Tibetan ancestors. Both the Tibetans and their ancestors had no prior descent to the lowlands. None of the Tibetan participants had previously suffered from mountain sickness. In addition, 65 Han Chinese students (36 males and 29 females; aged from 18 to 25 years with a mean = 20.4 ± 2.3) from Xiamen University (sea level, Xiamen, China) were recruited to determine HA adaptive genetic variants. All Han Chinese students were lowlanders who were born and lived in lowlands (<500 m) and had no prior exposure to HA. All participants were right‐handed and had no history of neurologic disorder, head injury, surgery, drug or alcohol dependence, smoking, or genetic disease. The experimental protocol was approved by the Research Ethics Review Board of Xiamen University, and each participant gave written informed consent before their participation in this study.

Genotyping

Whole blood was collected for each participant in the morning between 07:00 and 07:30 for DNA isolation and genotyping. All genotyping was conducted by the Beijing Genomics Institute. DNA samples were extracted from the blood of the 200 participants and were examined by laboratory personnel who were blinded to individual details. Genomic DNA was extracted from peripheral blood using a DNA extraction kit (TianGen Biotech Beijing, CO., LTD, China) and stored at −80°C.

Sixty candidate SNPs from three genes (EPAS1, EGLN1, and PPARA) were chosen a priori for genotyping in the current study according to previous studies. The genotyping was performed by matrix‐assisted laser desorption/ionization time‐of‐flight mass spectroscopy using a 384‐format high‐throughput Sequenom MassARRAY platform (MassARRAY Compact Analyzer, Sequenom, San Diego, CA). Sequenom's SpectroDesigner software was used for SNP assay design, and SpectroTyper 4.0 was used to call genotypes automatically, followed by manual review. To confirm the accuracy of the genotyping results, we randomly selected 10% of samples for each SNP for testing twice by different people, which resulted in 100% reproducibility. According to the genotyping results, one Tibetan participant was removed due to a sample detection rate < 90%. Of the 60 SNPs, seven were excluded from further analysis because of failure to meet a genotype call rate > 90% (3) or the Hardy–Weinberg equilibrium p > .001 (3) or to pass a minor allele frequency > 0.05 (1).

Physiological and neuropsychological assessments

All Tibetan participants underwent a comprehensive physiological and neuropsychological test battery. Specifically, the physiological assessments included percutaneous arterial oxygen saturation (SpO2), heart rate, blood pressure (systolic and diastolic), hematology (erythrocyte and hemoglobin), and pulmonary function (forced vital capacity [FVC], forced expiratory volume in the first second [FEV1], vital capacity [VC], and FEV1/FVC). The neuropsychological tests included the Chinese revised version of the Wechsler Memory Scale subset test (1 → 100, 100 → 1, accumulation, figure memory, figure recognition, figure construction, touch score, and digit span backward task) and the Rey‐Osterrieth Complex Figure test (immediate memory and long‐term memory).

Multimodal MRI data acquisition

In this study, we aim to relate HA adaptive genetic variants to interindividual variance in the topological organization of large‐scale brain networks in a single group of Tibetan highlanders. Therefore, we collected multimodal MRI data (structural, diffusion, and resting‐state functional) from each Tibetan participant on a 3.0 Tesla scanner (Siemens Magnetom Trio Tim system, Erlangen, Germany) at the MRI Center in Lhasa (Tibet Autonomous Region People's Hospital, Lhasa, China). We did not collect MRI data from the Han Chinese students, who were used only to screen out the HA adaptive SNPs.

Structural MRI

High‐resolution structural images were acquired using a magnetization‐prepared two rapid acquisition gradient echo (MP2RAGE) sequence. The parameters are as follows: slice number = 176, repetition time (TR) = 5,000 ms, echo time (TE) = 2.98 ms, field of view (FOV) = 240 × 256 mm2, matrix = 240 × 256, flip angle = 9°, voxel size = 1 × 1 × 1 mm3, slice thickness = 1. 0 mm, and no gap.

Diffusion MRI

The diffusion tensor imaging data were acquired using a single‐shot echo‐planar imaging‐based sequence with the following parameters: slice number = 40, TR = 5,400 ms, TE = 95 ms, FOV = 220 × 220 mm2, matrix = 128 × 128, flip angle = 90°, slice thickness = 3 mm, and gap = 0.3 mm. Thirty diffusion‐weighted directions with b = 1,000 s/mm2 and 1 b0 image (b = 0 s/mm2) were obtained.

Resting‐state functional MRI

Resting‐state functional images were collected using an echo‐planar imaging sequence: TR = 2,500 ms, TE = 30 ms, flip angle = 90°, slice thickness = 3 mm, gap = 1 mm, in‐plane resolution = 2.34 × 2.34 mm2, and matrix = 94 × 94. During the scanning, participants were instructed to keep their eyes closed, relax their minds, and remain motionless as much as possible but to avoid falling asleep.

All images were carefully inspected visually, and six participants were excluded from further analysis due to poor quality for at least one modality of the MRI data. All participants had normal MRI as determined by an experienced radiologist.

Multimodal MRI data preprocessing

The diffusion‐weighted MRI data were processed using the PANDA toolbox (Cui, Zhong, Xu, He, & Gong, 2013) based on FSL (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/). First, individual diffusion‐weighted images were coregistered to their corresponding b0 images using an affine transformation to correct for head motion and eddy current‐induced distortions. Then, a diffusion tensor model was estimated at each voxel, and diagonalization was performed to yield three eigenvalues and eigenvectors (Basser & Pierpaoli, 1996).

The resting‐state functional MRI data were preprocessed using the GRETNA toolbox (Wang, Wang, Xia, et al., 2015) based on SPM12 (https://www.fil.ion.ucl.ac.uk/spm/software/spm12/). After discarding the first five volumes, individual functional images were corrected for intravolume temporal offsets (Sinc interpolation) and intervolume head motion (rigid‐body transformation). Nine participants were excluded from further analysis of functional brain networks due to excess head motion in terms of the criterion of a displacement >2 mm, an angular rotation >2° or a mean framewise displacement >0.25 mm. The corrected images were then spatially normalized into the standard Montreal Neurological Institute (MNI) space (via the tissue segmentation of anatomical images) and resampled into 3‐mm isotropic voxels. The normalized images subsequently underwent removal of linear trend and temporal bandpass filtering (0.01–0.08 Hz). Finally, several nuisance signals were regressed out from each voxel's time series to reduce non‐neuronal sources, including head motion profiles, white matter signals, cerebrospinal fluid signals, and global signals. Specifically, the head motion profiles included 24 parameters to reduce residual head motion effects on subsequent network indices (Friston, Williams, Howard, Frackowiak, & Turner, 1996). Notably, all nuisance signals were also bandpass filtered (0.01–0.08 Hz; Hallquist, Hwang, & Luna, 2013).

Multimodal brain network construction

A network is composed of a set of nodes linked by edges. In the current study, we constructed weighted structural and functional brain networks for each participant with nodes representing brain regions, and edges representing interregional structural or functional connectivity.

Brain parcellation for network node definition

To define brain network nodes, we employed the automated anatomical labeling (AAL) parcellation atlas to divide the brain into 94 regions of interest (ROIs, 47 in each hemisphere; Table 1; Rolls, 2015). Specifically, the regional parcellation was performed in the MNI space for the functional data where the atlas was originally released while conducted in individual native diffusion spaces for the diffusion data where whole‐brain white matter fiber tracts were reconstructed. For the latter, an affine transformation (coregistration of individual fractional anisotropy [FA] images in native diffusion space to the corresponding T1 images) and a nonlinear transformation (registration of the individual T1 images to the ICBM152 T1 template in the MNI space) were first estimated for each participant. Based on these two transformations, an inverse transformation from the MNI space to each participant's native diffusion space was then calculated and applied to the AAL atlas.

Table 1

Regions of interest

Index	Regions	Abbreviations	Index	Regions	Abbreviations
1,2	Precentral gyrus	PreCG	49,50	Cuneus	CUN
3,4	Superior frontal gyrus, dorsolateral	SFGdor	51,52	Lingual gyrus	LING
5,6	Middle frontal gyrus	MFG	53,54	Superior occipital gyrus	SOG
7,8	Inferior frontal gyrus, opercular part	IFGoperc	55,56	Middle occipital gyrus	MOG
9, 10	Inferior frontal gyrus, triangular part	IFGtriang	57,58	Inferior occipital gyrus	IOG
11,12	Inferior frontal gyrus, orbital part	ORBinf	59,60	Fusiform gyrus	FFG
13,14	Rolandic operculum	ROL	61,62	Postcentral gyrus	PoCG
15,16	Supplementary motor area	SMA	63,64	Superior parietal gyrus	SPG
17,18	Olfactory cortex	OLF	65,66	Inferior parietal gyrus, excluding supramarginal and angular gyri	IPL
19,20	Superior frontal gyrus, medial	SFGmed	67,68	Supramarginal gyrus	SMG
21,22	Superior frontal gyrus, medial orbital	ORBsupmed	69,70	Angular gyrus	ANG
23,24	Gyrus rectus	REC	71,72	Precuneus	PCUN
25,26	Medial orbital gyrus	OFCmed	73,74	Paracentral lobule	PCL
27,28	Anterior orbital gyrus	OFCant	75,76	Caudate nucleus	CAU
29,30	Posterior orbital gyrus	OFCpost	77,78	Lenticular nucleus, putamen	PUT
31,32	Lateral orbital gyrus	OFClat	79,80	Lenticular nucleus, pallidum	PAL
33,34	Insula	INS	81,82	Thalamus	THA
35,36	Anterior cingulate and paracingulate gyri	ACC	83,84	Heschl gyrus	HES
37,38	Middle cingulate and paracingulate gyri	MCC	85,86	Superior temporal gyrus	STG
39,40	Posterior cingulate gyrus	PCC	87,88	Temporal pole: Superior temporal gyrus	TPOsup
41,42	Hippocampus	HIP	89,90	Middle temporal gyrus	MTG
43,44	Parahippocampal gyrus	PHG	91,92	Temporal pole: Middle temporal gyrus	TPOmid
45,46	Amygdala	AMYG	93,94	Inferior temporal gyrus	ITG
47,48	Calcarine fissure and surrounding cortex	CAL

Note: The regions are listed in terms of a prior AAL atlas. Regions of left and right hemispheres are indexed by odd and even numbers, respectively.

Connectivity estimation for network edge definition

For the diffusion data, whole‐brain white matter fiber tracts were reconstructed for each participant using a continuous streamline‐tracking algorithm (Mori, Crain, Chacko, & van Zijl, 1999). The tractography was terminated when it reached a voxel with an FA less than 0.2 or when the turning angle was greater than 45° between adjacent steps. For each region pair, the number of fiber streamlines with two end‐points located in these two regions was counted and defined as interregional connectivity weights. In addition to the fiber number (FN), we calculated the mean FA and fiber length (FL) over all fiber streamlines between each pair of regions to derive individual FA‐weighted and FL‐weighted networks, which allow a full characterization of genetic modulation of structural brain networks in Tibetan participants.

For the resting‐state functional data, the mean time series was extracted for each brain region by averaging the signals across voxels within that region. The resultant mean time series were then correlated with each other to generate a 94 × 94 correlation matrix and a 94 × 94 p‐value matrix for each participant. To exclude confounding effects of noisy interregional connectivity, a threshold of p < .05 (Bonferroni corrected across 94 × [94 – 1]/2 = 4,391 connections) was applied to each p‐value matrix to screen out significant connections. Finally, only correlation coefficients with p values that survived the p‐value threshold were retained in the functional brain networks. Such a significance‐based thresholding procedure has been used in previous studies (Bassett et al., 2011; Wang, Wang, He, et al., 2015). Additionally, negative correlations were excluded due to their ambiguous interpretation (Fox, Zhang, Snyder, & Raichle, 2009) and detrimental effects on test–retest reliability (Wang et al., 2011).

Of note, both structural and functional brain networks derived above differed in the number of connections (i.e., network density) across participants. Although previous studies have shown that different network densities can confound between‐group comparisons (Ginestet, Nichols, Bullmore, & Simmons, 2011; van Wijk, Stam, & Daffertshofer, 2010), this study aims to relate HA adaptive genetic variants to interindividual variance in topological organization of large‐scale brain networks in a single group of Tibetan highlanders. Accordingly, instead of enforcing the same network density across participants, we retained inherent network densities of both structural and functional brain networks for each participant.

Graph‐based network analysis

For the constructed structural and functional brain networks, we characterized their topological organizations at multiple levels using the GRETNA toolbox (Wang, Wang, Xia, et al., 2015). Details on the formulas, uses, and interpretations of the network measures employed in the current study are described elsewhere (Rubinov & Sporns, 2010).

Global whole‐brain organization

At the global level, we calculated five metrics for each network, including the clustering coefficient (C p), characteristic path length (L p), local efficiency (E loc), global efficiency (E glob), and modularity (Q). The clustering coefficient quantifies local interconnectivity or cliquishness of a network, the characteristic path length measures the mean distance or routing efficiency between any pair of nodes in a network, the local efficiency reflects the extent to which a network is fault tolerant, the global efficiency indicates how efficiently information can be exchanged in parallel over a network, and the modularity is an indicator of the extent to which nodes can be divided into several subsets with dense edges within but sparse edges between them.

Intermediate modular composition

Modular structure is an important organizational principle for human brain networks (Sporns & Betzel, 2016). To determine the modular compositions of structural and functional brain networks in this study, we first generated a group‐level matrix for each modality by counting the frequency of occurrence for each connection over all participants. Then, a backbone network that captured the multiscale structures of the group‐level matrix was extracted (Foti, Hughes, & Rockmore, 2011) and used to detect the modular architecture based on a spectral optimization algorithm (Newman, 2006). With the resultant group‐level modular composition as a representative reference, we further calculated the intramodule integration score (mean connectivity strength within a module), intermodule communication score (mean connectivity strength between two modules), and nodal participant coefficient (an index reflecting the ability of a node to maintain communication between its own module and other modules) for each participant.

Local nodal centrality

For each node, we calculated five centrality metrics, including nodal degree (k ), nodal efficiency (e ), nodal betweenness (b ), nodal eigenvector (v ), and nodal pagerank (p ) to characterize their roles in a network. The nodal degree is a simple metric to measure the overall connectivity of a node with other nodes in a network, the nodal efficiency quantifies the ability of information propagation of a node with other nodes in a network, the nodal betweenness captures the influence of a node over information flow between all other nodes in a network, the nodal eigenvector reflects the importance of node by taking the centralities of its neighbors into account and thus is sensitive to different layers in a network hierarchy, and the nodal pagerank is a variant of by introducing a damping factor to handle walking traps of a network.

Statistical analysis

Between‐group differences in demographic variables

Two‐sample t test was used to infer differences in age, and χ 2 tests were used to assess distribution differences in sex and grade between Tibetan and Han Chinese participants.

Identification of HA adaptive SNPs

For each of the 53 SNPs included in the final analysis, a χ 2 test was used to examine the allele frequency distribution between the two population samples to identify HA adaptive SNPs. A false discovery rate (FDR) procedure was utilized to correct for the multiple comparison issue at a level of q < 0.05 (53 comparisons). For each SNP showing a significant between‐group difference, the HA adaptive allele that exhibited a higher frequency in the Tibetans was identified further.

Effect of a single HA adaptive SNP

For each HA adaptive SNP, the 128 Tibetan participants were first divided into two groups, with one including homozygotes of HA adaptive alleles and the other including homozygotes of HA nonadaptive alleles. The two groups were hereafter termed the HA adaptive group and the HA nonadaptive group, respectively. Due to the unbalanced sample size between the two groups, the heterozygotes were included in the group with fewer participants to promote comparability. Then, a nonparametric permutation test was used to examine differences in the physiological characteristics, neuropsychological variables, and brain network attributes between the two groups. Briefly, for each metric, we initially calculated the between‐group difference in their mean values. An empirical distribution of the difference was subsequently obtained by randomly reallocating all of the values into two groups and recomputing the mean differences between the two randomized groups (100,000 permutations). The 95th percentile points of the empirical distribution were used as critical values in a one‐tailed test to determine whether the observed real‐group difference occurred by chance. Notably, prior to the permutation tests, multiple linear regressions were used to remove the effects of age, sex, grade, and altitude. The FDR procedure was separately used to correct for multiple comparisons for physiological characteristics (495 comparisons), neuropsychological variables (450 comparisons), and brain network attributes (27,090 comparisons for functional networks and 26,775 comparisons for each type of structural networks).

Holistic effect of all HA adaptive SNPs

Given the possibility that a single HA adaptive SNP has little effect on the global phenotype due to the lack of one‐to‐one mapping, we further examined the holistic effect of all HA adaptive SNPs as a whole. In doing so, we defined a genetic adaptive score (GAS) for each participant to assess the degree of cumulative genetic adaptation. GAS may show stronger association with phenotype than individual marker because they contain more information about the underlying functional genetic variants. The GAS is similar to the genic risk score that is frequently used to estimate cumulative genetic risk of a subject for a disease (Xu et al., 2018). Specifically, a weight of 1, 2, or 3 was first assigned to each HA adaptive SNP of each participant, depending on 0, 1, or 2 HA adaptive alleles carried by the SNPs. The GAS was then defined as the summation of the weights across all HA adaptive SNPs. To examine the relationships between the GAS and physiological characteristics, neuropsychological variables, and brain network attributes, we used a nonparametric Spearman correlation with age, sex, grade, and altitude as nuisance covariates. Multiple comparisons were corrected with the FDR procedure separately for physiological characteristics (11 comparisons), neuropsychological variables (10 comparisons), and brain network attributes (602 comparisons for functional networks and 595 comparisons for each type of structural networks).

To further test relationships between the GAS and interregional connectivity, we utilized a cluster‐based statistics procedure, which was originally proposed to identify a set of subnetworks consisting of multiple connected connections, rather than a single connection, that are strongly correlated with a specific cognitive behavior (Han, Yoo, Seo, Na, & Seong, 2013). In brief, a Spearman correlation test was first conducted for each connection with the GAS after controlling for the effects of age, sex, grade, and altitude. Then, an initial cluster forming threshold (p < .05) was used to identify suprathreshold connections among which all connected components and their associated sizes (i.e., number of connections) were identified and recorded. Finally, a null distribution of maximal connected component size was empirically derived using a nonparametric permutation approach (10,000 permutations) and used to estimate the significance of each component. Specifically, for a real connected component of size M, the corrected p value was determined by calculating the proportion of the 10,000 permutations for which the maximal connected component size was larger than M. Notably, only connections that existed in more than half of the participants were included in the correlation analyses.

Finally, we explored the possible genetic dependence of the correlation analyses above by calculating the GAS across the HA adaptive SNPs in each of EGLN1, EPAS1, and PPARA.

Relationships among physiology, neuropsychology, and brain network organization

In addition to the effects of HA adaptive SNPs, we examined the relationships between physiological characteristics, neuropsychological variables, and brain network attributes with the Spearman correlation after controlling for age, sex, grade, and altitude. Multiple comparisons were corrected with the FDR procedure separately for correlations between physiological characteristics and neuropsychological variables (110 comparisons), between physiological characteristics and network measures (6,622 comparisons for functional networks and 6,545 comparisons for each type of structural networks), between neuropsychological variables and network measures (6,020 comparisons for functional networks and 5,950 comparisons for each type of structural networks), and between network measures of any two types of brain networks (574 comparisons).

Validation analysis

We performed the following validation analyses to test the reproducibility of our findings of functional brain networks given the lack of consensus on the processing pipeline of resting‐state functional MRI data.

Head motion

In this study, we combined a series of strategies to minimize the effects of head motion on our results, including removal of participants with excess gross head motion (>2 mm in translation or >2° in rotation), removal of participants with large, rapid changes in head motion (i.e., spikes; >0.25 mm in mean framewise displacement), regression of 24‐parameter head motion profiles (Yan, Cheung, et al., 2013), and regression of global signals (Yan, Craddock, He, & Milham, 2013). It should be pointed out that there are numerous methods developed recently to control head motion for resting‐state functional MRI data. Therefore, we also re‐analyzed our data using a scrubbing approach (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012) to censor “bad” volumes (framewise displacement >0.25 mm) and volumes temporally close to the “bad” volumes (1 before and 2 after) for each participant before constructing individual functional networks.

Global signal

In this study, we regressed out global signals to control effect of head motion on functional brain networks (Yan, Craddock, et al., 2013). Currently, whether global signals should be regressed out is an ongoing topic of research, and different choices likely reveal complementary insights into the brain's functional organization (Murphy & Fox, 2017). Therefore, we re‐analyzed our data that did not undergo global signal regression.

Brain parcellation

How to define network nodes is an ongoing topic of research for brain network studies thus far. Currently, there are numerous brain parcellation methods available and there is no optimal choice (Arslan et al., 2018). Generally speaking, anatomical parcellation atlases are typically suitable for structural brain network studies while functionally defined atlases are more suitable for functional brain network studies. In this study, we constructed both structural and functional brain networks; it is thus hard to determine the “best” way for network node definition. Accordingly, we selected one of the most widely used atlases (i.e., AAL atlas) in previous literature. However, recent studies indicate that this atlas cannot reflect functional organization of the brain (Eickhoff, Constable, & Yeo, 2018). Therefore, we reconstructed functional brain networks using a functionally defined atlas (Dosenbach et al., 2010). This atlas includes 160 ROIs (5‐mm radius spheres) that are identified from meta‐analyses of task‐related fMRI studies.

RESULTS

Demographic characteristics

There were no significant differences in age, sex, or grade between Tibetan and Han Chinese participants (p > .05). The detailed demographic, physiological, and neuropsychological information for Tibetan participants are summarized in Table 2.

Table 2

Demographic, physiological, and neuropsychological information of all Tibetan participants

	N	Value
Demographic
Sex (M/F)	135	62/73
Age (years)	135	19.9 ± 1.1
Height (cm)	135	163.9 ± 7.2
Weight (kg)	135	54.0 ± 6.4
HA location (m)	135	3,994.5 ± 326.7
Education (years)	135	12.7 ± 0.5
Physiological
SpO2 (%)	88	92.7 ± 3.2
Erythrocyte (×1012/l)	128	5.1 ± 0.5
Hemoglobin (g/l)	128	152.9 ± 18.9
Heart rate (beats/min)	88	70.1 ± 9.9
Blood pressure (mmHg)
Systolic pressure	86	111.6 ± 9.9
Diastolic pressure	86	66.4 ± 7.8
Pulmonary function
VC (ml)	70	2,906.3 ± 896.1
FVC (ml)	90	2,870.9 ± 580.5
FEV1 (ml)	90	2,525.9 ± 636.1
FEV1/FVC (%)	90	88.1 ± 14.1
Neuropsychological
Wechsler memory scale test
1–100	123	7.4 ± 3.2
100–1	123	11.5 ± 1.7
Accumulation	123	10.9 ± 1.6
Figure memory	123	10.9 ± 1.6
Figure recognition	123	10.5 ± 2.2
Figure construction	123	11.3 ± 1.3
Touch score	123	11.7 ± 2.8
Digit span	123	13.2 ± 2.4
Rey‐Osterrieth complex figure test
Immediate memory	120	25.0 ± 6.3
Long‐term memory	120	25.2 ± 6.5

Note: Data are presented as mean ± SD (except for sex).

Abbreviations: F, female; FEV1, forced expired volume in 1 s; FVC, forced vital capacity; M, male; SpO2, percutaneous arterial oxygen saturation; VC, vital capacity.

HA adaptive SNPs

Out of the 53 SNPs, 45 showed significant differences in the allele frequency distribution between the two population samples (p < .05, FDR corrected). Among the 45 HA adaptive SNPs, 31 were located in EGLN1, 12 in EPAS1, and 2 in PPARA (Table 3 and Figure 1).

Table 3

HA adaptive genetic variants

SNPs	CHR	Gene	BP	Minor allele/major allele	MAF		HA_allele	χ 2	P_FDR
					TIB	HAN
rs6756667	2	EPAS1	46,352,270	G/A	0.2015	0.9231	A	185.2	1.88E‑40
rs4953388	2	EPAS1	46,486,062	G/A	0.2463	0.9615	A	179.4	1.74E‑39
rs13003074	2	EPAS1	46,444,281	T/A	0.2463	0.9385	A	168.3	3.16E‑37
rs4953396	2	EPAS1	46,509,655	C/A	0.2687	0.9538	A	164.4	1.65E‑36
rs4953354	2	EPAS1	46,348,249	A/G	0.188	0.8538	G	160.5	9.61E‑36
rs1447563	2	EPAS1	46,410,225	A/C	0.2463	0.9154	C	157.5	3.50E‑35
rs7583392	2	EPAS1	46,376,299	A/G	0.1903	0.8462	G	156.3	5.57E‑35
rs896210	2	EPAS1	46,519,606	A/G	0.2368	0.8538	G	135	2.28E‑30
rs1868092	2	EPAS1	46,387,063	G/A	0.3955	1.0000	A	132.5	6.89E‑30
rs13419896	2	EPAS1	46,329,206	G/A	0.1642	0.7231	A	120.7	2.38E‑27
rs4953353	2	EPAS1	46,340,137	T/G	0.1604	0.7154	G	119.8	3.36E‑27
rs6735530	2	EPAS1	46,621,307	T/C	0.3694	0.9308	C	111.6	1.98E‑25
rs2275279	1	EGLN1	231,591,348	A/T	0.2836	0.7422	T	74.26	2.85E‑17
rs12563076	1	EGLN1	231,191,920	T/A	0.2744	0.7231	A	72.07	8.00E‑17
rs11805033	1	EGLN1	231,410,706	T/C	0.2164	0.6385	C	68.16	5.43E‑16
rs11122250	1	EGLN1	231,199,851	G/C	0.2649	0.6462	C	53.5	8.71E‑13
rs12063614	1	EGLN1	231,188,795	G/A	0.2632	0.6385	A	51.88	1.88E‑12
rs2024878	1	EGLN1	231,392,669	A/C	0.2143	0.5692	C	49.7	5.16E‑12
rs1538664	1	EGLN1	231,419,549	T/C	0.209	0.5615	C	49.68	5.16E‑12
rs2244986	1	EGLN1	231,401,706	T/C	0.1917	0.5385	C	49.48	5.42E‑12
rs2437150	1	EGLN1	231,352,778	T/C	0.2052	0.5538	C	48.97	6.70E‑12
rs973254	1	EGLN1	231,385,403	A/T	0.2068	0.5538	T	48.29	9.04E‑12
rs2749710	1	EGLN1	231,399,366	T/G	0.209	0.5538	G	47.73	1.04E‑11
rs2486746	1	EGLN1	231,361,578	A/G	0.2164	0.5615	G	47.22	1.04E‑11
rs2066140	1	EGLN1	231,368,565	C/G	0.2164	0.5615	G	47.22	1.04E‑11
rs2486745	1	EGLN1	231,369,971	C/T	0.2164	0.5615	T	47.22	1.04E‑11
rs2491405	1	EGLN1	231,377,685	A/C	0.2164	0.5615	C	47.22	1.04E‑11
rs2486741	1	EGLN1	231,380,704	C/G	0.2164	0.5615	G	47.22	1.04E‑11
rs2491407	1	EGLN1	231,383,182	A/G	0.2164	0.5615	G	47.22	1.04E‑11
rs2739506	1	EGLN1	231,393,437	G/T	0.2164	0.5615	T	47.22	1.04E‑11
rs2486734	1	EGLN1	231,395,646	G/A	0.2164	0.5615	A	47.22	1.04E‑11
rs2486731	1	EGLN1	231,396,765	G/A	0.2164	0.5615	A	47.22	1.04E‑11
rs2790882	1	EGLN1	231,397,037	T/C	0.2164	0.5615	C	47.22	1.04E‑11
rs2808584	1	EGLN1	231,377,179	G/A	0.2164	0.5625	A	47.05	1.07E‑11
rs480902	1	EGLN1	231,395,881	C/T	0.2164	0.5625	T	47.05	1.07E‑11
rs2486729	1	EGLN1	231,399,838	T/C	0.1992	0.5385	C	46.92	1.11E‑11
rs1765805	1	EGLN1	231,462,108	T/A	0.1045	0.3769	A	41.64	1.60E‑10
rs12030600	1	EGLN1	231,469,633	G/A	0.25	0.5154	A	27.61	2.11E‑07
rs16854388	1	EGLN1	231,472,514	C/A	0.2463	0.4769	A	21.35	5.31E‑06
rs2790859	1	EGLN1	231,465,611	C/T	0.2388	0.4609	T	19.96	1.07E‑05
rs961154	1	EGLN1	231,462,787	C/T	0.2351	0.4538	T	19.71	1.19E‑05
rs12093061	1	EGLN1	231,433,945	C/T	0.2388	0.4538	T	18.96	1.72E‑05
rs1769795	1	EGLN1	231,466,456	T/C	0.2388	0.4524	C	18.37	2.28E‑05
rs4253778	22	PPARA	46,234,737	G/C	0.9216	1.0000	C	10.75	0.001249
rs4253712	22	PPARA	46,199,138	A/G	0.7910	0.9077	G	8.408	0.004386
rs1769813	1	EGLN1	231,445,839	T/A	0.007463	0.02308	A	1.72	0.2005
rs1765811	1	EGLN1	231,456,224	A/G	0.007463	0.02308	G	1.72	0.2005
rs2355865	1	EGLN1	231,458,051	G/T	0.007463	0.02308	T	1.72	0.2005
rs1614148	1	EGLN1	231,462,218	A/C	0.007463	0.02308	C	1.72	0.2005
rs1769792	1	EGLN1	231,462,872	G/A	0.007463	0.02308	A	1.72	0.2005
rs6520015	22	PPARA	46,067,551	T/C	0.7239	0.7846	C	1.694	0.2005
rs7292407	22	PPARA	46,057,832	C/A	0.7612	0.8047	A	0.9405	0.3384
rs9627403	22	PPARA	46,052,596	A/G	0.7424	0.7692	G	0.335	0.5627

Note: HA adaptive SNPs were identified via χ 2 tests after correcting for multiple comparisons with the FDR_BH procedure (bold in the last column).

Abbreviations: BP, base‐pair position; CHR, chromosome; HA_allele, alleles with higher frequencies in the Tibetan highlanders than Han Chinese lowlanders; MAF (HAN), frequency of minor allele in Han samples; MAF (TIB), frequency of minor allele in Tibetan samples; SNP, single‐nucleotide polymorphism.

Figure 1

Identification of the High‐altitude (HA) adaptive allele. A χ 2 test was used to examine the allele frequency distribution difference of all candidate single‐nucleotide polymorphisms (SNPs) between the two population samples. Forty‐five SNPs showing a significant between‐group difference (p < .05, FDR corrected), the allele that exhibited a higher frequency in the Tibetans was identified HA adaptive allele. (The candidate 53 SNPs meet sample detection rate > 90%, genotype call rate > 90%, Hardy–Weinberg equilibrium p > .001, and minor allele frequency > 0.05.)

Genetic modulation of individual HA adaptive SNPs

Physiological characteristics

A total of 25 SNPs were detected that significantly modulated specific physiological characteristics in Tibetan participants (p < .05, FDR; Table 4 and Figure 2). In detail, compared with the HA nonadaptive group, the HA adaptive group exhibited a lower erythrocyte level for three SNPs in EPAS1 (rs6756667, rs7583392, and rs4953354), higher FEV1/FVC for one SNP in PPARA (rs4253712), higher SpO2 for 20 SNPs in EGLN1 (rs11805033, rs1538664, rs2024878, rs2066140, rs2244986, rs2486741, rs2486745, rs2486746, rs2491405, rs2739506, rs2790882, rs973254, rs2437150, rs2486729, rs2486731, rs2486734, rs2749710, rs2808584, rs480902, and rs2491407), and higher heart rate for one SNP in EGLN1 (rs12563076).

Table 4

Summary of genetic modulation of individual HA adaptive SNPs

SNPs	Gene	Group (adaptive/nonadaptive)	Erythrocyte	FEV1/FVC	SpO2	Heart rate	Nodal degree (left OFCant)
rs6756667	EPAS1	AA/GA + GG	8.00E‑05 (↓)
rs7583392	EPAS1	GG/GA + AA	3.70E‑04 (↓)
rs4953354	EPAS1	GG/GA + AA	2.30E‑04 (↓)
rs4253712	PPARA	GG + AG/AA		8.00E‑04 (↑)
rs11805033	EGLN1	CC/CT + TT			3.70E‑04 (↑)		6.00E‑05 (↓)
rs1538664	EGLN1	CC/TC + TT			2.80E‑04 (↑)
rs2024878	EGLN1	CC/CA + AA			3.40E‑04 (↑)		4.00E‑05 (↓)
rs2066140	EGLN1	GG/CG + CC			3.00E‑04 (↑)		1.00E‑04 (↓)
rs2244986	EGLN1	CC/TC + TT			5.00E‑05 (↑)
rs2486741	EGLN1	GG/CG + CC			2.50E‑04 (↑)		1.10E‑04 (↓)
rs2486745	EGLN1	TT/TC + CC			2.30E‑04 (↑)		1.00E‑04 (↓)
rs2486746	EGLN1	GG/AG + AA			3.10E‑04 (↑)		8.00E‑05 (↓)
rs2491405	EGLN1	CC/CA + AA			2.50E‑04 (↑)		8.00E‑05 (↓)
rs2739506	EGLN1	TT/GT + GG			3.20E‑04 (↑)		7.00E‑05 (↓)
rs2790882	EGLN1	CC/TC + TT			1.80E‑04 (↑)		7.00E‑05 (↓)
rs973254	EGLN1	TT/TA + AA			3.30E‑04 (↑)
rs2437150	EGLN1	CC/CT + TT			3.30E‑04 (↑)
rs2486729	EGLN1	CC/TC + TT			1.10E‑04 (↑)
rs2486731	EGLN1	AA/AG + GG			3.60E‑04 (↑)		7.00E‑05 (↓)
rs2486734	EGLN1	AA/AG + GG			3.10E‑04 (↑)		1.40E‑04 (↓)
rs2749710	EGLN1	GG/GT + TT			2.90E‑04 (↑)
rs2808584	EGLN1	AA/GA + GG			2.60E‑04 (↑)		6.00E‑05 (↓)
rs480902	EGLN1	TT/CT + CC			2.00E‑04 (↑)		4.00E‑05 (↓)
rs2491407	EGLN1	GG/GA + AA			2.50E‑04 (↑)		8.00E‑05 (↓)
rs12563076	EGLN1	AA/TA + TT				2.30E‑03 (↑)

Note: Data denote significance level (i.e., p values) derived from nonparametric permutation tests.

Abbreviations: FVC, forced vital capacity; FEV1, forced expiratory volume in the first second; OFCant, the anterior orbital gyrus; SNPs, single‐nucleotide polymorphisms; SpO2, percutaneous arterial oxygen saturation; ↑, HA adaptive group > HA nonadaptive group; ↓, HA adaptive group < HA nonadaptive group.

Figure 2

Violin plots showing the effects of single single‐nucleotide polymorphism on physiological variables in Tibetan participants. ** p < .01; *** p < .001

Neuropsychological variables

No significant effects were observed for any SNP on any neuropsychological variable in Tibetan participants (p > .05, FDR corrected).

Structural network attributes

No significant effects were observed for any SNP on any global or intermediate metric from any type of structural network in Tibetan participants (p > .05, FDR corrected). However, for local nodal centrality metrics from the FN‐weighted structural networks, 14 SNPs were determined to significantly modulate nodal degree of the left anterior orbital gyrus in Tibetan participants, as characterized by lower values in the HA adaptive group than in the HA nonadaptive group (p < .05, FDR corrected; Table 4 and Figure 3). All the SNPs were located in EGLN1, including rs11805033, rs2024878, rs2066140, rs2486741, rs2486745, rs2486746, rs2491405, rs2739506, rs2790882, rs2486731, rs2486734, rs2808584, rs480902, and rs2491407.

Figure 3

Violin plots showing the effects of single single‐nucleotide polymorphism on nodal centrality from the FN‐weighted structural networks in Tibetan participants. k , nodal degree; *** p < .001

Functional network attributes

No significant effects were observed for any SNP on any functional network attribute (global, intermediate, or local) in Tibetan participants (p > .05, FDR corrected).

Genetic modulation of all HA adaptive SNPs

No significant correlations were observed between the GAS over all 45 HA adaptive SNPs and any physiological characteristic in Tibetan participants (p > .05, FDR corrected).

No significant correlations were observed between the GAS over all 45 HA adaptive SNPs and any neuropsychological variable in Tibetan participants (p > .05, FDR corrected).

At the global level, the GAS over all 45 HA adaptive SNPs was negatively correlated with the clustering coefficients of the FN‐weighted structural networks in Tibetan participants (r = −.227, p = .011; Figure 4a). The negative correlations held true when the clustering coefficients were calculated from individual FA‐weighted (r = −.179, p = .046) and FL‐weighted (r = −.211, p = .019) structural networks (Figure 4a). No significant correlations were observed for the other global metrics (p > .05).

Figure 4

Relationship between the GAS and global network measures from structural networks in Tibetan participants. (a) Scatter plots showing the correlations for the GAS over all HA adaptive single‐nucleotide polymorphisms (SNPs). (b) Scatter plots showing the correlations for the GAS over HA adaptive SNPs in EGLN1. * p < .05; ** p < .01

At the intermediate level, no significant correlations were observed between the GAS over all 45 HA adaptive SNPs and any module‐related metric from the FN‐weighted structural networks in Tibetan participants (p > .05, FDR corrected). However, analyses of the FL‐weighted structural networks revealed a significantly positive correlation for the GAS over all 45 HA adaptive SNPs with intermodule communication scores between the right frontal and parietal module and the left occipital module (r = .283, p < .05, FDR corrected; Figure 5).

Figure 5

Relationship between the GAS and module‐related measures from structural networks in Tibetan participants. (a) Cortical surfaces showing modular organization of the structural brain networks. Six modules were identified for the group‐level structural brain network with a maximum Q = 0.641, which was significantly larger than those from matched random networks (n = 1,000; z = 12.180, p < .001). (b) Scatter plot showing the correlation between the GAS and intermodule communication scores. A significantly positive correlation was found between the GAS over all 45 HA adaptive single‐nucleotide polymorphisms and intermodule communication scores between the right frontal and parietal module and the left occipital module from the FL‐weighted structural networks. *** p < .001

At the local level, no regions showed significant correlations with the GAS over all 45 HA adaptive SNPs regardless of the centrality measures and types of structural networks (p > .05, FDR corrected).

At the connectional level, no connected components exhibited significant correlations with the GAS over all 45 HA adaptive SNPs in Tibetan participants regardless of the types of structural networks (p > .05, corrected).

No significant correlations were observed between the GAS over all 45 HA adaptive SNPs and any functional network attribute (global, intermediate, local, or connectional) in Tibetan participants (p > .05, corrected).

Genetic modulation of HA adaptive SNPs in different gene loci

Regardless of the gene in which the GAS was calculated over the HA adaptive SNPs, no significant correlations with any physiological characteristic were observed in Tibetan participants (p > .05, FDR corrected).

Regardless of the gene in which the GAS was calculated over the HA adaptive SNPs, no significant correlations with any neuropsychological variable were observed in Tibetan participants (p > .05, FDR corrected).

At the global level, the GAS over the 31 HA adaptive SNPs in EGLN1 was negatively correlated with the clustering coefficient from the FN‐weighted structural networks (r = −.247, p = .006; Figure 4b). The negative correlation held true when the clustering coefficients were calculated from individual FA‐weighted (r = −.211, p = .018) and FL‐weighted (r = −.254, p = .005) structural networks (Figure 4b). In addition, the GAS was negatively correlated with the local efficiency from the FA‐weighted structural networks (r = −.207, p = .021; Figure 4b). No significant correlations were observed for the GAS over either the 12 HA adaptive SNPs in EPAS1 or the 2 HA adaptive SNPs in PPARA with any global metric from any type of structural network (p > .05).

At the intermediate level, regardless of the gene in which the GAS was calculated over the HA adaptive SNPs, no significant correlations with any module‐related metric were observed regardless of the types of structural networks (p > .05, FDR corrected).

At the local level, the GAS over the 31 HA adaptive SNPs in EGLN1 was significantly correlated with the nodal degree (r = −.318) and nodal pagerank (r = −.313) of the left anterior orbital gyrus from the FN‐weighted structural networks (p < .05, FDR corrected). The GAS over the 12 HA adaptive SNPs in EPAS1 was significantly correlated with the nodal betweenness of the right anterior orbital gyrus (r = .358) from the FL‐weighted structural networks (p < .05, FDR corrected). The GAS over the two HA adaptive SNPs in PPARA was significantly correlated with the nodal eigenvector of the right middle cingulate and paracingulate gyri (r = .299), the right middle frontal gyrus (r = −.323), and the triangular part of the right inferior frontal gyrus (r = −.309) from the FN‐weighted structural networks (p < .05, FDR corrected; Figure 6).

Figure 6

Scatter plots showing the correlations between the GAS and nodal centrality from the structural networks in Tibetan participants. *** p < .001

At the connectional level, a subnetwork containing 33 regions linked by 34 connections was observed for the FN‐weighted structural networks to show a significantly negative correlation with the GAS over the 2 HA adaptive SNPs in PPARA in Tibetan participants (p = .013, corrected; Figure 7a). When superimposed on the modular architecture (Figure 5a), the 33 regions were mainly located in the right temporal and occipital module (15, 45.5%) and left temporal and parietal module (7, 21.2%), and the 34 connections included 19 (55.9%) intramodule edges predominantly linking regions in the right temporal and occipital module (14/19, 73.7%) and 15 (44.1%) intermodule edges linking regions in different modules. Given the heterogeneous sizes (number of nodes) among modules, we further tested whether specific modules were susceptible to the subnetwork. First, we randomly selected 33 regions from all 94 regions and then recorded the number of regions belonging to each module. This process was implemented 10,000 times to generate an empirical distribution for each module; according to this approach, the right temporal and occipital module was selectively involved in the subnetwork (Figure 7b). No connected components were observed for either the FA‐weighted or FL‐weighed structural networks (p > .05, corrected). No connected components showed significant correlations with the GAS over the HA adaptive SNPs in either EGLN1 or EPAS1 in Tibetan participants regardless of the types of structural networks (p > .05, corrected).

Figure 7

Relationship between the GAS and interregional connectivity from the FN‐weighted structural networks in Tibetan participants. (a) Cortical surfaces showing the subnetwork that was related to the GAS over HA adaptive single‐nucleotide polymorphisms in PPARA. The subnetwork contained 33 regions linked by 34 connections, with the regions mainly located in the right temporal and occipital module (15, 45.5%) and left temporal and parietal module (7, 21.2%), and the connections predominantly linking regions in different modules 15 (44.1%) and regions in the right temporal and occipital module (14, 41.2%). (b) Empirical frequency distributions of 33 randomly selected regions (10,000 times) in each of the six modules as depicted in Figure 7a. The right temporal and occipital module was selectively involved in the subnetwork (p < .001). The colors are the same as those in Figure 7a. For regional abbreviations, see Table 1

Regardless of the gene in which the GAS was calculated over the HA adaptive SNPs, no significant correlations were observed with any functional network attribute (global, intermediate, local, or connectional) in Tibetan participants (p > .05, FDR corrected).

Physiology–neuropsychology–brain relationships

Physiological–neuropsychological relationships

No significant correlations were observed between any physiological and neuropsychological variables in Tibetan participants (p > .05, FDR corrected).

Physiological–brain network relationships

Significant negative correlations were observed between diastolic blood pressure and nodal degree (r = −.435) and pagerank (r = −.469) of the right rolandic operculum for the FN‐weighted structural networks (p < .05, FDR corrected; Figure 8).

Figure 8

Scatter plots showing the correlations between nodal centrality from the FN‐weighted structural networks and physiological variables in Tibetan participants. ROL, the rolandic operculum; k , nodal degree; p , nodal pagerank; *** p < .001

Neuropsychological–brain network relationships

No significant correlations were observed between any neuropsychological variable and brain network measure in Tibetan participants (p > .05, FDR corrected).

Structural–functional brain network relationships

No significant correlations were observed for any global or nodal metric between structural (FN‐weighted, FA‐weighted, or FL‐weighted) and functional networks in Tibetan participants (p > .05, FDR corrected).

Relationships among different types of structural networks

Very strong positive correlations were observed for global metrics among the three types of structural networks in Tibetan participants (FN‐weighted and FA‐weighted: r = .698 ± .161; FN‐weighted and FL‐weighted: r = .773 ± .136; and FA‐weighted and FL‐weighted: r = .874 ± .063). At the local level, positive correlations that depended on the region locations and centrality metrics were observed (Figure 9).

Figure 9

Cortical surfaces showing the correlations of nodal centrality among different types of structural networks in Tibetan participants. Positive correlations that depended on region locations and centrality metrics were observed. k , nodal degree; e , nodal efficiency; b , nodal betweenness; v , nodal eigenvector; p , nodal pagerank

Reproducibility

In general, the main results of functional networks were largely reproducible when different strategies were employed to control head motion, to deal with global signals, and to define network nodes. That is, no SNP was identified to significantly modulate any global or nodal network measure (p > .05, FDR corrected), and no significant correlations were found between the GAS over all 45 HA adaptive SNPs and any global or nodal network measure (p > .05, FDR corrected). However, when the resting‐state functional MRI data did not undergo global signal regression, significantly positive correlations were observed between the GAS over the 12 HA adaptive SNPs in EPAS1 and the nodal eigenvector (r = .322) and pagerank (r = .324) of the left posterior orbital gyrus (p < .05, FDR corrected), and when the functionally defined atlas was used for network definition, a significantly positive correlation was observed between the GAS over the 12 HA adaptive SNPs in EPAS1 and the nodal eigenvector of a left parietal region (MNI coordinate = [−47 –18 50]; r = .333, p < .05, FDR corrected).

DISCUSSION

In this pilot study, we investigated the modulation of HA adaptive genetic variants on a variety of phenotypic traits, particularly large‐scale brain network organization, in indigenous Tibetan highlanders. The results showed that HA adaptive genetic variants confer more optimal and efficient cardio‐pulmonary function to ensure higher oxygen saturation levels even at lower erythrocyte levels for better adaption to HA environments. Moreover, the HA adaptive genetic variants regulated the topological organization of large‐scale structural brain networks in a gene locus, weight type, or topological scale‐dependent manner. Together, these findings provide novel insights into the genetic substrate and dynamic evolution of HA adaptive phenotypes in native Tibetans.

Modulation of HA adaptive genetic variants on physiological traits

Oxygen saturation is crucial for cellular metabolic processes and is thus considered an important indicator of human adaption to HA hypoxia. Numerous studies have reported that Tibetan highlanders exhibit higher levels of oxygen saturation than lowlanders do both at rest and during submaximal and/or maximal exercise (Gilbert‐Kawai et al., 2014). Moreover, higher levels of oxygen saturation are even found for Tibetan infants (Niermeyer, Andrade, Vargas, & Moore, 2015) and when Tibetans reascend to HA after a long period of stay at low altitude (Gonggalanzi, Bjertness, Wu, Stigum, & Nafstad, 2017). These findings suggest that higher levels of oxygen saturation have evolved into an intrinsic trait in Tibetan highlanders (Beall et al., 1997; Beall, Blangero, Williams‐Blangero, & Goldstein, 1994), presumably for better adaption to the HA hypoxic environment. In line with this speculation, our results showed that the HA adaptive group exhibited higher SpO2 values than the nonadaptive group did in terms of 20 SNPs in EGLN1. This finding is consistent with a recent study showing that two SNPs around the EGLN1 regulated SpO2 responses under acute hypoxia exposure (Yasukochi, Nishimura, Motoi, & Watanuki, 2018). Oxygen saturation is complex and depends on functional interactions among respiratory, cardiovascular, and hematological systems. Previous evidence has shown that certain EGLN1 variants affect hemoglobin–oxygen affinity (Oliveira et al., 2018), ventilatory and pulmonary vascular responses to hypoxia (Talbot et al., 2017), and hemoglobin concentration (Bhandari et al., 2017; Simonson, Huff, Witherspoon, Prchal, & Jorde, 2015). Thus, all these processes may contribute to the higher levels of oxygen saturation observed in Tibetan highlanders, although how they are directly related to oxygen saturation is not fully understood (Eaton, Skelton, & Berger, 1974; Storz, 2016). Further studies are required to unveil the common and specific contributions of these physiological processes to the modulation of HA adaptive EGLN1 variants on oxygen saturation.

Previous studies have shown that the HA environment induces polycythemia, which further elevates the risk of heart attack, stroke, and fetal loss (Semenza, 2012). However, Tibetan highlanders are protected from this risk due to a blunted erythropoietic response to the HA environment (Beall & Goldstein, 1987; Petousi & Robbins, 2014). Accumulating evidence indicates that the beneficial adaption results from HA adaptive genetic variants in Tibetan highlanders. Specifically, HA adaptive genetic variants in EPAS1 are frequently related to low erythrocyte concentrations (Beall et al., 2010; Peng et al., 2017). Consistent with these findings, we further showed that the HA adaptive group exhibited lower erythrocyte levels than the nonadaptive group did in terms of three SNPs in EPAS1 in indigenous Tibetan highlanders. It suggests that HA adaptive EPAS1 variants are a key regulatory factor involved in the low erythrocyte concentration in Tibetan highlanders, although the complete regulation pathways remain to be determined. HA adaptive response in erythropoiesis is likely a complicated process in which multiple genes are implicated. Further studies can comprehensively examine this issue by employing more sophisticated experimental designs and analysis models. It should be emphasized that although hemoglobin and erythrocytes are usually proportional, we found HA adaptive changes in the erythrocyte concentration but not the hemoglobin level in Tibetan highlanders. The separation implies that Tibetans have evolved into a better tradeoff between the oxygen diffusion efficiency and risk of blood stickiness under an HA hypoxic environment by maintaining sufficient hemoglobin levels on the one hand and lowering erythrocyte concentrations on the other hand.

The respiratory system is an organ system that directly interacts with HA hypoxia. Previous studies have consistently found that FEV1 and FVC values are higher in Tibetan highlanders than in Han Chinese born and raised at HA, mainly due to an accelerated pattern of lung growth (Weitz, Garruto, & Chin, 2016). Here, we further showed that the HA adaptive group had a higher FEV1/FVC ratio than the nonadaptive group did in terms of an SNP in PPARA in indigenous Tibetan highlanders. The higher ratio in the HA adaptive group suggests that Tibetan highlanders have evolved more efficient respiratory patterns for better adaption to the HA hypoxic environment, which ensures sustained high levels of ventilation while reducing the work of breathing. The efficient respiratory patterns may be a consequence of a higher growth of airways relative to lung size in Tibetan highlanders (Horner et al., 2017) since the FEV1/FVC ratio quantifies the degree of airflow limitation. According to previous studies, the FEV1/FVC ratio has high heritability, as revealed by both SNP‐based and pedigree‐based methods (Klimentidis et al., 2013; Wilk et al., 2009). Although genetic loci associated with FEV1/FVC are well documented (Hancock et al., 2010), the genetic basis underlying HA adaptive alterations in the respiratory system has not yet been established. Our results suggest that genetic variants in PPARA may be involved in this process. Evidence from Sherpas further shows that lower muscle PPARA expression is associated with enhanced efficiency of oxygen utilization and improved muscle energetics (Horscroft et al., 2017; Murray et al., 2018). Based on these findings, a possible regulatory pathway underlying our finding is that HA‐selective PPARA variants affect pulmonary function and thus respiratory patterns by modulating oxygen metabolism in the respiratory muscle.

Previous studies have indicated that Tibetan highlanders have evolved cardiac adaptive phenotypes different from those of lowlanders, including reduced rates of pathological right ventricular hypertrophy, increased maximum heart rate and cardiac output, and increased myocardial glucose uptake (Gilbert‐Kawai et al., 2014). Here, the HA adaptive group exhibited a higher heart rate than the nonadaptive group did among indigenous Tibetan highlanders. This finding is in accordance with a previous study showing that heart rate was associated with EGLN1 variants in Han Chinese patients with mountain sickness (Buroker et al., 2012). EGLN1 encodes prolyl hydroxylase 2, a principal negative regulator of HIFs. Based on previous studies, hypoxia could lead to the inactivation of EGLN1 and thereby increased HIF activity, which regulates gene expression involved in glucose metabolism (Aggarwal et al., 2010; Lanikova et al., 2017). Furthermore, metabolic adaptation has been associated with the preservation of cardiac function in a chronic hypoxic environment (Holden et al., 1995; Murray, 2016). Combining our results, we speculate that cardiac metabolism may be an important biological and molecular process that links HA adaptive EGLN1 variants and cardiac phenotypes.

Modulation of HA adaptive genetic variants on brain network organization

Recently, joint analyses of genes and brain imaging data have greatly enriched and deepened our understanding of the genetic bases of brain organization by linking individual SNPs with brain imaging phenotypes (Elliott et al., 2018). Here, we found limited effects of HA adaptive genetic variants of a single SNP on whole‐brain network topology in indigenous Tibetan highlanders. This result may be due to a relatively small sample size or reflect a lack of one‐to‐one mapping between them. The latter seems plausible intuitively since network measures characterize the roles of individual regions as part of an interconnected entity, and was supported by our analyses of the GAS that summarized the overall degree of genetic adaption to HA. Specifically, the GAS was correlated with the topological organization of large‐scale brain networks at multiple topological levels in indigenous Tibetan highlanders. These findings suggest that HA adaptive genetic variants may as a whole exert influence on large‐scale brain network organization. Notably, genetic modulations were observed for structural but not for functional brain networks. A possible explanation of this separation is that functional brain networks are dynamic and flexible and are thus more susceptible to environmental and environmental factors than structural brain networks. This characteristic implies the need for greater statistical power and better experimental design to detect the effects of HA adaptive genetic variants on functional brain networks.

We found a negative correlation between the GAS and the clustering coefficient of structural brain networks at the global level, a finding that was independent of network weight types. Clustering coefficient measures, on average, the extent of connectivity among a node's first‐degree neighbors in a network and reflects the ability of the network to support segregated information processing. Thus, the observed negative correlations indicate that the more HA adaptive alleles Tibetan highlanders carry, the less local interconnectivity their structural brain networks have. HA exposure is reportedly associated with modifications in white matter microstructure, such as decreased FA (Paola et al., 2008; Zhang, Zhang, et al., 2013). Evidence from laboratory rodents further shows that hypoxia exposure causes disrupted myelination in the brain (Kanaan, Farahani, Douglas, Lamanna, & Haddad, 2006). The microstructural degeneration of white matter can reduce the fidelity of communication and conduction velocity between brain regions (Bartzokis, Lu, & Mintz, 2007; Salami, Itami, Tsumoto, & Kimura, 2003). Interestingly, the observed negative correlation between the GAS and clustering coefficient was mainly driven by HA adaptive genetic variants in EGLN1. Substantial evidence from indigenous Tibetans has shown that HA adaptive EGLN1 variants are related to lower hemoglobin concentrations, which is pivotal for brain oxygen consumption and cellular energy supply (Beall et al., 2010; Lorenzo et al., 2014). Furthermore, local loss of cerebral blood flow, such as that caused by a stroke, reportedly leads to damage to oligodendrocytes and demyelination (Tekkok & Goldberg, 2001), which further results in disruptions in axonal integrity of the lesioned area (Saab & Nave, 2016). Accordingly, a plausible hypothesis regarding the observed negative correlation is that HA adaptive EGLN1 variants reduce hemoglobin concentrations in Tibetans, which leads to insufficient energy supply to the brain, followed by further impairments in the microstructure of white matter and ultimately, an adaptive decrease in network‐level communication. Notably, the mechanisms underlying any step in this speculative process, such as damage to myelin owing to oxygen and glucose deficits, may be far more complex than anticipated (Hamilton, Kolodziejczyk, Kougioumtzidou, & Attwell, 2016). Future studies that combine neuroimaging with cellular and molecular biology could help clarify this speculation.

At the intermediate level, we found a positive correlation between the GAS and FL between the left occipital module and right frontal and parietal module in indigenous Tibetan highlanders. These two modules covered most regions of the dorsal visual pathway that connects the primary visual cortex to the posterior parietal cortex. Functionally, the dorsal visual streams are involved in visual processing of spatial localization and visuomotor control (Galletti & Fattori, 2018), visuospatial construction (Morris, 2010; Sarpal et al., 2008), visual memory (Dimond, Perry, Iaria, & Bray, 2019; Hirel et al., 2017), and attention (Atkinson & Braddick, 2011; Tonks et al., 2019). In addition, several frontal regions have been implicated in visual working memory, visual search, and visual attention, such as the right inferior triangular frontal gyrus and frontal eye field (a region around the intersection of the right middle frontal gyrus with the right precentral gyrus), both of which were included in the right frontal and parietal module in this study (Offen, Gardner, Schluppeck, & Heeger, 2010; van Ewijk et al., 2015). Indeed, previous findings have suggested the existence of a right fronto‐parietal network for sustained attention in favor of rapid visual information processing (Coull, Frith, Frackowiak, & Grasby, 1996). Regarding the interaction between occipital regions and the fronto‐parietal network, a previous study found that transcranial magnetic stimulation to the right frontal eye field can affect activity in the left visual cortex (Marshall, O'Shea, Jensen, & Bergmann, 2015). Given accumulating evidence that the long‐distance effect on remote sites by transcranial magnetic stimulation is through the modulation of targeted functional networks (Wang, Rogers, et al., 2014), the white matter tract bundles between the left occipital module and right frontal and parietal module may thus be a structural basis underlying visual information flow for both top‐down control and bottom‐up feedback between the two modules because functional connectivity patterns are largely shaped by structural layouts (Greicius, Supekar, Menon, & Dougherty, 2009; van den Heuvel, Mandl, Kahn, & Hulshoff Pol, 2009). Returning to our findings, indigenous Tibetan highlanders with larger GAS and thus better adaption to HA were associated with longer white matter tract bundles between these two modules. This result is comparable to a previous study showing that long‐term acclimatization to HA hypoxia can lead to lengthened commissural fibers connecting the bilateral visual cortex (Chen, Li, et al., 2016). Currently, interpreting the HA adaptive increase in vision cortex‐related FL is difficult in the absence of a full understanding of underlying cellular and molecular events. Presumably, this adaption may result from plastic fiber rewiring or the generation of new projections due to axonal sprouting (Johansen‐Berg, 2007). Since longer fibers are related to more decay and delay in neuronal signal conduction (Hugdahl, 2011), the increased FL may imply low efficiency for visual information processing in Tibetan highlanders. Although direct evidence is lacking for this speculation, we noted that significantly lower values from the pursuit aiming tests (reflect the ability to make quick and accurate movements and is thus related to visuomotor control) were reported in participants native to HA of 3,700, 4,500, and 5,100 m than in participants at sea level and sea‐level participants acclimated for 5 days at 3,700 m (Zhang, Zhou, et al., 2013). Moreover, accumulating evidence from behavioral and neuropsychological studies has shown that HA immigrants and their descendants as well as mountain climbers exhibit poor performance in visual working memory, visual construction, and attention (Chen et al., 2017; Wang, Ma, et al., 2014). In the future, clarifying the extent to which the structural connectivity alterations of the visual cortex could account for visual deficits due to short‐ and long‐term exposure to HA hypoxia would be interesting. Notably, our previous study found that the visual cortex‐related FL increase due to long‐term acclimatization to HA hypoxia was related to higher functional connectivity in HA immigrants (Chen, Li, et al., 2016). This finding suggests that other restorative processes occurred in the brain to functionally compensate for the detrimental structural alterations with respect to signal conduction. Since human brain networks are optimized not only for minimal global wiring cost but also for maximal efficiency of information propagation, there may be an alternative interpretation: the HA adaptive FL increase is a compromise to maintain an optimal balance between specialized and integrated information processing but at a relatively low wiring cost under an HA hypoxic environment.

At the node level, several frontal regions showed effects of HA adaptive genetic variants, including the anterior orbital gyrus, middle frontal gyrus, inferior frontal triangular gyrus and middle cingulate, and paracingulate gyri. These regions show HA adaptive alterations in both local morphology and interregional connectivity (Chen et al., 2017; Wei et al., 2017). Interestingly, most of the regions belonged to the right frontal and parietal module, which plays a key role in the top‐down modulation of various visual processes as discussed above. For example, a previous study showed that reduced gray matter volume in the inferior frontal triangular gyrus was related to poor visuospatial memory in patients with chronic obstructive pulmonary disease (characterized by decreased oxygen supply; Zhang, Wang, et al., 2013). Thus, the nodal results suggest that frontal regions regulating and integrating visual information flow exhibit dynamic adaption to HA in indigenous Tibetan highlanders. Notably, different mechanisms governing the involvement of the frontal regions in HA adaptive genetic variants may exist because opposite correlation directions were observed among them. Further studies are required in the future to clarify this point.

At the connectivity level, interregional FNs within a subnetwork that was predominantly involved in regions in the right temporal and occipital module were negatively correlated with the GAS. Anatomically, this module overlaps, to a great extent, with the ventral visual pathway that courses through the occipitotemporal cortex to the anterior part of the inferior temporal gyrus. As a key region in the ventral visual streams, the middle temporal gyrus shows altered morphology in native Tibetans compared with that in lowlanders, and the alterations are positively correlated with altitude (Wei et al., 2017). Here, our finding implies a decrease in the FNs interconnecting different components of the ventral visual pathway with increasing adaption to HA. The ventral visual streams are mainly involved in the visual search process (Zhou & Desimone, 2011) and visual perception (Eldridge et al., 2018). In previous studies, compared with lowlanders, HA residents and immigrants showed worse performance in the visual search process (Chen, Li, et al., 2016) and visual perception (Zhang, Liu, Yan, & Weng, 2011). According to our findings, the visual deficits may be due to decreased interconnectivity in the ventral visual pathway under the HA hypoxia environment. Presumably, the decrease in connectivity may be attributed to a complex array of changes in astrocyte numbers and morphology and the proliferation and differentiation of oligodendrocyte precursor cells into myelinating cells (Hofstetter et al., 2019; Wake et al., 2015).

Associations among phenotypes

According to previous studies, compared with other ethnicities in China, Tibetans have a higher essential hypertension prevalence (Zhao et al., 2012). This difference might be due to environmental cold temperature influence (Aryal, Weatherall, Bhatta, & Mann, 2016), high dietary salt intake (Sehgal, Krishan, Malhotra, & Gupta, 1968), or hypertension‐related genetic variants (Guo et al., 2017; Shi et al., 2018). Here, the nodal strength of the right rolandic operculum was negatively correlated with diastolic blood pressure in indigenous Tibetan highlanders. This finding suggests a vital role played by the rolandic operculum in the higher essential hypertension prevalence in Tibetans. The rolandic operculum is involved in the processing and awareness of interoceptive signals (Blefari et al., 2017). This structure is anatomically adjacent to and functionally interacts with the insula (Maliia et al., 2018). Accumulating evidence indicates that insular cortices are involved in the regulation of autonomic nervous system activity that further affects blood pressure (Oppenheimer, Kedem, & Martin, 1996; Shoemaker, Wong, & Cechetto, 2012). Accordingly, we speculate that the modulation of the rolandic operculum on diastolic blood pressure is achieved through bidirectional exchange of information flow with the insula. In addition to the relationship between brain network organization and blood pressure, we found that different types of structural brain networks correlated with each other with respect to both global and local network architecture. Nonetheless, our results demonstrated that the effects of HA adaptive genetic variants on structural brain networks were dependent on how they were constructed. These findings indicate the mutually complementary nature of different types of structural brain networks and suggest the need to develop integrative models to summarize the unique information carried by them.

Limitations and future directions

First, with respect to the genes involved in HA adaption, we chose 60 SNPs in three genes a priori based on previous literature. However, this candidate variant study may underestimate the associations between HA adaptive genetic variants and interindividual variance in phenotypes. Thus, genome‐wide association studies are required in the future to comprehensively map phenotypic correlates of HA adaptive genetic variants. Second, given the small sample size and numerous analyses, we did not use any strategy to control overall multiple comparisons. Therefore, our findings may suffer from a high risk for false positives. Future studies are required to test the reproducibility of our fining on an independent, large cohort of participants with more rigorous correlation methods. Third, confined by the nature of the cross‐sectional design, the mapped genotype–phenotype associations should be explained with caution in the current study. In particular, the modulation of phenotypes by HA adaptive genetic variants may involve multiple cellular and molecular processes. Mechanistic insights into this question could benefit from future long‐term follow‐up studies that allow examining the sequential emergence and dynamic development of HA adaptive phenotypes by combining disciplines such as cell and molecular biology. Fourth, “best” practices to determine an optimal analytical strategy for conducting graph‐based brain network studies are still lacking. Several studies have shown that both functional and structural brain networks vary considerably in their quantitative topological descriptions across different choices of brain parcellation schemes for network node definition (Wang et al., 2009; Zalesky et al., 2010), similarity measures for interregional connectivity estimation (Sarwar, Ramamohanarao, & Zalesky, 2019), and thresholding procedures for denoising (van den Heuvel et al., 2017). Examining the extent to which the current findings are reproducible over different choices of these factors is important. Fifth, the age distribution of the participants in this study was narrow, which may reduce the possibility of observing positive results due to small interindividual variance. The extent to which the genotype–phenotype mapping is underestimated could be determined in the future by recruiting participants with large age spans. Finally, this study utilized the two most popular neuroimaging techniques (i.e., diffusion and resting‐state functional MRI) with increasingly well‐understood methodology to investigate effects of HA adaptive genetic variants on brain network origination. Recently, methodological advances allow constructing and charactering morphological and metabolic brain networks at the individual level (Tijms, Series, Willshaw, & Lawrie, 2012; Wang, Jin, Zhang, & Wang, 2016; Yao, Hu, Nan, Zheng, & Xie, 2016). It is interesting in the future to combine different types of brain networks for comprehensively mapping efforts of HA adaptive genetic variants.

CONCLUSION

In conclusion, this pilot study associated HA adaptive genetic variants with interindividual variance in a variety of phenotypic traits in indigenous Tibetan highlanders. The results showed that HA adaptive genetic variants selectively modulate physiological characteristics and large‐scale structural brain networks in a complicated, gene loci‐dependent manner. This unique genotype–phenotype mapping is helpful for explaining HA adaptive phenotypic alterations in previous literature on the one hand and has important implications for understanding dynamic human evolution in response to the HA environment on the other hand. It should be noted that these findings are from a very small sample, and thus this study should be regarded as exploratory in nature. Future studies are required to test the reproducibility of our findings in an independent sample.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

J.Z., Q.G. and J.W. designed experiments; Z.G., C.F., T.L., W.Y., L.G. carried out experiments; Z.G., N.W. and J.W. analyzed experimental data; Z.G. and J.W. wrote the manuscript.