In this work we investigate numerically the dynamics of sheared suspensions in the limit of vanishingly small fluid and particle inertia. The numerical model we used is able to handle the multi-body hydrodynamic interactions between thousands of particles embedded in a linear shear flow. The presence of the particles is modelled by momentum source terms spread out on a spherical envelop forcing the Stokes equations of the creeping flow. Therefore all the velocity perturbations induced by the moving particles are simultaneously accounted for. The statistical properties of the sheared suspensions are related to the velocity fluctuation of the particles. We formed averages for the resulting velocity fluctuation and rotation rate tensors. We found that the latter are highly anisotropic and that all the velocity fluctuation terms grow linearly with particle volume fraction. Only one off-diagonal term is found to be non zero (clearly related to trajectory symmetry breaking induced by the non-hydrodynamic repulsion force). We also found a strong correlation of positive/ negative velocities in the shear plane, on a time scale controlled by the shear rate (direct interaction of two particles). The time scale required to restore uncorrelated velocity fluctuations decreases continuously as the concentration increases. We calculated the shear induced self-diffusion coefficients using two different methods and the resulting diffusion tensor appears to be anisotropic too. The microstructure of the suspension is found to be drastically modified by particle interactions. First, the probability density function of velocity fluctuations showed a transition from exponential to Gaussian behaviour as particle concentration varies. Second, the probability of finding close pairs while the particles move under shear flow is strongly enhanced by hydrodynamic interactions when the concentration increases.