This paper investigates the three-dimensional instabilities and the optimal perturbations on a pair of horizontal counter-rotating Lamb-Oseen vortices in a vertically stably stratified flow. Two-dimensional (2D) simulations are first performed, showing that while the dipole moves vertically against the stratification the vortex parameters: the radius a∗, the separation distance b∗, and the circulation Γ∗ are solely function of the time rescaled by the Brunt-Väisälä frequency N, independently of the Froude number, when the Reynolds number is large enough. Here, the Froude number is Fr = W0/Nb0 with W0 the initial advection velocity of the dipole and b0 the initial separation distance between the two vortices. For weak and moderate stratifications (large Fr), the stratification acts on a long time scale compared to the advection time of the dipole implying that the 2D flow can be considered as quasi-steady. In that case, when three dimensional instabilities are added, a linear stability analysis of this 2D flow at different instants retrieves the instability peaks corresponding to the Crow instability for the long wavelengths and to the elliptic instability for the short wavelengths showing that the dynamics is almost unaffected by buoyancy effects. The Crow and elliptic instabilities scale with the instantaneous dipole parameters showing in particular that stratification promotes instability by reducing the distance b∗ between vortices [K. K. Nomura et al., “Short-wavelength instability and decay of a vortex pair in stratified fluid,” J. Fluid Mech. 553, 283-322 (2006)]. For strong stratifications (Froude numbers of order unity or smaller), the quasi-steady approximation is not valid, and the question of stability should be formulated in a different way, by, for example, searching for the transient growth of the energy of perturbation that may be computed for steady or unsteady base flow. Then, for each time horizon τ, we should determine the critical perturbation leading to the largest energy growth by the time τ. Presently, we compute the optimal perturbations at two time horizons τ = 4 and τ = 10 dimensionalized by 2πb02/Γ0 with a direct-adjoint technique which takes into account the evolution of the base flow. In the homogeneous case, this technique allows to investigate the effect of the weak unsteadiness of the flow due to viscous diffusion which induces a growth of the vortex core radius a∗. Both Crow and elliptic instabilities are retrieved in the optimal response and in the energy gain curves. Even if very slow, the viscous diffusion is found to increase the gain of the antisymmetric elliptic perturbation compared to the symmetric one. When the fluid is stratified, peaks at small wavenumber and at wavenumber of the order of the vortex core size are found for all Froude numbers with optimal responses strongly resembling, respectively, the Crow and the elliptic modes with optimal gains corresponding to mean growth rates larger than in the homogeneous case for both modes. However, as the strength of stratification increases (Froude numbers smaller than 2), optimal perturbations start departing from their homogeneous counterpart with large perturbation in the wake of the dipole associated with density effects.