This paper is concerned with the estimation of unknown drift functions of stochastic differential equations (SDEs) from observations of their sample paths. We propose to formulate this as a non-parametric Gaussian process regression problem and use an Itô-Taylor expansion for approximating the SDE. To address the computational complexity problem of Gaussian process regression, we cast the model in an equivalent state-space representation, such that (non-linear) Kalman filters and smoothers can be used. The benefit of these methods is that computational complexity scales linearly with respect to the number of measurements and hence the method remains tractable also with large amounts of data. The overall complexity of the proposed method is O log(N), where N is the number of measurements, due to the requirement of sorting the input data. We evaluate the performance of the proposed method using simulated data as well as with real-data applications to sunspot activity and electromyography.