The creation of non invasive biomarkers from blood DNA methylation profiles is a cutting-edge achievement in personalized medicine: DNAm epimutations have been demonstrated to be tightly related to lifestyle and environmental risk factors, ultimately providing an unbiased proxy of an individual state of health. At present, the creation of DNAm surrogates relies on univariate penalized regression model, with elastic net being the standard way to-go when accomplishing the task. Nonetheless, more advanced modeling procedures are required when the response is multivariate in nature and the samples showcase a structured dependence pattern. In this work, with the aim of developing a multivariate DNAm biomarker from a multi-centric study, we propose a general framework for high-dimensional, mixed-effects multitask learning. A penalized estimation scheme based on an EM algorithm is devised, in which any penalty criteria for fixed-effects models can be conveniently incorporated in the fitting process. The methodology is then employed to create a novel surrogate of cardiovascular and high blood pressure comorbidities, showcasing better results, both in terms of predictive power and epidemiological interpretation, than state-of-the-art alternatives.