In this work, the shape of a bluff body is optimized to mitigate velocity fluctuations of turbulent wake flows based on large-eddy simulations (LES). The Reynolds-averaged Navier–Stokes method fails to capture velocity fluctuations, while direct numerical simulations are computationally prohibitive. This necessitates using the LES method for shape optimization given its scale-resolving capability and relatively affordable computational cost. However, using LES for optimization faces challenges in sensitivity estimation as the chaotic nature of turbulent flows can lead to the blowup of the conventional adjoint-based gradient. Here, we propose using the regularized ensemble Kalman method for the LES-based optimization. The method is a statistical optimization approach that uses the sample covariance between geometric parameters and LES predictions to estimate the model gradient, circumventing the blowup issue of the adjoint method for chaotic systems. Moreover, the method allows for the imposition of smoothness constraints with one additional regularization step. The ensemble-based gradient is first evaluated for the Lorenz system, demonstrating its accuracy in the gradient calculation of the chaotic problem. Further, with the proposed method, the cylinder is optimized to be an asymmetric oval, which significantly reduces turbulent kinetic energy and meander amplitudes in the wake flows. The spectral analysis methods are used to characterize the flow field around the optimized shape, identifying large-scale flow structures responsible for the reduction in velocity fluctuations. Furthermore, it is found that the velocity difference in the shear layer is decreased with the shape change, which alleviates the Kelvin–Helmholtz instability and the wake meandering.