This paper investigates numerically and through an asymptotic approach the three-dimensional stability of steady vertical vortex arrays in a stratified and rotating fluid. Three classical vortex arrays are studied: the Kármán vortex street, the symmetric double row and the single row of co-rotating vortices. The asymptotic analysis assumes well-separated vortices and long-wavelength bending perturbations following Billant (J. Fluid Mech., vol. 660, 2010, p. 354) and Robinson & Saffman (J. Fluid Mech., vol. 125, 1982, p. 411). Very good agreement with the numerical stability analysis is found even for finite wavelength and relatively close vortices. For a horizontal Froude number Fh ≤ 1 and for a non-rotating fluid, it is found that the Kármán vortex street for a street spacing ratio (the distance h between the rows divided by the distance b between vortices in the same row) κ ≤ 0.41 and the symmetric double row for any spacing ratio are most unstable to a three-dimensional instability of zigzag type that vertically bends the vortices. The most amplified vertical wavenumber scales like 1/(bFh) and the growth rate scales with the strain Γ/(2πb2), where Γ is the vortex circulation. For the Kármán vortex street, the zigzag instability is symmetric with respect to the midplane between the two rows while it is antisymmetric for the symmetric double row. For the Kármán vortex street with well-separated vortex rows κ > 0.41 and the single row, the dominant instability is two-dimensional and corresponds to a pairing of adjacent vortices of the same row. The main differences between stratified and homogeneous fluids are the opposite symmetry of the dominant three-dimensional instabilities and the scaling of their most amplified wavenumber. When Fh > 1, three-dimensional instabilities are damped by a viscous critical layer. In the presence of background rotation in addition to the stratification, symmetric and antisymmetric modes no longer decouple and cyclonic vortices are less bent than anticyclonic vortices. However, the dominant instability remains qualitatively the same for the three vortex arrays, i.e. quasi-symmetric or quasi-antisymmetric and three-dimensional or two-dimensional. The growth rate continues to scale with the strain but the most unstable wavenumber of three-dimensional instabilities decreases with rotation and scales like Ro/(bFh) for small Rossby number Ro, in agreement with quasi-geostrophic scaling laws.