Hybrid zones of ecologically divergent populations are ideal systems to study the interaction between natural selection and gene flow during the initial stages of speciation. Here, we perform an amplified fragment length polymorphism (AFLP) genome scan in parallel hybrid zones between divergent ecotypes of the marine snail Littorina saxatilis, which is considered a model case for the study of ecological speciation. Ridged-Banded (RB) and Smooth-Unbanded (SU) ecotypes are adapted to different shore levels and microhabitats, although they present a sympatric distribution at the mid-shore where they meet and mate (partially assortatively). We used shell morphology, outlier and nonoutlier AFLP loci from RB, SU and hybrid specimens captured in sympatry to determine the level of phenotypic and genetic introgression. We found different levels of introgression at parallel hybrid zones and nonoutlier loci showed more gene flow with greater phenotypic introgression. These results were independent from the phylogeography of the studied populations, but not from the local ecological conditions. Genetic variation at outlier loci was highly correlated with phenotypic variation. In addition, we used the relationship between genetic and phenotypic variation to estimate the heritability of morphological traits and to identify potential Quantitative Trait Loci to be confirmed in future crosses. These results suggest that ecology (exogenous selection) plays an important role in this hybrid zone. Thus, ecologically based divergent natural selection is responsible, simultaneously, for both ecotype divergence and hybridization. On the other hand, genetic introgression occurs only at neutral loci (nonoutliers). In the future, genome-wide studies and controlled crosses would give more valuable information about this process of speciation in the face of gene flow.