Distinct multimodal biological and functional profiles of symptom-based subgroups in recent-onset psychosis

Symptom heterogeneity characterizes psychotic disorders and hinders the delineation of underlying biomarkers. Here, we identify symptom-based subtypes of recent-onset psychosis (ROP) patients from the multi-center PRONIA (Personalized Prognostic Tools for Early Psychosis Management) database and explore their multimodal biological and functional signatures. We clustered N = 328 ROP patients based on their maximum factor scores in an exploratory factor analysis on the Positive and Negative Syndrome Scale items. We assessed inter-subgroup differences and compared to N = 464 healthy control (HC) individuals regarding gray matter volume (GMV), neurocognition, polygenic risk scores, and longitudinal functioning trajectories. Finally, we evaluated factor stability at 9- and 18-month follow-ups. A 4-factor solution optimally explained symptom heterogeneity, showing moderate longitudinal stability. The ROP-MOTCOG (Motor/Cognition) subgroup was characterized by GMV reductions within salience, control and default mode networks, predominantly throughout cingulate regions, relative to HC individuals, had the most impaired neurocognition and the highest genetic liability for schizophrenia. ROP-SOCWD (Social Withdrawal) patients showed GMV reductions within medial fronto-temporal regions of the control, default mode, and salience networks, and had the lowest social functioning across time points. ROP-POS (Positive) evidenced GMV decreases in salience, limbic and frontal regions of the control and default mode networks. The ROP-AFF (Affective) subgroup showed GMV reductions in the salience, limbic, and posterior default-mode and control networks, thalamus and cerebellum. GMV reductions in fronto-temporal regions of the salience and control networks were shared across subgroups. Our results highlight the existence of behavioral subgroups with distinct neurobiological and functional profiles in early psychosis, emphasizing the need for refined symptom-based diagnosis and prognosis frameworks.


