Context. M dwarfs are active stars that exhibit variability in chromospheric emission and photometry at short and long timescales, including long cycles that are related to dynamo processes. This activity also impacts the search for exoplanets because it affects the radial velocities. Aims. We analysed a large sample of 177 M dwarfs observed with HARPS during the period 2003–2020 in order to characterise the long-term variability of these stars. We compared the variability obtained in three chromospheric activity indices (Ca II H & K, the Na D doublet, and Hα) and with ASAS photometry. Methods. We focused on the detailed analysis of the chromospheric emission based on linear, quadratic, and sinusoidal models. We used various tools to estimate the significance of the variability and to quantify the improvement brought by the models. In addition, we analysed complementary photometric time series for the most variable stars to be able to provide a broader view of the long-term variability in M dwarfs. Results. We find that most stars are significantly variable, even the quietest stars. Most stars in our sample (75%) exhibit a long-term variability, which manifests itself mostly through linear or quadratic variability, although the true behaviour may be more complex. We found significant variability with estimated timescales for 24 stars, and estimated the lower limit for a possible cycle period for an additional 9 stars that were not previously published. We found evidence of complex variability because more than one long-term timescale may be present for at least 12 stars, together with significant differences between the behaviour of the three activity indices. This complexity may also be the source of the discrepancies observed between previous publications. Conclusions. We conclude that long-term variability is present for all spectral types and activity level in M dwarfs, without a significant trend with spectral type or mean activity level.