Amplified fragment length polymorphism analysis was employed to investigate the population structure of 32 Iranian Silybum marianum populations along with two commercial varieties. A total of 415 polymorphic marker loci were produced by 27 primer combinations with an average of 15.37 markers per combination. Polymorphic information content ranged from 0.24 to 0.44 with an average of 0.35 per primer combination, and marker index was in the range of 2.56–9.50 with an average of 5.37. The average Nei's genetic diversity (HE) and Shannon's diversity index (I) were 0.201 and 0.296, respectively. The coefficient of differentiation among populations (GST) was 0.44, indicating that 44% of the total molecular diversity resulted from differences between populations. We identified three major groups based on cluster analysis and principal coordinate analysis, which were mostly in concordance with the geographical grouping of the populations. The molecular diversity estimate could be useful for selecting appropriate populations to improve S. marianum through conventional and molecular breeding.