By employing the QCD factorization approach, we calculated the next-to-leading order new physics contributions to the branching ratios, CP asymmetries, isospin and U-spin symmetry breaking of the exclusive decays $B \to V \gamma$ ( $V = K^*, \rho$ ), induced by the charged Higgs penguins in general two-Higgs-doublet models. Within the considered parameter space, we found that (a) the new physics corrections to the observables are generally small in model I and model III-A, moderate in model II, but large in model III-B; (b) from the well measured branching ratios and upper limits, a lower bound of M H > 200 GeV in model II was obtained, while the allowed range of M H in model III-B is $ 226 \leq M_{H} \leq 293 $ GeV; these bounds are comparable with those from the inclusive $B \to X_s \gamma$ decay; (c) the NLO Wilson coefficient C 7(m b ) in model III-B is positive and disfavored by the measured value of isospin symmetry breaking $ \Delta_{0-}^{\mathrm{{exp}}} (K^*\gamma) = (3.9 \pm 4.8)\%$ , but it still cannot be excluded if we take the large errors into account; (d) the CP asymmetry $\mathcal{A}_{CP}(B \to \rho \gamma)$ in model III-B has an opposite sign to the one in the standard model (SM), which may be used as a good observable to distinguish the SM from model III-B; (e) the isospin symmetry breaking $ \Delta(\rho\gamma)$ is less than $10\%$ in the region of $\gamma = [ 40 \sim 70]^\circ$ preferred by the global fit result, but it can be as large as 20 to $40\%$ in the regions of $\gamma \leq 10^\circ$ and $\gamma \geq 120^\circ$ . The SM and model III-B predictions for $ \Delta(\rho\gamma)$ are opposite in sign for small or large values of the CKM angles; (f) the U-spin symmetry breaking $\Delta U(K^*,\rho)$ in the SM and the general two-Higgs-doublet models is generally small in size: $\sim 10^{-7}$ .