Coupled-channel effects are taken into account for the study of scalar form factors in semi-leptonic $D\to\pi\bar{\ell}\nu_\ell$ and $D\to \bar{K}\bar{\ell}\nu_\ell$ decays, by solving the Muskhelishvili-Omn\`es integral equations. As inputs, we employ the unitarized amplitudes taken from chiral effective theory for the region not far away from thresholds, while, at higher energies of the Goldstone bosons, proper asymptotic conditions are employed. Within Muskhelishvili-Omn\`es formalism, the scalar form factors are represented in terms of Omn\`es matrix multiplied by a vector of polynomials. We reduce the number of subtraction constants by matching to the scalar form factors derived in chiral perturbation theory up to next-to-leading order. The recent lattice QCD data by ETM collaboration for ${D\to\pi}$ and ${D\to\bar{K}}$ scalar form factors are simultaneously well described. The scalar form factors for $D\to\eta$, $D_s\to \bar{K}$ and $D_s\to \eta$ transitions are predicted in their whole kinematical regions. Using our fitting parameters, we also extract the following Cabibbo--Kobayashi--Maskawa elements: $|V_{cd}|=0.243(_{-12}^{+11})_{\rm sta.}({3})_{\rm sys.}(3)_{\rm exp.}$ and $|V_{cs}|=0.950(_{-40}^{+39})_{\rm sta.}(1)_{\rm sys.}(7)_{\rm exp.}$. The approach used in this work can be straightforwardly extended to the semileptonic $B$ decays.