The linear stability of three-dimensional (3D) vortices in rotating, stratified flows has been studied by analyzing the non-hydrostatic inviscid Boussinesq equations. We have focused on a widely-used model of geophysical and astrophysical vortices, which assumes an axisymmetric Gaussian structure for pressure anomalies in the horizontal and vertical directions. For a range of Rossby number (−0.5 < Ro < 0.5) and Burger number (0.02 < Bu < 2.3) relevant to observed long-lived vortices, the growth rate and spatial structure of the most unstable eigenmodes have been numerically calculated and presented as a function of Ro − Bu. We have found neutrally-stable vortices only over a small region of the Ro − Bu parameter space: cyclones with Ro ∼ 0.02 − 0.05 and Bu ∼ 0.85 − 0.95. However, we have also found that anticyclones in general have slower growth rates compared to cyclones. In particular, the growth rate of the most unstable eigenmode for anticyclones in a large region of the parameter space (e.g., Ro < 0 and 0.5 Bu 1.3) is slower than 50 turn-around times of the vortex (which often corresponds to several years for ocean eddies). For cyclones, the region with such slow growth rates is confined to 0 < Ro < 0.1 and 0.5 Bu 1.3. While most calculations have been done for f /N = 0.1 (where f andN are the Coriolis and background Brunt-Väisälä frequencies), we have numerically verified and explained analytically, using non-dimensionalized equations, the insensitivity of the results to reducing f /N to the more ocean-relevant value of 0.01. The results of our stability analysis of Gaussian vortices both support and contradict findings of earlier studies with QG or multi-layer models or with other families of vortices. The results of this paper provide a steppingstone to study the more complicated problems of the stability of geophysical (e.g., those in the atmospheres of giant planets) and astrophysical vortices (in accretion disks).