During the Cosmic Dawn (CD), the HI 21-cm optical depth (
τ ) in the intergalactic medium can become significantly large. Consequently, the second and higher-order terms of
τ appearing in the Taylor expansion of the HI 21-cm differential brightness temperature (
δTb ) become important. This introduces additional non-Gaussianity into the signal. We study the impact of large
τ on statistical quantities of HI 21-cm signal using a suite of standard numerical simulations that vary X-ray heating efficiency and the minimum halo mass required to host radiation sources. We find that the higher order terms suppress statistical quantities such as skewness, power-spectrum and bispectrum. However, the effect is found to be particularly strong on the non-Gaussian signal. We find that the change in skewness can reach several hundred percent in low X-ray heating scenarios, whereas for moderate and high X-ray heating models changes are around
∼40% and
60%, respectively, for
Mh,min=109M⊙. This change is around
∼75%,
25% and
20% for low, moderate and high X-ray heating models, respectively, for
Mh,min=1010M⊙. The change in bispectrum in both the halo cutoff mass scenarios ranges from
∼10% to
∼300% for low X-ray heating model. However, for moderate and high X-ray heating models the change remains between
∼10% to
∼200% for both equilateral and squeezed limit triangle configuration. Finally, we find that up to third orders of
τ need to be retained to accurately model
δTb, especially for capturing the non-Gaussian features in the HI 21-cm signal.