Introduction
Psychosis is an umbrella term for diverse and dynamically evolving clinical phenotypes that affect the processing of perceived information about reality.This heterogeneity hinders the development of optimal biologically informed early diagnostic and therapeutic tools in line with the scope of precision psychiatry [1], [2].As an alternative to the classic psychiatric nosology (e.g., Diagnostic And Statistical Manual Of Mental Disorders, [3]), approaches such as the Research Domain Criteria (RDoC, [4]) or the Hierarchical Taxonomy Of Psychopathology (HiTOP, [5]) propose focusing on lower-order, symptom-based disease constructs, which may be more closely related to speci c neurobiological pathways.
In schizophrenia, the distinction between positive (e.g., hallucinations, delusions), negative (e.g., blunted affect, avolition, anhedonia), and cognitive/general symptoms (e.g., disorganized speech and thought) has proven helpful in delineating the symptom complexity observed in clinical practice [6], [7], based on widely used scales such as the Positive and Negative Syndrome Scale (PANSS, [8]).In the early stages of psychosis, four to ve PANSS symptoms-based subgroups have been commonly reported [7], [9]- [14], with symptom dimensions shown to dynamically evolve [15], [16]) and associate with functional ( [10], [14]) and treatment outcomes (Martinuzzi et al., 2019).Importantly, psychotic symptom dimensions have been associated with different layers of biological variability, suggesting their proximity to pathophysiological processes, in line with RDoC.Speci cally, neurocognitive functions closely linked to speci c brain networks [17]- [19] are mainly associated with negative symptoms in early-stage psychosis [20]- [23].Furthermore, symptom-dependent structural brain abnormalities have been reported in patients with established schizophrenia [24]- [27].For instance, Koutsouleris et al. (2008) [24] found that patients with predominantly positive, negative, or disorganized symptoms, as identi ed using a PANSS-factorization approach, showed distinct patterns of gray matter density reductions but also cross-subtype alterations within prefrontal-perisylvian regions.Evidence for relationships between gray matter volume (GMV) alterations and speci c symptom categories also exists in patients with rst-episode psychosis [28]- [34], pointing towards multiple potential neurobiological pathways that pleiotropically map to psychosis.Moreover, the genetic liability for schizophrenia has been differentially associated with symptom dimensions both in patients with chronic schizophrenia [35], [36] and in rst-episode psychosis [37], [38].Speci cally, signi cant associations between polygenic risk scores for schizophrenia and both negative [38], positive [38], and general PANSS-based symptoms [37] were reported in rst-episode psychosis by different studies, highlighting the complex and still insu ciently delineated associations between genetic factors and distinct symptom constellations.
Despite the existing evidence that psychotic symptom dimensions are differentially associated with biological signatures, we currently lack a comprehensive multi-layer characterization of early-stage symptom-based psychosis subgroups.Here, we aimed to address this gap by exploring multimodal biological signatures of symptom-based subtypes of patients with recent-onset psychosis (ROP) from the multi-site European PRONIA study (Personalized Prognostic Tools for Early Psychosis Management; https://www.proniapredictors.eu/pronia/index.html).Concretely, we examined (1) the PANSS-based symptom structure within ROP patients and its longitudinal stability over 9-and 18-month follow-up time points using a data-driven factorization approach, (2) characterized patient subgroups de ned by their most pronounced symptom dimension in terms of neurocognition, GMV, and polygenic risk scores (PRS) to obtain a deep multimodal delineation of disease heterogeneity in early illness stages, and (3) compared longitudinal functioning trajectories of the obtained patient subgroups to assess their prognostic value.

Study Participants
The total sample consisted of 328 ROP patients and 464 HC individuals (see consort chart, Figure S1), recruited across nine sites in Finland, Germany, Italy, Switzerland, and the United Kingdom as part of the longitudinal, multi-site PRONIA study.For assessing the longitudinal stability of the factor analysis solution, we used data collected at three time points: baseline (T0), 9-month (T1), and 18-month follow-up (T2).For further information on the PRONIA project, including inclusion and exclusion criteria, see Table S1.Adult participants gave written informed consent, and patients younger than 18 and their guardians gave written informed assent and consent, respectively.The study was registered at the German Clinical Trials Register (DRKS00005042) and approved by the research ethics committees of all sites.

Exploratory factor analysis and factor model selection
We analyzed patterns of co-occurring symptoms at baseline (T0) by employing a maximum-likelihood-based exploratory factor analysis (EFA) of the PANSS [8] items using R [39], R Studio [40], and the psych (2.1.6)package in R [41].EFA is a dimensionality reduction technique that aims to nd latent dimensions underlying the observed data [42].As factor loadings depend on the speci ed number of factors, we compared solutions with one to ve factors.To account for correlated symptom dimensions, we used the oblique promax rotation method [43].We selected the besttting model based on the lowest average Bayesian Information Criterion (BIC) after jackknife resampling (n = 327) [44].The BIC is a measure of model t that introduces a penalty term relative to the number of parameters in the model [45], addressing the issue that more complex models are more prone to over tting [46].Jackknife resampling was employed to evaluate the stability of the factor solutions and control for sampling biases on model t estimates [47].We tested signi cant differences between models using paired t-tests and the Mann-Whitney U test between the BIC distributions derived during jackknife resampling.To investigate the stability of the factors over time, we computed the Pearson correlation coe cients (r) between item factor loadings across time points and assessed the percentage of items loading highest on the same factor at the different time points using the dice coe cient.
ROP patients were assigned to subgroups via majority vote during jackknife resampling based on their maximum score on the symptom dimensions of the best-tting model.Assignment certainty was quanti ed by how often a subject was assigned to a group relative to another group during jackknife resampling.We assessed the longitudinal stability of the identi ed symptom dimensions (factor scores) through the relative number of patients assigned to the same symptom group at T1 and T2 respectively, by applying the optimal factor model found at T0, and the Pearson correlations of factor dimensions across time points.

Acquisition and preprocessing of structural MRI data
Neuroimaging data was available for 328 ROP patients and 464 HC individuals.The PRONIA study aimed to capture the natural heterogeneity of MRI sequences encountered in clinical settings and, therefore, required minimal harmonization of magnetic resonance imaging (MRI) scanners across sites.The T1-weighted structural sequence involved (1) acquisition in 1mm 3 resolution, (2) enforcement of Field of View parameters that ensured a full 3D coverage of the brain, (3) maximization of the signal-to-noise ratio by choice of the relaxation-and echo-time and other related parameters (for details refer to Table S2).
At each study site, acquired T1-weighted images were visually inspected, anonymized, and defaced using an in-house script based on the Freesurfer toolbox.Subsequently, MRI scans were pre-processed with the open-source CAT12 toolbox (version r1207; http://dbm.neuro.unijena.de/cat12/),an extension of SPM12 using a standardized pipeline.Images were segmented into white matter, gray matter, and cerebrospinal uid and registered to the MNI-152 space using the DARTEL algorithm [48].Gray matter tissue maps were modulated using the Jacobian determinants from the registration step, producing GMV maps in standard space.These GMV maps were smoothed using a 4mm 3 full-width-athalf-maximum Gaussian kernel (for details regarding these steps, see Supplementary Methods).Image quality checks were done by applying CAT12's quality assurance framework [49], which produces a score from excellent (A) to failed (F) for each image.Images of nine participants with a score of C (satisfactory) or lower were excluded from the analysis.

Subgroup characterization: Neurocognition, clinical and sociodemographic data
Neurocognitive performance was assessed by mapping several cognitive batteries available in the PRONIA study onto the six domains from the MATRICS consensus cognitive battery: social cognition, working memory, speed of processing, verbal learning, reasoning, and attention [50], [51], as described in the Supplementary Methods and Table S3.
We tested for neurocognitive, as well as sociodemographic and clinical group-level differences between HC individuals and ROP patients (regardless of subgroup) using, for continuous variables, two-sample t-tests or Mann Whitney U test in case the assumption of normality was not ful lled, and, for binary variables,  2 test.We further tested for group differences between HC individuals and each subgroup, employing the same tests with Benjamini-Hochberg correction for the number of tests conducted for each variable (i.e., the number of subgroups).Lastly, we compared the ROP subgroups against each other in terms of the same variables, and the distribution of different ICD-10 diagnoses within each group, using ANOVA or Kruskal-Wallis and  2 tests.Post-hoc comparisons for ROP-subgroups were conducted using Tukey HSD or Dunn tests.

Subgroup characterization: Neuroimaging data analysis
Spatially smoothed GMV images were corrected for site/batch effects using a MATLAB-based implementation of the ComBat algorithm [52], [53] while preserving the effect of the patient subgroups.Then, corrected images were compared using a non-parametric permutation-based inference strategy as implemented in the Permutation Analysis of Linear Models software (PALM, a119, [54] running in MATLAB R2020b [55]).Statistical contrasts included bidirectional tests between HC individuals and all ROP patients, between HC individuals and each ROP factor subgroup individually, and between each pair of ROP subgroups.Additionally, we conducted a conjunction test computed as the maximum of the P-values of individual contrasts to identify signatures common for all ROP subgroups relative to the HC individuals.Mean-centered total-intracranial volume, age, and sex were entered as covariates into the design matrix.To ensure the tests' validity, we de ned exchangeability blocks for each site and performed within-site permutations for each contrast [56].Family-wise error correction across the whole brain was performed using synchronized permutations [57], [58] of the threshold-free-cluster-enhancement [59] of the unpooled variance Welch t-statistic.We performed 1000 permutations and approximated the tail of the maximum permutation distribution using a generalized Pareto distribution [57].The signi cance level was de ned at P FWE <.05.Lastly, we quanti ed the percentage of statistically signi cant voxels within large-scale brain networks as de ned based on the YEO 7 Networks atlas [60] and the Buckner atlas for cerebellum parcellation [61].

Subgroup characterization: Genotyping and Polygenic Risk Score calculation
Associations of PRS with overall ROP status and each ROP subgroup individually were estimated using binary logistic regression models (with HC individuals as reference) adjusted for sex and the rst ten ancestry-informative principal components.Nagelkerke's pseudo-R 2 was derived and corrected for the covariates by substituting the null model in Nagelkerke's equation for the model, including the covariates.The corrected pseudo-R 2 obtained was then re-scaled to the liability scale [62], obtaining a value directly comparable with heritability and robust against ascertainment bias.A linear transformation on the liability scale was based on lifetime risk (K = .03for overall ROP status [63], and for each subgroup, K = .03was proportionally adjusted to the relative group size with a factor of .25.Within each set of analysis, Benjamini-Hochberg correction was employed to account for the assessed PRS P-thresholds.

Subgroup characterization: Longitudinal functioning trajectories
We employed linear mixed-effects models to account for the correlated nature of repeated measures within participants to investigate the association between the ROP subgroup and functioning across the follow-up time points.The models included xed effects for subgroup membership, time point, age, and sex, with subject ID incorporated as a random effect.Functioning outcomes, including current Global Functioning Scale [64] ratings for Role (GF:Role) and Social (GF:Social), were considered as dependent variables.Primary outcomes focused on the coe cients associated with subgroups, indicating their impact on the speci ed dependent variables.All models were estimated using the lmer function from the lme4 package in R [65] after assessing the assumptions of the mixed-effects models, including checks for collinearity, residuals normality, and homoscedasticity.Additionally, pairwise contrasts between the ROP subgroups were analyzed using the emmeans package in R [66] with Tukey p-value adjustment for multiple comparisons.To further examine differential associations of subgroup membership and functioning over time, we added an interaction effect of ROP subgroup and timepoint to the models.

Gray matter volume signatures of ROP subgroups
Independent of the PANSS-based symptom subgroup, all ROP patients showed signi cant GMV reductions across large-scale brain networks such as the control, default mode, and salience networks (Fig. 1A, 2) compared to HC individuals.
Regarding differences between each ROP subgroup and HC individuals, the ROP-MOTCOG subgroup was characterized by extended GMV reductions mainly within cingulate regions of the salience and control networks, as well as regions of the default mode, dorsal attention and limbic (orbitofrontal cortex) networks relative to HC individuals (Fig. 1A, 2).The ROP-SOCWD subgroup showed overlapping GMV decreases to the ROP-MOTCOG subgroup relative to HC; however, outside the cingulate gyrus and more extended in the medial prefrontal cortex (Fig. 1A, 2).The ROP-POS subgroup was characterized by GMV decreases mainly within the thalamus and anterior regions of the limbic (orbitofrontal cortex) and salience networks compared to HC individuals (Fig. 1A, 2).Lastly, the ROP-AFF subgroup showed GMV decreases in the thalamus, cerebellum, insula, and posterior regions of the control and default-mode network (Fig. 1A, 2).Furthermore, the conjunction analysis revealed GMV reductions within the salience (insula), default-mode (parahippocampus, PFC), limbic and frontotemporal regions within the control networks that were shared among all ROP subgroups relative to the HC individuals (Fig. 1B and 2).Finally, no signi cant GMV differences were found when comparing the ROP subgroups against each other.

Polygenic risk score characterization of ROP subgroups
The binary logistic regression analyses estimating the association between PRS for schizophrenia and a diagnosis of ROP, regardless of subgroup, revealed signi cant (P < .05)PRS effects across all 10 PRS thresholds (see Table S11).Irrespective of the threshold, a higher PRS increased the odds of ROP diagnosis.The PRS explained between 1.05-4.82% of the variance (re-scaled, corrected pseudo-R2) in ROP liability across thresholds (see Fig. 3).The PRS, including SNPs associated with SCZ at P < .2,explained 4.82% of the liability variance of ROP (OR = Inf, P < .001).The same PRS explained a similar amount of liability variance of the ROP-MOTCOG (5.32% variance explained, OR = Inf, P < .001)and ROP-POS subgroup (4.86% variance explained, OR = Inf, P < .001),and relatively smaller proportions of liability variance for ROP-SOCWI (2.86% variance explained, OR = Inf, P < .001),and ROP-AFF subgroup (2.74% variance explained, OR = Inf, P < .001)(see Fig. 3).

Longitudinal functioning trajectories of ROP subgroups
The mixed-effects model with GF:Social as the dependent variable without interaction between group and time point revealed a signi cant effect of the ROP (F = 3.53, P = .015,see Table S12), highlighting overall differences among the four groups.Similarly, the effect of time point was signi cant (F = 73.23,P < .001).When looking at the pairwise contrasts between the ROP subgroup throughout the whole follow-up period, ROP-SOCWD showed signi cantly lower scores in GF:Social than ROP-AFF (estimate=-0.68,P = .008,see Table S13).In the GF:Social model including the subgroup and time interaction, while signi cant ROP subgroup effects (F = 3.00, P = .030),time point effects (F = 70.52,P < .001)were found, the interaction between ROP group and time point was not statistically signi cant (F = .59,P = .740,see Table S12), suggesting that the impact of ROP subgroup on GF:Social did not change across different time points (see Fig. 4).
Neither of the GF:Role models revealed a signi cant effect of the ROP subgroup membership, indicating no overall differences among the four groups in GF:Role over the period (see Table S12).When including the interaction term between group and time, this effect was not signi cant, suggesting that the impact of the ROP subgroup on GF:Role did not change across time points (see Fig. 4).All model coe cients are depicted in Table S13.

Discussion
psychosis stages are characterized by diverse clinical manifestations, complicating the development of biologically informed diagnostic and treatment approaches [67].In this context, a focus on clinical symptoms as proxies for potentially divergent pathophysiological pathways holds promise for advancing psychosis psychiatric care.Here, we provide novel evidence for distinct multimodal biological signatures underlying symptom-based subgroups within a multi-site longitudinal recent-onset psychosis (ROP) cohort.
Our data-driven factorization approach using symptom scores (PANSS) revealed four ROP subgroups with distinct clinical pro les characterized by Motor/Cognition, Social Withdrawal, Positive, and Affective symptoms.Our results are partially aligned with previous PANSS factorization studies, which identi ed between three and seven distinct factors in the early stages of psychosis [7], [9]- [13], [68], [69].The variability across studies could result from different methodological approaches, sample composition, and statistical power.In this context, the optimal granularity of symptom dimensions for clinical translation may require evaluation relative to their mapping onto distinct biological pathways, and their relationship to prognostic and theragnostic outcomes.Furthermore, we show moderate stability of the predominant symptom dimensions for the ROP patients over 18 months but also signi cant positive associations between the dimensions, which become stronger at the latest follow-up time point.These results support the idea of a crystallization of phenotypes along the disease trajectory [58], [59], [61]-[66] and a potential aggregation of the symptomatologic load over time, driven by the complex interactions and reciprocal loops between the brain systems underlying distinct symptoms.Importantly, we found distinct biological signatures characterizing the ROP subgroups.The ROP-MOTCOG subgroup showed the lowest taskbased processing speed and working memory performances, and extensive GMV decreases throughout the salience, default-mode, fronto-parietal networks, especially pronounced within the cingulate cortex.This pattern aligns with previous research showing gray matter volume reductions within the salience and control networks being associated with cognitive impairments in rst-episode psychosis [28], [70]- [72].The ROP-MOTCOG patients also showed the highest genetic liability for schizophrenia and the highest longitudinal stability of their subgroup membership from baseline to the follow-ups, suggesting a more "trait"-like symptom dimension in comparison to the other groups.Collectively, ndings for the ROP-MOTCOG subtype align with the concept of "de cit schizophrenia", characterized by severe and persistent cognitive de cits and structural alterations within frontotemporal brain systems [73]- [75].
The ROP-SOCWD subgroup showed overlapping GMV reductions with the ROP-MOTCOG subgroup outside the cingulate cortex and more extended to the medial prefrontal cortex, consistent with previous associations between medial prefrontal regions and emotion regulation and social cognition [76]- [78].Additionally, this group showed the highest impairment of social functioning across time points, consistent with the socio-emotional nature of the negative symptoms included in this dimension.These results highlight the potential prognostic value of the ROP-SOCWD for poor social outcomes and encourage its further investigation within personalized strati cation and intervention frameworks.
For the ROP-POS patients, we found GMV decreases in the thalamus, as well as limbic (orbitofrontal cortex) and salience networks, in line with previous reports of associations between positive psychotic symptoms and GMV reductions in the thalamus [24], [33], [34], and both functional and structural abnormalities in regions of the limbic [79], [80] and salience networks [81], [82].Lastly, the ROP-AFF subgroup showed GMV decreases in the thalamus, cerebellum and posterior regions of the default-mode, salience and control networks, which have been previously linked to mood disorders [83]-[87].
Additionally, our results revealed GMV decreases in the fronto-temporal regions of the salience and control networks common among all ROP patient subgroups relative to the HC individuals, highlighting a potential "core" psychosis biomarker.This is consistent with previous reports of salience network abnormalities as a key nding across psychosis-spectrum disorders [81], [82], as well as studies highlighting structural frontotemporal alterations independent of symptom-based subtypes in chronic schizophrenia patients [24].
While our study provides valuable insights, it's important to acknowledge certain limitations.First, despite the suitability of subtyping approaches for clinical translation, the categorical assignment of subjects to factors based on their highest loading scores may impose hard boundaries to the continuous nature of psychiatric symptoms as proposed by modern nosological frameworks such as the RDoC [4], [88].Future research may bene t from exploring the dimensional nature of the factors identi ed here, such that the exact clinical manifestations of the individual emanate from the interactions and aggregation of these distinct dimensions at each time point [89], [90].Second, we have not explored subgroupdependent longitudinal brain changes in the current analyses.Investigating whether early psychosis clinical subtypes are linked to divergent trajectories of brain volumetric degeneration/recovery would thus provide further insight into their prognostic value, as well as underlying pathophysiological mechanisms.Along similar lines, further research is required to establish the clinical validity and utility of such bio-behavioral subtypes by elucidating their individual-level prognostic value relative to clinically relevant outcomes such as daily functioning, treatment response or remission rates.
In conclusion, our ndings support the existence of behavioral ROP subgroups characterized by distinct brain structural, genetic, neurocognitive and longitudinal functioning pro les.Our results highlight the potential of symptom-based subtyping for facilitating the development of biologically informed diagnostic, prognostic, and therapeutic strategies in early-stage psychosis.

Table 1
Statistical comparison of sociodemographic, cognitive, and clinical characteristics between the HC individuals and ROP patients, as well as between the four ROP subgroups at baseline.