High-dimensional biomarkers, such as gene expression profiles, are often collected longitudinally to monitor disease progression in clinical studies, where the primary endpoint of interest is often a survival outcome. It is of great interest to study the associations between high-dimensional longitudinal biomarkers and the survival outcome as well as to identify biomarkers related to the survival outcome. Joint models, which have been extensively studied in the past decades, are commonly used to study the associations between longitudinal biomarkers and the survival outcome. However, existing joint models only consider one or a few longitudinal biomarkers and cannot deal with high-dimensional longitudinal biomarkers. In this paper we propose a novel penalized joint model that can handle high-dimensional longitudinal biomarkers. Specifically, we impose an adaptive lasso penalty on the parameters for the effects of the longitudinal biomarkers on the survival outcome, which allows for variable selection. We also develop a computationally efficient algorithm for model estimation based on the Gaussian variational approximation method, which can be implemented using the HDJM package in R. Furthermore, based on the penalized joint model, we propose a two-stage selection procedure that can reduce the estimation bias, due to the penalization, and allows for inference. We conduct extensive simulation studies to evaluate the performance of our proposed method. The performance of our proposed method is further demonstrated on a longitudinal gene expression dataset of patients with idiopathic pulmonary fibrosis.