We studied hybrid interactions of Lilium meleagrinum, Lilium gongshanense, and Lilium saluenense using an integrative approach combining population genetics, fieldwork, and phenological research. These three species occur along an elevational gradient, with L. meleagrinum occurring at lower elevations, L. saluenense at higher elevations, and L. gongshanense between them. The species show strong morphological differentiation despite there being no clear environmental barriers to gene flow among them. Lilium gongshanense is likely to have a hybrid origin based on our prior work, but its progenitors remain uncertain. We sought to determine whether gene flow occurs among these three parapatric species, and, if so, whether L. gongshanense is a hybrid of L. meleagrinum and/or L. saluenense. We analyzed data from multiple chloroplast genes and spacers, nuclear internal transcribed spacer (ITS), and 18 nuclear Expressed Sequence Tag-Simple Sequence Repeat (EST-SSR) microsatellites for accessions of the three species representing dense population-level sampling. We also inferred phenology by examining species in the field and using herbarium specimens. We found that there are only two types of chloroplast genomes shared among the three species and that L. gongshanense forms two distinct groups with closest links to other species of Lilium based on ITS. Taken together, L. gongshanense is unlikely to be a hybrid species resulting from a cross between L. meleagrinum and L. saluenense, but gene flow is occurring among the three species. The gene flow is likely to be rare according to evidence from all molecular datasets, and this is corroborated by detection of only one putative hybrid individual in the field and asynchronous phenology. We suspect that the rarity of hybridization events among the species facilitates their continued genetic separation.