Marine heatwaves (MHWs) are ocean temperature extremes that can occur at any ocean depth. Surface characteristics and physical mechanisms of MHWs have been extensively explored based on satellite observations; however, their subsurface features and drivers remain unclear. This study investigates the characteristics and drivers of subsurface MHWs near the thermocline in the Bay of Bengal (BoB) from 1993 to 2024 using high-resolution ocean reanalysis datasets. Annual subsurface MHW days exhibit a dipole pattern in response to the El Niño-Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD). This dipole structure is characterized by spatiotemporal asymmetries, with a later (earlier) peak and wider (smaller) coverage in the eastern (western) BoB. During the developing and mature phases of pure positive IOD, the anticyclonic mesoscale eddies in the western BoB are favorable for MHW generation in this region, related to the anomalous anticyclonic winds and currents over the BoB. Meanwhile, the equatorial easterly anomalies drive upwelling Kelvin waves that propagate eastward into the eastern and even western BoB, inhibiting MHW formation over extensive areas during the mature and decay phases of pure positive IOD. Differently, during the developing phase of ENSO, its direct impact on subsurface MHWs in the BoB is limited. It is likely to exert indirect forcing by triggering IOD. The direct impacts of ENSO emerge during the mature phase and reach peaks during the decay phase. This study highlights the inconsistent spatial responses of subsurface MHWs to distinct ocean dynamics induced by ENSO and IOD.