We consider the problem of distributedly estimating Gaussian processes in multi-agent frameworks. Each agent collects few measurements and aims to collaboratively reconstruct a common estimate based on all data. Agents are assumed with limited computational and communication capabilities and to gather M noisy measurements in total on input locations independently drawn from a known common probability density. The optimal solution would require agents to exchange all the M input locations and measurements and then invert an M×M matrix, a non-scalable task. Differently, we propose two suboptimal approaches using the first E orthonormal eigenfunctions obtained from the Karhunen-Loève (KL) expansion of the chosen kernel, where typically E≪M. The benefits are that the computation and communication complexities scale with E and not with M, and computing the required statistics can be performed via standard average consensus algorithms. We obtain probabilistic non-asymptotic bounds that determine a priori the desired level of estimation accuracy, and new distributed strategies relying on Stein's unbiased risk estimate (SURE) paradigms for tuning the regularization parameters and applicable to generic basis functions (thus not necessarily kernel eigenfunctions) and that can again be implemented via average consensus. The proposed estimators and bounds are finally tested on both synthetic and real field data.