Participants
A total of 157 healthy young adults were recruited by advertisement. All participants met the inclusion criteria of Chinese Han, right handedness, and within a restricted age range of 18-30 years. Exclusion criteria included neuropsychiatric or severe somatic disorder, a history of alcohol or drug abuse, regular smoker, current medication (e.g., antibiotics or sedative hypnotics) within a month, pregnancy, MRI contraindications, and a family history of psychiatric illness among first-degree relatives. The MINI-International Neuropsychiatric Interview (M.I.N.I.) and Alcohol Use Disorders Identification Test (AUDIT) were used in the process of excluding participants. This study was approved by the ethics committee of The First Affiliated Hospital of Anhui Medical University. Written informed consent was obtained from all participants after they had been given a complete description of the study.
Go/No-Go task
The Go/No-Go task was conducted on a computer to assess the ability of behavioral inhibition using E-Prime 2.0 (http://www.pstnet.com/eprime.cfm) [27]. During the task, the letter X or Y was presented at a frequency of 1 Hz on the screen. In “Go” conditions, the current letter is different from the previous one and participants should respond quickly by pressing the button within 900 ms. In “No-Go” conditions (10% of all trials), the current letter is the same as the previous one and participants cannot press the button; if one presses the button, it would be counted as an error. The Go/No-Go task consisted of a practice test and a formal test. There were 20 trials (15 “Go” trials and 5 “No-Go” trials) in the practice test. If a participant responds correctly in 3 “No-Go” trials, he or she can shift to the formal test; otherwise, the participant needs to restart the practice test. The formal test was divided into two groups with 210 trials in each group and 30 s break between the two groups. It took about 12 min for the Go/No-Go task. The primary variable of interest is the accuracy in “No-Go” conditions (Acc_No-Go) that reflects behavioral inhibition.
MRI data acquisition
MRI scans were obtained using a 3.0-Tesla MR system (Discovery MR750w, General Electric, Milwaukee, WI, USA) with a 24-channel head coil. Earplugs were used to reduce scanner noise, and tight but comfortable foam padding was used to minimize head motion. High-resolution 3D T1-weighted structural images were acquired by employing a brain volume (BRAVO) sequence with the following parameters: repetition time (TR) = 8.5 ms; echo time (TE) = 3.2 ms; inversion time (TI) = 450 ms; flip angle = 12°; field of view (FOV) = 256 mm × 256 mm; matrix size = 256 × 256; slice thickness = 1 mm, no gap; 188 sagittal slices. The perfusion imaging was performed using a pseudo-continuous ASL sequence with a 3D fast spin-echo acquisition and background suppression (TR = 5070 ms, TE = 11.5 ms; post-label delay = 2025 ms; spiral in readout of eight arms with 512 sample points; flip angle = 111°; FOV = 240 mm × 240 mm; reconstruction matrix = 128 × 128; slice thickness = 3 mm, no gap; 50 axial slices; number of excitation = 3). The label and control whole-brain image volumes required 8 TRs, respectively. A total of three pairs of label and control volumes were acquired. Resting-state blood-oxygen-level-dependent (BOLD) fMRI data were acquired using a gradient-echo single-shot echo planar imaging (GRE-SS-EPI) sequence with the following parameters: TR = 2000 ms; TE = 30 ms; flip angle = 90°; FOV = 220 mm × 220 mm; matrix size = 64 × 64; slice thickness = 3 mm, slice gap = 1 mm; 35 interleaved axial slices; 185 volumes. DTI data were acquired using a spin-echo single-shot echo planar imaging (SE-SS-EPI) sequence with the following parameters: TR = 10000 ms; TE = 74 ms; flip angle = 90°; FOV = 256 mm × 256 mm; matrix = 128 × 128; slice thickness = 3 mm without gap; 50 axial slices; 64 diffusion gradient directions (b = 1000 s/mm2) plus five b = 0 reference images. All images were visually inspected to ensure that only images without visible artifacts were included in subsequent analyses.
Gray matter volume analysis
Voxel-based morphometry (VBM) analysis was performed using the CAT12 toolbox (http://www.neuro.uni-jena.de/cat) implemented in the Statistical Parametric Mapping software (SPM12, http://www.fil.ion.ucl.ac.uk/spm). First, all the structural T1-weighted images were corrected for bias-field inhomogeneities. Second, these images were segmented into gray matter, white matter and cerebrospinal fluid density maps using the “new-segment” approach [62]. Third, a diffeomorphic anatomical registration through the exponentiated Lie algebra (DARTEL) technique was used to generate a custom, study-specific template [63]. Fourth, each participant’s gray matter density image was warped to the customized template; then the resultant images were affine registered to the Montreal Neurological Institute (MNI) space and resampled to a voxel size of 1.5 mm × 1.5 mm × 1.5 mm. Fifth, the modulation was applied by multiplying the transformed gray matter density maps with the non-linear components of Jacobian determinants, which resulted in the normalized GMV maps representing the local native-space GMV after correcting the confounding effect of variance induced by individual whole-brain size. Finally, the resultant GMV images were smoothed with a 6 mm full-width at half maximum (FWHM) Gaussian kernel.
Cerebral blood flow analysis
Three ASL difference images were calculated by subtracting the label images from the control images and then averaged. Next, CBF was quantified by applying a single-compartment model [64] to the mean ASL difference and proton-density-weighted reference images [65-67]. SPM12 software was used to normalize the CBF images into the MNI space using the following steps: (1) individual structural images were firstly co-registered with the CBF images; (2) the transformed structural images were segmented and normalized to the MNI space; and (3) the CBF image of each subject was written into the MNI space using the deformation parameter derived from the prior step and was resliced into a 2-mm cubic voxel. For the purpose of standardization, the CBF value of each voxel was divided by the global mean CBF value. Finally, the CBF images were smoothed with a 6 mm FWHM Gaussian kernel.
Functional connectivity strength analysis
Resting-state BOLD data were preprocessed using SPM12 and Data Processing & Analysis for Brain Imaging (DPABI, http://rfmri.org/dpabi) [68]. The first 10 volumes for each participant were discarded to allow the signal to reach equilibrium and the participants to adapt to the scanning noise. The remaining volumes were corrected for the acquisition time delay between slices. Then, realignment was performed to correct the motion between time points. Head motion parameters were computed by estimating the translation in each direction and the angular rotation on each axis for each volume. All participants’ BOLD data were within the defined motion thresholds (i.e., translational or rotational motion parameters less than 2 mm or 2°). We also calculated frame-wise displacement (FD), which indexes the volume-to-volume changes in head position. Several nuisance covariates (the linear drift, the estimated motion parameters based on the Friston-24 model, the spike volumes with FD > 0.5, the white matter signal, and the cerebrospinal fluid signal) were regressed out from the data. The datasets were then band-pass filtered using a frequency range of 0.01 to 0.1 Hz. In the normalization step, individual structural images were firstly co-registered with the mean functional image; then the transformed structural images were segmented and normalized to the MNI space using a high-level nonlinear warping algorithm, that is, the DARTEL technique. Finally, each filtered functional volume was spatially normalized to MNI space using the deformation parameters estimated during the above step and resampled into a 3-mm cubic voxel.
Functional connectivity strength (FCS) is a graph theory measure that evaluates functional connectivity of each voxel with all other voxels across the whole gray matter [69-71]. Firstly, we computed Pearson’s correlation coefficients between the BOLD time courses of all pairs of voxels and obtained a whole gray matter functional connectivity matrix for each participant. For a given voxel, FCS was computed as the sum of positive functional connectivity above a threshold of 0.25 between that voxel and all other voxels within the whole gray matter. Then, we normalized the FCS value of each voxel by dividing it by the global mean FCS value. Finally, the FCS maps were smoothed with a 6 mm FWHM Gaussian kernel.
White matter integrity analysis
For DTI data, standard processing steps were performed by using the FMRIB Software Library (FSL, www.fmrib.ox.ac.uk/fsl). First, eddy current distortion and head motion were corrected by registering the diffusion-weighted images to the first b0 image through the affine transformations. Second, the data were skull-stripped by using the FMRIB Brain Extraction Tool. Finally, diffusion parameters including fractional anisotropy (FA), axial diffusivity (AD), radial diffusivity (RD), and mean diffusivity (MD), were calculated by using the DTIFIT toolbox. Then, tract-based spatial statistics (TBSS) pipeline was conducted [51, 72]. Briefly, individual FA images were firstly non-linearly registered to the MNI space. After transformation into the MNI space, mean FA image was created and thinned to generate a mean FA skeleton. Then, each subject’s FA image was projected onto the skeleton via filling the mean FA skeleton with FA values from the nearest relevant tract center by searching perpendicular to the local skeleton structure for maximum FA value. Finally, the registration and projection information derived from the FA analysis was applied to the other diffusion parameters to project AD, RD, and MD images onto this common skeleton.
Fecal samples collection and gut microbiota analysis
Fecal samples were collected in sterilized tubes and stored immediately in a -80 ℃ freezer within 1 day before or after MRI examination. Microbial genome DNA was extracted from the fecal samples using a QIAamp DNA Stool Mini Kit (Qiagen Inc., Hilden, Germany). To construct the Polymerase Chain Reaction (PCR)-based 16S rDNA gene amplicon library for sequencing, PCR enrichment of the V4 hypervariable region of 16S rDNA gene was performed with the forward primer 515F (5’-GTGCCAGCMGCCGCGGTAA-3’) and reverse primer 806R (5’-GGACTACHVGGGTWTCTAAT-3’). The qualified amplicon mixture was then sequenced on the MiSeq platform with the PE250 sequencing strategy. Before the 16S rDNA data analysis, raw reads were filtered to remove adaptors and low-quality and ambiguous bases, and then paired-end reads were added to tags by the Fast Length Adjustment of Short reads program (FLASH, v1.2.11) [73]. The tags were clustered into operational taxonomic units (OTUs) with a cutoff value of 97% using UPARSE software (v9.1.13) [74] and chimera sequences were compared with the Gold database using UCHIME (v4.2.40) [75] to detect. Then, the representative sequence from each OTU cluster was obtained. These OTU representative sequences were taxonomically classified using Ribosomal Database Project (RDP) Classifier (v.2.2) [76] with a minimum confidence threshold of 0.8, and the training database was the Greengene Database (v201305) [77]. The USEARCH_global [78] was used to compare all tags back to OTU to get the OTU abundance statistics table of each sample. Alpha diversity was assessed using the species richness indices (Sobs, Chao, and Ace) and species diversity indices (Shannon and Simpson that reflect both species richness and species evenness) [79, 80], which were calculated by MOTHUR (v1.31.2) [81] and QIIME (v1.8.0) [82] at the OTU level. The species accumulation curves were plotted in supplementary Fig. S1, which indicated that the sampling amount was sufficient.
Statistical analysis
Demographic, cognitive, and gut microbial variables were compared between males and females using two sample t-tests in the SPSS 26.0 software (SPSS Inc., Chicago, IL, United States).
In the male and female groups separately, voxel-based partial correlation analyses between alpha diversity and brain imaging measures (GMV, CBF, and FCS) were performed using multiple regression analyses in the SPM12 software. For CBF analyses, age was included as a nuisance covariate, with total intracranial volume (TIV) and FD as additional covariates for GMV and FCS analyses respectively. Multiple comparisons were corrected using a cluster-level family-wise error (FWE) method, resulting in a cluster defining threshold of p = 0.001 and a corrected cluster significance of p < 0.05. For the TBSS analysis, non-parametric permutation testing (permutation number = 5000) and threshold-free cluster enhancement (TFCE) in the FSL software were used for statistical inference of the partial correlations between alpha diversity and diffusion parameters (AD, RD and MD) controlling for age. The FWE method was also used to correct for multiple comparisons with a corrected significance threshold of p < 0.05. In case of significant correlations identified for any brain regions in either males or females, these significant regions were defined as regions of interest (ROIs) and mean imaging values within these ROIs were extracted to further examine whether there were significant sex differences in the correlations. That is, ROI-based partial correlation coefficients between imaging measures and alpha diversity were transformed into Fisher’s Z scores and then compared between males and females [83]. To estimate the effect sizes, we also calculated Cohen’s q (no effect: q < 0.1, small effect: 0.1 < q < 0.3, intermediate effect: 0.3 < q < 0.5, large effect: q > 0.5) [84].
For brain imaging parameters showing correlations with gut microbial diversity, we further examined their associations with the ability of behavioral inhibition (Acc_No-Go) using partial correlation analyses. To test whether the association between variables was mediated by other variables, mediation analysis was performed using the PROCESS macro (http://www.processmacro.org/) [85]. In the mediation analysis model (Fig. 8A), all paths were reported as unstandardized ordinary least squares regression coefficients, namely, total effect of X on Y (c) = indirect effect of X on Y through M (a × b) + direct effect of X on Y (c’). The significance analysis was based on 5000 bootstrap realizations and a significant indirect effect is indicated when the bootstrap 95% confidence interval (CI) does not include zero. In this study, only variables that showed a significant correlation with others in the correlation analyses were considered independent (alpha diversity), dependent (cognition), or mediating variables (neuroimaging parameters) in the mediation analysis.