Chemotaxonomy presents various challenges that need to be overcome in order to obtain valid and reliable results. Individual genetic and environmental variations can give a false picture and lead to wrong conclusions. Applying a holistic approach, based on multivariate data analysis, these challenges can be overcome. Thus, a metabolomics approach has to be optimized depending on the subject of research. We used 1H NMR-based metabolomics as a potential chemotaxonomic tool on the selected Euphorbia species growing wild in Serbia. Principal components analysis (PCA), soft independent modeling by class analogy (SIMCA) and Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA) were used to analyze obtained NMR data in order to reveal chemotaxonomic biomarkers. The standard protocol for plant metabolomics was optimized aiming to extract more specific metabolites, which are characteristic for the Euphorbia genus. The obtained models were validated, which revealed that variables unique for each species were associated with certain classes of molecules according to literature data. In E. salicifolia, acacetin-7-O-glycoside (not found before in the species) was detected, and the structure of the aglycone part was solved based on 2D NMR data. In the presented paper, we have shown that metabolomics can be successfully used in Euphorbia chemotaxonomy.