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 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. RB (Ridged-Banded) and SU (Smooth-Unbanded) 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 non-outlier 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 non-outlier 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 QTLs 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 hybridisation. On the other hand, genetic introgression occurs only at neutral loci (non-outliers). 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.