We use SHARK, a semi-analytic galaxy formation model, to investigate the physical processes involved in dictating the shape, scatter and evolution of the HI-halo mass relation at $0leq z leq 2$. We compare SHARK with HI clustering and spectral stacking of the HI-halo mass relation derived from observations finding excellent agreement with the former and a deficiency of HI in SHARK at $M_{rm vir}approx 10^{12-13} M_{odot}$ in the latter, but otherwise great agreement below and above that mass threshold. In SHARK, we find that the HI mass increases with the halo mass up to a critical mass of $approx 10^{11.8} M_{odot}$; between $sim 10^{11.8}-10^{13}M_{odot}$, the scatter in the relation increases by 0.7 dex and the HI mass decreases with the halo mass on average; at $M_{rm vir} geq 10^{13} M_{odot}$, the HI content continues to increase with halo mass. We find that the critical halo mass of $approx 10^{12} M_{odot}$ is largely set by feedback from Active Galactic Nuclei (AGN), and the exact shape and scatter of the HI-halo mass relation around that mass is extremely sensitive to how AGN feedback is modelled, with other physical processes playing a less significant role. We determine the main secondary parameters responsible for the scatter of the HI-halo mass relation, namely the halo spin parameter at $M_{rm vir}leq 10^{11.8} M_{odot}$, and the fractional contribution from substructure to the total halo mass for $M_{rm vir}geq 10^{13} M_{odot}$. The scatter at $10^{11.8}<M_{rm vir}<10^{13} M_{odot}$ is best described by the black-hole-to-stellar mass ratio of the central galaxy, reflecting the AGN feedback relevance. We present a numerical model to populate dark matter-only simulations with HI at $0leq z leq 2$ based solely on halo parameters that are measurable in such simulations.