We computed the leading order Wilson coefficients relevant to all the exclusive b→sℓ+ℓ− decays in the framework of the two Higgs doublet model (2HDM) with a softly broken Z2 symmetry by including the O(mb) corrections. We elucidate the issue of appropriate matching between the full and the effective theory when dealing with the (pseudo-)scalar operators for which keeping the external momenta different from zero is necessary. We then make a phenomenological analysis by using the measured B(Bs→μ+μ−) and B(B→Kμ+μ−)high−q2 , for which the hadronic uncertainties are well controlled, and we discuss their impact on various types of 2HDM. A brief discussion of the decays with τ -leptons in the final state is provided too.