  In Horndeski theories containing a scalar coupling with the Gauss-Bonnet (GB)
curvature invariant $R_{\rm GB}^2$, we study the existence and linear stability
of neutron star (NS) solutions on a static and spherically symmetric
background. For a scalar-GB coupling of the form $\alpha \xi(\phi) R_{\rm
GB}^2$, where $\xi$ is a function of the scalar field $\phi$, the existence of
linearly stable stars with a nontrivial scalar profile without instabilities
puts an upper bound on the strength of the dimensionless coupling constant
$|\alpha|$. To realize maximum masses of NSs for a linear (or dilatonic) GB
coupling $\alpha_{\rm GB}\phi R_{\rm GB}^2$ with typical nuclear equations of
state, we obtain the theoretical upper limit $\sqrt{|\alpha_{\rm GB}|}<0.7~{\rm
km}$. This is tighter than those obtained by the observations of gravitational
waves emitted from binaries containing NSs. We also incorporate cubic-order
scalar derivative interactions, quartic derivative couplings with nonminimal
couplings to a Ricci scalar besides the scalar-GB coupling and show that NS
solutions with a nontrivial scalar profile satisfying all the linear stability
conditions are present for certain ranges of the coupling constants. In
regularized 4-dimensional Einstein-GB gravity obtained from a Kaluza-Klein
reduction with an appropriate rescaling of the GB coupling constant, we find
that NSs in this theory suffer from a strong coupling problem as well as
Laplacian instability of even-parity perturbations. We also study NS solutions
with a nontrivial scalar profile in power-law $F(R_{\rm GB}^2)$ models, and
show that they are pathological in the interior of stars and plagued by ghost
instability together with the asymptotic strong coupling problem in the
exterior of stars.
