The participants included in this study were 6 stress participants submitted to prolonged psychological stress exposure (3 males, 3 females; mean age, 23.83 ± 0.37), the same 6 stress participants, 6 weeks after the end of the exposure to stress and 6 controls (3 males, 3 females; mean age, 24.33 ± 1.24). Control participants included a cohort of medical students under their normal academic activities, whereas the stress group included participants that had just finished their long period of preparation for the medical residence selection exam (stress group). Participants responded to a laterality test and to a self-administered questionnaire regarding stress assessment (Perceived Stress Scale—PSS) (Cohen et al., 1983). Participants were further assessed with the Hamilton anxiety scale—HAS (Hamilton, 1959) and the Hamilton depression scale—HDS (Hamilton, 1967) by a certified psychologist. Upon filling of the questionnaires, and immediately before the imaging acquisitions, participants collected saliva samples with the help of Salivette (Sarstedt, Germany) collection devices. Collection took place between 9 and 5 p.m. in all participants. Samples were stored at −20°C until the biologically active, free fraction of the stress hormone cortisol was analyzed using an immunoassay (IBL, Hamburg).

The study was conducted in accordance with the principles expressed in the Declaration of Helsinki and was approved by the Ethics Committee of Hospital de Braga (Portugal). The study goals and tests were explained to all participants and all gave informed written consent.

Participants were scanned on a clinical approved Siemens Magnetom Avanto 1.5 T (Siemens Medical Solutions, Erlangen, Germany) on Hospital de Braga using the Siemens 12-channel receive-only head coil. The imaging sessions, including one structural T1, one resting state functional, and two task related functional acquisitions, were conducted in the same day and the Siemens Auto Align scout protocol was used to minimize variations in head positioning. For structural analysis, a T1 high-resolution anatomical sequence, 3D MPRAGE (magnetization prepared rapid gradient echo) was performed with the following scan parameters: repetition time (TR) = 2.4 s, echo time (TE) = 3.62 ms, 160 sagittal slices with no gap, field-of-view (FoV) = 234 mm, flip angle (FA) = 8°, in-plane resolution = 1.2 × 1.2 mm² and slice thickness = 1.2 mm. During resting-state fMRI acquisition, using gradient echo T2^* weighted echo-planar images (EPIs), participants were instructed to keep the eyes closed and to think about nothing in particular. The imaging parameters were: 100 volumes, TR = 3 s, TE = 50 ms, FA = 90°, in-plane resolution = 3.4 × 3.4 mm², 30 interleaved slices, slice thickness = 5 mm, imaging matrix 64 × 64 and FoV = 220 mm. fMRI paradigm acquisition was acquired using: TR = 2 s, TE = 20 ms, FA = 90°, in-plane resolution and slice thickness 3.3 mm, 38 ascending interleaved axial slices with no gap and FoV = 212 mm. The functional paradigm acquisitions were previously described (Soares et al., 2012) and the paradigm was presented using the fully integrated fMRI system IFIS-SA.

Before any data processing and analysis, all acquisitions were visually inspected and confirmed that they were not affected by critical head motion and that participants had no brain lesions.

