Rhodiola is widely consumed in traditional folk medicine and nutraceuticals. To establish a procedure for the hydrogen (¹H)-NMR spectroscopic fingerprinting of secondary metabolites from three different Rhodiola species, the variation among three Rhodiola species were studied using ¹H-NMR metabolomics combined with multivariate data analysis. Gene expression programming (GEP) was used to generate a formula to distinguish Rhodiola crenulata from two other Rhodiola species. Finally, HPLC was used to demonstrate the results. Same metabolites were compared by quantitative ¹H-NMR (qNMR). Three Rhodiola species were clearly discriminated by ¹H-NMR fingerprinting involved 22 nuclear magnetic signals of chemical constituents. y = d₁₆₆ × 2 + C₁ + d₅₆ + d₂₃₆ − d₁₂₈ × C₂ can be used to distinguish R. crenulata from two other Rhodiola species by GEP. The gallic acid concentration in R. crenulata was significantly higher than in the other. Rhodiola species as was the level of salidroside. R. crenulata also exhibited substantially higher levels of α-glucose. The fatty acid level in Rhodiola kirilowii was lower than the other species. These findings demonstrated that ¹H-NMR fingerprinting combined with principal component analysis (PCA), partial least squares discriminant analysis (PLS-DA), hierarchical cluster analysis (HCA) and GEP can be used to distinguish different Rhodiola species and these methods were applicable and effective approaches for metabolic analysis, species differentiation, and quality assessment. In addition, gallic acid, salidroside, α-D-glucose, glycine, alanine, caffeic acid and tyrosol and are the discriminators.