  We present the first large-scale, full radiative transfer simulations of the
reionization of the intergalactic medium in the presence of Lyman-limit systems
(LLSs). To illustrate the impact of LLS opacity, possibly missed by previous
simulations, we add either a uniform or spatially-varying hydrogen bound-free
opacity. This opacity, implemented as the mean free path (mfp) of the ionizing
photons, extrapolates the observed, post-reionization redshift dependence into
the epoch of reionization. In qualitative agreement with previous studies, we
find that at late times the presence of LLSs slows down the ionization fronts,
and alters the size distribution of H II regions. We quantitatively
characterize the size distribution and morphological evolution of H II regions
and examine the effects of the LLSs on the redshifted 21-cm signal from the
patchy reionization. The presence of LLSs extends the ionization history by
$\Delta z \sim 0.8$. The LLS absorbers significantly impede the late-time
growth of the H II regions. The position dependent LLS distribution slows
reionization further and additionally limits the late growth of the ionized
regions. However, there is no "freeze out" of the H II regions and the largest
regions grow to the size of the simulation volume. The 21-cm power spectra show
that at large scales the power drops by a factor of 2 for 50% and 75%
ionization stages (at $k = 0.1$ $\text{h} \, \text{Mpc}^{-1} $) reflecting the
limiting effect of the LLSs on the growth of ionized patches. The statistical
observables such as the RMS of the brightness temperature fluctuations and the
peak amplitudes of the 21-cm power spectra at large-scales ($k = 0.05 - 0.1$
$\text{h} \, \text{Mpc}^{-1} $) are diminished by the presence of LLS.
