In recent years, coupled with traditional turbulence models, the second-order gas-kinetic scheme (GKS) has been used in the turbulent flow simulations. At the same time, high-order GKS has been developed, such as the two-stage fourth-order scheme (S2O4) GKS, and used for laminar flow calculations. In this paper, targeting on the high-Reynolds number engineering turbulent flows, an implicit high-order GKS with Lower-Upper Symmetric Gauss-Seidel (LU-SGS) technique is developed under the S2O4 framework. Based on Vreman-type LES model and k−ω SST model, a turbulent relaxation time is obtained and used for the determination of an enlarged particle collision time for the high-Reynolds number turbulent flow simulation. Numerical experiments include incompressible decaying homogeneous isotropic turbulence, incompressible high-Reynolds number flat plate turbulent flow, incompressible turbulence around NACA0012 airfoil, transonic turbulence around RAE2822 airfoil, and transonic high-Reynolds number ARA M100 wing-body simulation. Comparisons among the numerical solutions from current implicit high-order GKS, the explicit high-order GKS, the implicit second-order GKS, and experimental measurements have been conducted. Through these examples, it is concluded that the high-order GKS has high accuracy in space and time, especially for smooth flows, and obtains more accurate turbulent flow fields on coarse grids than the second-order GKS. In addition, significant acceleration on computational efficiency, as well as super robustness in simulating complex flows are confirmed from the current implicit high-order GKS. This study also indicates that turbulence modeling plays a dominant role in capturing physical solution, such as in the transonic three-dimensional complex RANS simulation, which is beyond the numerical discretization error, such as the differences between the second and fourth-order GKS.