The 21 cm intensity mapping experiments promise to obtain the large-scale distribution of HI gas at the post-reionization epoch. In order to reveal the underlying matter density fluctuations from the HI mapping, it is important to understand how HI gas traces the matter density distribution. Both nonlinear halo clustering and nonlinear effects modulating HI gas in halos may determine the scale below which the HI bias deviates from linearity. We employ three approaches to generate the mock HI density from a large-scale N-body simulation at low redshifts, and demonstrate that the assumption of HI linearity is valid at the scale corresponding to the first peak of baryon acoustic oscillations, but breaks down at $k gtrsim 0.1,h, {rm Mpc}^{-1}$. The nonlinear effects of halo clustering and HI content modulation counteract each other at small scales, and their competition results in a model-dependent sweet-spot redshift near $z$=1 where the HI bias is scale-independent down to small scales. We also find that the linear HI bias scales approximately linearly with redshift for $zle 3$.