Non-linear gyrokinetic simulations of the interchange instability are discussed. The semi-Lagrangian numerical scheme allows one to address two critical points achieved with simulations lasting several confinement times: an accurate statistical analysis of the fluctuations and the back reaction of the turbulence on equilibrium profiles. Zonal flows are found to quench a 2D + 1D interchange turbulence when one of the species has a vanishing response to zonal modes. Conversely, when streamers dominate, the equilibrium profiles are found to be stiff. In the non-linear regime and steady-state turbulence, the distribution function exhibits a significant departure from a Maxwellian distribution. This property is characterized by an expansion on generalized Laguerre functions with a slow decay of the series of moments. This justifies the use of gyrokinetic simulations since a standard fluid approach, based on a limited number of moments, would certainly require a complex closure so as to take into account the impact of these non-vanishing high order moments.