In this paper, the development of three-dimensional fluid interfacial Rayleigh-Taylor (RT) instability coupled with the surface tension was numerically studied using the mesoscopic lattice Boltzmann method. We mainly analyzed the influence of surface tension on fluid interfacial dynamics and spike/bubble late-time growth. The numerical experiments show that there exists the critical surface tension (
$\sigma_{\rm{c}}$
) in the three-dimensional RT instability, above which the RT phenomenon does not appear and below which it would take place. It is found that the critical surface tension increases with the fluid Atwood number and the corresponding numerical predictions show good agreements with those of the theoretical analysis
${\sigma_{\rm{c}}}= {{({{\rho_{\rm{h}}}-{\rho_{\rm{l}}}})g}}/{{{k^2}}}$
. In addition, we can find that increasing surface tension reduces the roll-up of the interface and the complexity of interfacial structure, also preventing the breakup of the interface into the individual droplets. The late-time dynamics of phase interface change from the asymmetric development to the symmetry with respect to the middle axis. When the surface tension is sufficiently low, the spike and bubble amplitudes almost no longer change with it, and further increasing the surface tension can slow down the growth of the spike and bubble amplitudes. Furthermore, we can observe that the development of the high-Reynolds-number RT instability under different surface tensions can also be divided into four distinct stages, including the linear growth, saturated velocity growth, reacceleration, and chaotic mixing. The spike and bubble grow with approximately constant velocities at the saturated stage and their asymptotic values are consistent with those of the modified potential flow theory. In the following, the spike and bubble driven by the increasing Kelvin-Helmholtz vortices are accelerated such that their evolutional velocities exceed the solutions of the potential flow model at the reacceleration stage. The reacceleration stage cannot last infinitely and the spike and bubble velocities at the late time fluctuate with time, implying that the growth of the RT instability enters into the chaotic mixing stage. By numerical analysis, we demonstrate that the three-dimensional RT instability at the chaotic mixing stage has a quadratic growth and also report that the spike and bubble growth rates decrease with the surface tension in general.