We compute multi gluon production at mid rapidity in proton-nucleus collisions within the Color Glass Condensate framework. We show that, in the dilute-dense limit valid for such collisions, the terms responsible for the multi gluon correlation have two origins. On the one hand, the Hanbury-Brown-Twiss interference in the final state. On the other hand, the Bose enhancement of gluons in the projectile and target wave functions, with the latter being suppressed by the number of colors with respect to the former. We also demonstrate that such correlations come from the highest order relevant correlator of Wilson lines in the target wave function, i.e., the quadrupole and sextupole for two and three gluon correlations respectively. We develop a general method for the computation of such high order correlators that captures the bulk of their contribution to the multi gluon production cross section but does not employ the approximation of a large number of colours.