  We study the dependence of galaxy clustering on atomic gas mass using a
sample of $\sim$16,000 galaxies with redshift in the range of $0.0025<z<0.05$
and HI mass of $M_{\rm HI}>10^8M_{\odot}$, drawn from the 70% complete sample
of the Arecibo Legacy Fast ALFA survey. We construct subsamples of galaxies
with $M_{\rm HI}$ above different thresholds, and make volume-limited
clustering measurements in terms of three statistics: the projected two-point
correlation function, the projected cross-correlation function with respect to
a reference sample selected from the Sloan Digital Sky Survey, and the
redshift-space monopole moment. In contrast to previous studies, which found
no/weak HI-mass dependence, we find both the clustering amplitude on scales
above a few Mpc and the bias factors to increase significantly with increasing
HI mass for subsamples with HI mass thresholds above $10^9M_{\odot}$. For HI
mass thresholds below $10^9M_{\odot}$, while the measurements have large
uncertainties caused by the limited survey volume and sample size, the inferred
galaxy bias factors are systematically lower than the minimum halo bias factor
from mass-selected halo samples. The simple halo model, in which galaxy content
is only determined by halo mass, has difficulties in interpreting the
clustering measurements of the HI-selected samples. We extend the simple model
by including the halo formation time as an additional parameter. A model that
puts HI-rich galaxies into halos that formed late can reproduce the clustering
measurements reasonably well. We present the implications of our best-fitting
model on the correlation of HI mass with halo mass and formation time, as well
as the halo occupation distributions and HI mass functions for central and
satellite galaxies. These results are compared with the predictions from
semi-analytic galaxy formation models and hydrodynamic galaxy formation
simulations.