To achieve signal stabilization and allow participants to adjust to the scanner noise, the first 5 resting state fMRI volumes (15 s) were discarded. Data preprocessing was performed using SPM8 (Statistical Parametrical Mapping, version 8, http://www.fil.ion.ucl.ac.uk) analysis software. Images were firstly corrected for slice timing using first slice as reference and SPM8's Fourier phase shift interpolation. To correct for head motion, images were realigned to the mean image with a six-parameter rigid-body spatial transformation and estimation was performed at 0.9 quality, 4 mm separation, 5 mm FWHM smoothing kernel using 2nd degree B-Spline interpolation. No participants exceed head motion higher than 2 mm in translation or 1° in rotation. Images were then spatially normalized to the MNI (Montreal Neurological Institute) standard coordinate system using SPM8 EPI template and trilinear interpolation. Data were then re-sampled to 3 × 3 × 3 mm3 using sinc interpolation, smoothed to decrease spatial noise with a 8 mm, full-width at half-maximum (FWHM), Gaussian kernel, temporally band-pass filtered (0.01–0.08 Hz) and the linear trend was removed. The pre-processing of fMRI paradigm images was previously described (Soares et al., 2012).

Spatial independent component analysis was conducted for using the Group ICA 2.0d of fMRI Toolbox (GIFT, http://www.icatb.sourceforge.net) (Calhoun et al., 2001; Correa et al., 2005). Concisely, spatial ICA analysis is a fully data-driven approach that consists in extracting the non-overlapping spatial maps with temporally coherent time courses that maximize independence. The methodology employed by GIFT can be summarized in three main stages: dimensionality reduction, estimation of the group independent components, and back-reconstruction of each subject's corresponding independent components. The reduction of dimensionality of the functional data and computational load was performed with Principal Component Analysis (PCA) in the concatenated dataset over all subjects, independently of the groups. Then, 20 independent components were estimated, based on a good trade-off (clustering/splitting) between preserving the information in the data while reducing its size (Beckmann et al., 2005; Zuo et al., 2010), using the iterative Infomax algorithm. The ICASSO tool was used to assess the ICA reliability, and 20 computational runs were performed on the dataset, during which the components were being recomputed and compared across runs and the robustness of the results was ensured (Himberg et al., 2004). The previous steps result in the estimation of a mixing matrix with partitions, unique to each subject. The individual independent components were then back-reconstructed from the group-level components. This back-reconstruction step is accomplished by projecting each subject's data onto the inverse of the partition of the calculated matrix corresponding to that subject. The obtained independent components were expressed in t-statistic maps, which were finally converted to a z-statistic. Z-statistic describes the voxels that contributed more intensely to a specific independent component, providing a degree of functional connectivity within the network (Bartels and Zeki, 2005; Beckmann et al., 2005). The final components were visually inspected, sorted, and spatially correlated with resting state functional networks from (Shirer et al., 2012). Each subject's map corresponding to the best-fit component of each RSN was used to perform group statistical analyses.

The fMRI decision-making paradigm analyzed to investigate the task-induced deactivations consisted of two different event-related jittered design sessions. First session of valued actions with reward delivery and, after 30 min break, the second session consisted of the devalued actions with the outcome devaluation and extinction. Both sessions had 150 trials, each with 1.5 s for decision, 4 s with the choice highlighted, and 2 s for reward delivery, followed by the interstimulus interval with mean duration of 4 s [please see Soares et al. (2012), for further details].

fMRI paradigm was analyzed by creating a set of regressors at resting and decision making periods, which were convolved with the hemodynamic response function. In order to reliably map task-induced deactivations, we combined all the resting periods (resting baseline condition) and all the decision periods (decision condition), given that decision periods were equally demanding. The contrast used to assess task-induced deactivations was the resting baseline condition minus decision condition. Resulting functional patterns were masked with the previously described RSNs masks (Shirer et al., 2012).

Structural analysis based on segmentation of brain structures from T1 high-resolution anatomical data was performed using the freely available Freesurfer toolkit version 5.0 (http://surfer.nmr.mgh.harvard.edu). Intracranial volume (ICV) was used to correct the volumes and the processing pipeline was the same as previously described (Soares et al., 2012). DMN was defined by the summed volume of the angular gyrus of inferior parietal lobe, the posterior cingulate, the precuneus, and the frontopolar region (Raichle et al., 2001; Buckner et al., 2008). The summed volume of the middle frontal gyrus (dorsolateral and prefrontal region) and the posterior parietal region constituted the DAN (Seeley et al., 2007; Sridharan et al., 2008). VAN was constituted by the sum of the temporal-parietal junction and the ventral frontal cortex volumes (Fox et al., 2006). SMN was defined by the summed volume of the paracentral, precentral postcentral, and the cerebellum (Shirer et al., 2012). The summed volume of the cuneus, pericalcarine, and the lingual region constituted the primary VN (Shirer et al., 2012). Auditory network (AN) was defined by the summed volume of the temporal transverse and the temporal superior (Shirer et al., 2012).

Results of the psychological scales, cortisol levels, and regional volumes were analyzed in the IBM SPSS Statistics software, v.19 (IBM, New York). Comparisons between the stress recovered and stress were done with paired samples t-test and between stress recovered and control with two-tailed independent-samples t-test. For all these comparisons significance level was set at 0.05. Values are presented as mean ± standard error of the mean.

Group analysis of the resting state fMRI and task induced deactivations were performed using the second level random effect analyses in SPM8. Initially, within group analyses were performed only to confirm the functional connectivity of the RSNs in the different groups, using one-sample t-tests. Therefore, between group analyses were implemented with directional two-sample t-tests. Functional results for all RSNs were considered significant at p < 0.05 corrected for multiple comparisons using a combination of an uncorrected height threshold of p < 0.025 with a minimum cluster size. The cluster size was determined over 1000 Monte Carlo simulations using AlphaSim program distributed with REST software tool (http://resting-fmri.sourceforge.net/). AlphaSim input parameters were the following: individual voxel probability threshold = 0.025, cluster connection radius = 3 mm, gaussian filter width (FWHM) = 8 mm, number of Monte Carlo simulations = 1000 and mask was set to the corresponding RSN template mask. Anatomical labeling was defined by a combination of visual inspection and Anatomical Automatic Labeling atlas (AAL) (Tzourio-Mazoyer et al., 2002).

Stress impact was confirmed in several parameters: PSS [Figure 1A; t₍₁₀₎ = 2.52, P < 0.05; Stress Group M = 35.50; SD = 2.59; Control Group M = 30.17; SD = 4.49] and in the HAS [Figure 1A; t₍₁₀₎ = 2.37, P < 0.05 Stress Group M = 11.00; SD = 7.95; Control Group M = 3.00; SD = 2.28], and depression scores [HAD, Figure 1A; t₍₁₀₎ = 3.65, P < 0.01; Stress Group M = 7.50; SD = 2.59; Control Group M = 3.17; SD = 1.33], but only by a trend when it regards to salivary cortisol levels [Figure 1B; t₍₁₀₎ = 1.69, P = 0.12; Stress Group M = 0.44; SD = 0.29; Control Group M = 0.23; SD = 0.10]. After a stress-free period of 6 weeks after the end of the stress exposure, we observed that all the psychological changes were restored [Figure 1A; SPP: t₍₅₎ = 3.72, P < 0.05; Stress Group M = 35.50; SD = 2.59; Stress-Recovered Group M = 30.00; SD = 3.03; anxiety score: t₍₅₎ = 2.86, P < 0.05, Stress Group M = 11.00; SD = 7.95; Stress-Recovered Group M = 2.17; SD = 2.71]; depression score: HAD: [t₍₄₎ = 4.84, P < 0.01; Stress Group M = 7.50; SD = 2.59; Stress-Recovered Group M = 2.5; SD = 0.35], except salivary cortisol levels [t₍₅₎ = 0.67, P = 0.53; Stress Group M = 0.44; SD = 0.29; Stress-Recovered Group M = 1.38; SD = 0.8]. Importantly, stress-recovered group did not differ with the control group in all psychological and salivary cortisol measures [PSS: t₍₁₀₎ = −0.08, P = 0.94; HAS: t₍₁₀₎ = −0.58, P = 0.58; HAD: t₍₁₀₎ = −0.85, P = 0.41; (Figure 1B) Cortisol: t₍₁₀₎ = 1.42, P = 0.19].

The ICA analysis revealed the typical spatial pattern of functional connectivity and deactivation in DMN, DAN, VAN, SMN, VN, and AN in all experimental conditions (results not shown).

Whole brain analysis for relative ICVs did not differ between experimental groups. We showed in a previous study (Soares et al., 2013) that exposure to stress triggered a significant reduction in total DMN volume (corrected for ICV) with specific contraction in the left pCC, and bilateral parietal inferior brain regions. Herein, however, we did not find any significant differences in the volume of any of the RSNs between stress participants before and after recovery from stress. No significant areas of expansion or constriction were found in the dorsal and ventral attention networks, SMN, AN, and primary VN between stress and stress-recovered participants (p = 0.99, p = 0.98, p = 0.87, p = 0.84, and p = 0.99, respectively) and between stress-recovered and control groups (p = 0.89, p = 0.54, p = 0.18, p = 0.47, and p = 0.87, respectively).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.