This study investigates the impact of relative submergence, defined as the ratio of water depth to the diameter of boulders (k = H/D), on turbulence structures in flow through boulder arrays. The large-eddy simulation method is employed to simulate flow through boulder arrays across a range of k values from 0.25 to 3.50. Within this range, three distinct flow regimes are identified: low (k = 0.25), intermediate (k = 0.75 and 1.25), and high (k = 2.0 and 3.5) relative submergence regimes. Across these three regimes, distributions of time-averaged velocities, secondary flow, turbulent kinetic energy, and dominant turbulence structures in the wakes of boulders exhibit significant variations. The wake of boulders, characterized by recirculation flow, only manifests at k >= ${\ge} $ 0.75 and is more pronounced at higher k values. Particularly at k = 3.5, funnel vortices in the wake and secondary flow at the sides of boulders develop, enhancing vertical momentum exchange. Three types of coherent structures are identified within the wake: (a) the near-bed hairpin vortex with a wavelength (lambda $\lambda $) of 0.8D at the lowest k, (b) the lateral flapping of boulder wakes with lambda=2.1D $\lambda =2.1D$ intermediate k, and (c) the meandering of high-speed streaks at the side of boulders with lambda=10D $\lambda =10D$ at high k. These structures alter the distribution of the near-bed Reynolds shear stresses (RSS) and contribute up to 20% of the near-bed RSS. At k <= ${\le} $ 1.25, a region of low near-bed shear stress appears upstream of boulders, while it shifts to the wake of boulders at k = 3.5, contributing the observed variations in deposition patterns at different k values as reported by Papanicolaou et al. (2018, ). In addition, the two bedload periodicities reported in the experiment of Papanicolaou et al. (2018, ) are justified by the ratio of the wavelength of lateral flapping of boulder wakes to that of meandering of low- and high-speed streaks.