Genome-wide association studies (GWAS) have identified over 170 common breast cancer susceptibility variants using standard GWAS methods. Many of these variants have differential associations by estrogen receptor (ER), but how these variants relate with other tumor features and intrinsic molecular subtypes is unclear. We used two-stage polytomous logistic regression models to evaluate variants in relation to multiple tumor features (ER, progesterone receptor (PR), human epidermal growth factor receptor 2 (HER2) and grade) adjusting for each other, and to intrinsic-like subtypes. Eighty-five of 178 variants were associated with at least one tumor feature (false discovery rate <5%), most commonly ER and grade followed by PR and HER2. Case-control comparisons among these 85 variants identified 65 variants strongly or exclusively associated (P<0.05) with luminal-like subtypes, 5 variants associated with all subtypes at differing strengths and 15 variants primarily associated with non-luminal subtypes, especially triple-negative (TN) disease. Five variants were associated with risk of Luminal A-like and TN subtypes in opposite directions. This report demonstrates a high level of complexity in the etiology heterogeneity of breast cancer susceptibility loci and can inform investigations of subtype-specific risk prediction.