The variability of surface and deep layers in the southern Gulf of Mexico and their predictability and stochastic origin are studied. Considering separated and coupled layers analyses, the most important variability modes were estimated via Empirical Orthogonal Functions using daily isopycnic layer-thickness anomalies from a 21-year free-running simulation of the Gulf hydrodynamics performed with the HYbrid Coordinate Ocean Model. There is a separation between the principal and higher-order coupled variability. The deep layer strongly determines the variability throughout the water column for the principal coupled variability: the timescales and long-term persistence are mainly associated with deep dynamics. Higher-order coupled variability has no clear association with surface or deep dynamics. Deep dynamics is likely to influence the subsequent evolution of surface dynamics; however, an evident causality relationship between them was not found. No vertical correspondence between surface and deep isopycnal fluctuations was found. The principal coupled variability mode is described by a surface region in the southwest where the Campeche Gyre occurs and a deep region in the center of the basin extending to the north. The predictability was estimated through the decorrelation times of the variability modes. The predictability of deep variability is three times that of surface variability, with 30.5-month predictability for the principal deep mode. Layer coupling evinced the role of the deep ocean in generating long-term variability by extending the predictability of the principal surface mode 2.6-fold, from 10.6 to 27.2 months. Strong evidence is provided for the stochastic origin of the principal variability, suggesting it can be described using linear dynamics in terms of a fast and a slow component.