We present constraints on Horndeski gravity from a combined analysis of cosmic shear, galaxy-galaxy lensing and galaxy clustering from $450,mathrm{deg}^2$ of the Kilo-Degree Survey (KiDS) and the Galaxy And Mass Assembly (GAMA) survey. The Horndeski class of dark energy/modified gravity models includes the majority of universally coupled extensions to $Lambda$CDM with one scalar field in addition to the metric. We study the functions of time that fully describe the evolution of linear perturbations in Horndeski gravity. Our results are compatible throughout with a $Lambda$CDM model. By imposing gravitational wave constraints, we fix the tensor speed excess to zero and consider a subset of models including e.g. quintessence and $f(R)$ theories. Assuming proportionality of the Horndeski functions $alpha_B$ and $alpha_M$ (kinetic braiding and the Planck mass run rate, respectively) to the dark energy density fraction $Omega_{mathrm{DE}}(a) = 1 - Omega_{mathrm{m}}(a)$, we find for the proportionality coefficients $hat{alpha}_B = 0.20_{-0.33}^{+0.20} ,$ and $, hat{alpha}_M = 0.25_{-0.29}^{+0.19}$. Our value of $S_8 equiv sigma_8 sqrt{Omega_{mathrm{m}}/0.3}$ is in better agreement with the $Planck$ estimate when measured in the enlarged Horndeski parameter space than in a pure $Lambda$CDM scenario. In our joint three-probe analysis we report a downward shift of the $S_8$ best fit value from the $Planck$ measurement of $Delta S_8 = 0.016_{-0.046}^{+0.048}$ in Horndeski gravity, compared to $Delta S_8 = 0.059_{-0.039}^{+0.040}$ in $Lambda$CDM. Our constraints are robust to the modelling uncertainty of the non-linear matter power spectrum in Horndeski gravity. Our likelihood code for multi-probe analysis in both $Lambda$CDM and Horndeski gravity is publicly available at http://github.com/alessiospuriomancini/KiDSHorndeski .