In this paper, we introduce the incompressible discrete Maxwellian equilibrium distribution function and external forces into the two-stage third-order Discrete Unified Gas-Kinetic Scheme (DUGKS) for simulating low-speed incompressible turbulent flows with forcing term. The Wall-Adapting Local Eddy-viscosity (WALE) and Vreman sub-grid models for Large-Eddy Simulations (LES) of wall-bounded turbulent flows are coupled within the present framework. In order to simulate the three-dimensional turbulent flows associated with great computational cost, a parallel implementation strategy for the present framework is developed, and is validated by three canonical wall-bounded turbulent flows, viz., the fully developed turbulent channel flow at a friction Reynolds number (Re) about 180, the turbulent plane Couette flow at a friction Re number about 93 and three-dimensional lid-driven cubical cavity flow at a Re number of 12000. The turbulence statistics are computed by the present approach with both WALE and Vreman models, and their predictions match precisely with each other. Especially, the predicted flow physics of three-dimensional lid-driven cavity are consistent with the description from abundant literatures. While, they have small discrepancies in comparison to the Direct Numerical Simulation (DNS) due to the relatively low grid resolution. The present numerical results verify that the present two-stage third-order DUGKS-based LES method is capable for simulating inhomogeneous wall-bounded turbulent flows and getting reliable results with relatively coarse grids.