The progression of Parkinson's disease (PD) is heterogeneous across patients, affecting counseling and inflating the number of patients needed to test potential neuroprotective treatments. Moreover, disease subtypes might require different therapies. This work uses a data-driven approach to investigate how observed heterogeneity in PD can be explained by the existence of distinct PD progression subtypes. To derive stable PD progression subtypes in an unbiased manner, we analyzed multimodal longitudinal data from three large PD cohorts and performed extensive cross-cohort validation. A latent time joint mixed-effects model (LTJMM) was used to align patients on a common disease timescale. Progression subtypes were identified by variational deep embedding with recurrence (VaDER). In each cohort, we identified a fast-progressing and a slow-progressing subtype, reflected by different patterns of motor and non-motor symptoms progression, survival rates, treatment response, features extracted from DaTSCAN imaging and digital gait assessments, education, and Alzheimer's disease pathology. Progression subtypes could be predicted with ROC-AUC up to 0.79 for individual patients when a one-year observation period was used for model training. Simulations demonstrated that enriching clinical trials with fast-progressing patients based on these predictions can reduce the required cohort size by 43%. Our results show that heterogeneity in PD can be explained by two distinct subtypes of PD progression that are stable across cohorts. These subtypes align with the brain-first vs. body-first concept, which potentially provides a biological explanation for subtype differences. Our predictive models will enable clinical trials with significantly lower sample sizes by enriching fast-progressing patients.