diff --git a/jdaviz/configs/cubeviz/plugins/tests/test_regions.py b/jdaviz/configs/cubeviz/plugins/tests/test_regions.py
index d2a9bcdc6b..7caf95c0a2 100644
--- a/jdaviz/configs/cubeviz/plugins/tests/test_regions.py
+++ b/jdaviz/configs/cubeviz/plugins/tests/test_regions.py
@@ -42,7 +42,7 @@ def test_regions_pixel(self):
             [my_reg], return_bad_regions=True)
         assert len(bad_regions) == 0
         self.verify_region_loaded('Subset 1', count=1)
-        assert len(self.cubeviz.get_interactive_regions()) == 1
+        assert len(self.cubeviz.plugins['Subset Tools']._obj.get_subsets_as_regions()) == 1
 
     def test_regions_sky_has_wcs(self):
         sky = SkyCoord(205.4397, 27.0035, unit='deg')
@@ -63,10 +63,14 @@ def test_spatial_spectral_mix(self):
                                                                           4.896 * unit))
         self.cubeviz.app.session.edit_subset_mode.edit_subset = None
 
-        # Get interactive spatial regions only.
-        spatial_subsets = self.cubeviz.get_interactive_regions()
-        assert list(spatial_subsets.keys()) == ['Subset 1'], spatial_subsets
-        assert isinstance(spatial_subsets['Subset 1'], EllipsePixelRegion)
+        # Get spatial regions only.
+        st = self.cubeviz.plugins['Subset Tools']._obj
+        spatial_subsets_as_regions = st.get_subsets_as_regions(region_type='spatial')
+        assert list(spatial_subsets_as_regions.keys()) == ['Subset 1'], spatial_subsets_as_regions
+        assert isinstance(spatial_subsets_as_regions['Subset 1'], EllipsePixelRegion)
+        # ensure agreement between app.get_subsets and subset_tools.get_subsets_as_regions
+        ss = self.cubeviz.app.get_subsets()
+        assert ss['Subset 1'][0]['region'] == spatial_subsets_as_regions['Subset 1']
 
         # NOTE: This does not test that spectrum from Subset is actually scientifically accurate.
         # Get spectral regions only.
diff --git a/jdaviz/configs/default/plugins/subset_tools/subset_tools.py b/jdaviz/configs/default/plugins/subset_tools/subset_tools.py
index bdd38bcafa..0707a4c236 100644
--- a/jdaviz/configs/default/plugins/subset_tools/subset_tools.py
+++ b/jdaviz/configs/default/plugins/subset_tools/subset_tools.py
@@ -9,7 +9,7 @@
                                         ReplaceMode, XorMode, NewMode)
 from glue.core.roi import CircularROI, CircularAnnulusROI, EllipticalROI, RectangularROI
 from glue.core.subset import (RoiSubsetState, RangeSubsetState, CompositeSubsetState,
-                              MaskSubsetState)
+                              MaskSubsetState, Subset)
 from glue.icons import icon_path
 from glue_jupyter.widgets.subset_mode_vuetify import SelectionModeMenu
 from glue_jupyter.common.toolbar_vuetify import read_icon
@@ -165,12 +165,111 @@ def __init__(self, *args, **kwargs):
 
     @property
     def user_api(self):
-        expose = ['subset', 'combination_mode',
-                  'recenter_dataset', 'recenter',
-                  'get_center', 'set_center',
-                  'import_region']
+        expose = ['subset', 'combination_mode', 'recenter_dataset', 'recente',
+                  'get_center', 'set_center', 'import_region', 'get_subsets_as_regions']
         return PluginUserApi(self, expose)
 
+    def get_subsets_as_regions(self, region_type=None, list_of_subset_labels=None,
+                               use_display_units=False):
+        """
+        Return spatial and/or spectral subsets of ``region_type`` (spatial or
+        spectral, default both) as ``regions`` or ``SpectralRegions`` objects,
+        respectivley.
+
+        Parameters
+        ----------
+        region_type : str or None, optional
+            Specifies the type of subsets to retrieve. Options are:
+            - ``spatial``: Retrieve only spatial subsets.
+            - ``spectral``: Retrieve only spectral subsets.
+            - ``None`` (default): Retrieve both spatial and spectral subsets,
+              determined by the current configuration.
+
+        list_of_subset_labels : list of str or None, optional
+            If specified, only subsets matching these labels will be included.
+            If not specified, all subsets matching the `region_type` will be
+            returned.
+
+        use_display_units : bool, optional
+            (For spectral subsets) If False (default), subsets are returned in
+            the native data unit. If True, subsets are returned in the spectral
+            axis display unit set in the  Unit Conversion plugin.
+
+        Returns
+        -------
+        regions : dict
+            A dictionary mapping subset labels to their respective ``regions``
+            objects (for spatial regions) or ``SpectralRegions`` objects
+            (for spectral regions).
+        """
+
+        if region_type is not None:
+            region_type = region_type.lower()
+            if region_type not in ['spectral', 'spatial']:
+                raise ValueError("`region_type` must be 'spectral' or 'spatial'.")
+            if ((self.config == 'imviz' and region_type == 'spectral') or
+               (self.config == 'specviz' and region_type == 'spatial')):
+                raise ValueError(f"No {region_type} subests in {self.config}.")
+            region_type = [region_type]
+
+        else:  # determine which subset types should be returned by config, if type not specified
+            if self.config == 'imviz':
+                region_type = ['spatial']
+            elif self.config == 'specviz':
+                region_type = ['spectral']
+            else:
+                region_type = ['spatial', 'spectral']
+
+        regions = {}
+
+        if 'spatial' in region_type:
+            from glue_astronomy.translators.regions import roi_subset_state_to_region
+
+            failed_regs = set()
+            to_sky = self.app._align_by == 'wcs'
+
+            # Subset is global, so we just use default viewer.
+            for lyr in self.app._jdaviz_helper.default_viewer._obj.layers:
+                if (not hasattr(lyr, 'layer') or not isinstance(lyr.layer, Subset)
+                        or lyr.layer.ndim not in (2, 3)):
+                    continue
+
+                subset_data = lyr.layer
+                subset_label = subset_data.label
+
+                # TODO: Remove this when Jdaviz support round-tripping, see
+                # https://github.com/spacetelescope/jdaviz/pull/721
+                if not subset_label.startswith('Subset'):
+                    continue
+                try:
+                    if self.app.config == "imviz" and to_sky:
+                        region = roi_subset_state_to_region(subset_data.subset_state,
+                                                            to_sky=to_sky)
+                    else:
+                        region = subset_data.data.get_selection_definition(
+                            subset_id=subset_label, format='astropy-regions')
+                except (NotImplementedError, ValueError):
+                    failed_regs.add(subset_label)
+                else:
+                    regions[subset_label] = region
+
+            if len(failed_regs) > 0:
+                self.app.hub.broadcast(SnackbarMessage(
+                    f"Regions skipped: {', '.join(sorted(failed_regs))}",
+                    color="warning", timeout=8000, sender=self.app))
+
+        if 'spectral' in region_type:
+            spec_regions = self.app.get_subsets(spectral_only=True,
+                                                use_display_units=use_display_units)
+            if spec_regions:
+                regions.update(spec_regions)
+
+        # filter by list_of_subset_labels
+        if list_of_subset_labels is not None:
+            regions = {key: regions[key] for key in list_of_subset_labels}
+
+        return regions
+
     def _on_link_update(self, *args):
         """When linking is changed pixels<>wcs, change display units of the
         subset plugin from pixel (for pixel linking) to sky (for WCS linking).
diff --git a/jdaviz/configs/default/plugins/subset_tools/tests/test_subset_tools.py b/jdaviz/configs/default/plugins/subset_tools/tests/test_subset_tools.py
index 01c75bd658..2672767a4c 100644
--- a/jdaviz/configs/default/plugins/subset_tools/tests/test_subset_tools.py
+++ b/jdaviz/configs/default/plugins/subset_tools/tests/test_subset_tools.py
@@ -4,6 +4,7 @@
 import numpy as np
 from astropy.nddata import NDData
 import astropy.units as u
+from regions import CirclePixelRegion, PixCoord
 from specutils import SpectralRegion
 from glue.core.roi import EllipticalROI, CircularROI, CircularAnnulusROI, RectangularROI
 from glue.core.edit_subset_mode import ReplaceMode, OrMode
@@ -226,3 +227,28 @@ def test_import_spectral_regions_file(cubeviz_helper, spectrum1d_cube, tmp_path)
 
     with pytest.raises(ValueError, match='\'test\' not one of'):
         plg.combination_mode = 'test'
+
+
+def test_get_subsets_as_regions(cubeviz_helper, spectrum1d_cube):
+
+    cubeviz_helper.load_data(spectrum1d_cube)
+    plg = cubeviz_helper.plugins['Subset Tools']
+
+    # load one spectral region, which will become 'Subset 1'
+    plg.import_region(SpectralRegion(1* u.um, 2 * u.um))
+
+    # load one spatial region, which will become 'Subset 2'
+    spatial_reg = CirclePixelRegion(center=PixCoord(x=2, y=2), radius=2)
+    plg.import_region(spatial_reg)
+
+    # call get_subsets_as_regions, which by default for cubeviz will return both
+    # spatial and spectral regoins
+    all_regions = plg.get_subsets_as_regions()
+    # assert len(all_regions) == 2
+
+    # now specify region type
+    spatial_regions = plg.get_subsets_as_regions(region_type='spatial')
+    assert len(spatial_regions) == 1
+
+    spectral_regions = plg.get_subsets_as_regions(region_type='spectral')
+    assert len(spectral_regions) == 1
diff --git a/jdaviz/configs/imviz/tests/test_linking.py b/jdaviz/configs/imviz/tests/test_linking.py
index 9e0b98e152..81a597faec 100644
--- a/jdaviz/configs/imviz/tests/test_linking.py
+++ b/jdaviz/configs/imviz/tests/test_linking.py
@@ -120,7 +120,17 @@ def test_wcslink_affine_with_extras(self):
 
         # Ensure subsets are still there.
         all_labels = [layer.layer.label for layer in self.viewer.state.layers]
-        assert sorted(self.imviz.get_interactive_regions()) == ['Subset 1', 'Subset 2']
+        # Retrieved subsets as sky regions from Subset plugin, and ensure they
+        # match what was loaded and that they are in sky coordinates.
+        subset_as_regions = self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions()
+        assert sorted(subset_as_regions) == ['Subset 1', 'Subset 2']
+        assert_allclose(subset_as_regions['Subset 1'].center.ra.deg, 337.519449)
+        assert_allclose(subset_as_regions['Subset 2'].center.ra.deg, 337.518498)
+        # ensure agreement between app.get_subsets and subset_tools.get_subsets_as_regions
+        ss = self.imviz.app.get_subsets(include_sky_region=True)
+        assert ss['Subset 1'][0]['sky_region'] == subset_as_regions['Subset 1']
+        assert ss['Subset 2'][0]['sky_region'] == subset_as_regions['Subset 2']
+
         assert 'MaskedSubset 1' in all_labels
         assert 'MaskedSubset 2' in all_labels
 
diff --git a/jdaviz/configs/imviz/tests/test_parser.py b/jdaviz/configs/imviz/tests/test_parser.py
index 9289d83c03..d7eaa772ef 100644
--- a/jdaviz/configs/imviz/tests/test_parser.py
+++ b/jdaviz/configs/imviz/tests/test_parser.py
@@ -9,7 +9,7 @@
 from astropy.wcs import WCS
 from gwcs import WCS as GWCS
 from numpy.testing import assert_allclose, assert_array_equal
-from regions import CirclePixelRegion, RectanglePixelRegion
+from regions import CirclePixelRegion, PixCoord, RectanglePixelRegion
 from skimage.io import imsave
 from stdatamodels import asdf_in_fits
 
@@ -259,10 +259,17 @@ def test_parse_jwst_nircam_level2(self, imviz_helper):
         imviz_helper._apply_interactive_region('bqplot:rectangle',
                                                (982, 1088),
                                                (1008, 1077))  # Background
-        subsets = imviz_helper.get_interactive_regions()
+        subsets = imviz_helper.plugins['Subset Tools']._obj.get_subsets_as_regions()
         assert list(subsets.keys()) == ['Subset 1', 'Subset 2'], subsets
+        # check that retrieved subsets-as-regions from subset plugin match what was loaded.
         assert isinstance(subsets['Subset 1'], CirclePixelRegion)
+        assert subsets['Subset 1'].center == PixCoord(970.95, 1116.05)
         assert isinstance(subsets['Subset 2'], RectanglePixelRegion)
+        assert subsets['Subset 2'].center == PixCoord(995.0, 1082.5)
+        # ensure agreement between app.get_subsets and subset_tools.get_subsets_as_regions
+        ss = imviz_helper.app.get_subsets()
+        assert ss['Subset 1'][0]['region'] == subsets['Subset 1']
+        assert ss['Subset 2'][0]['region'] == subsets['Subset 2']
 
         # Test simple aperture photometry plugin.
         phot_plugin = imviz_helper.app.get_tray_item_from_name('imviz-aper-phot-simple')
diff --git a/jdaviz/configs/imviz/tests/test_regions.py b/jdaviz/configs/imviz/tests/test_regions.py
index 4d8c53093a..4ea606093b 100644
--- a/jdaviz/configs/imviz/tests/test_regions.py
+++ b/jdaviz/configs/imviz/tests/test_regions.py
@@ -70,7 +70,7 @@ def test_regions_invalid(self):
 
         # Make sure nothing is actually loaded
         self.verify_region_loaded('MaskedSubset 1', count=0)
-        assert self.imviz.get_interactive_regions() == {}
+        assert self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions() == {}
 
     def test_regions_fully_out_of_bounds(self):
         """Glue ROI will not error when out of bounds."""
@@ -78,7 +78,7 @@ def test_regions_fully_out_of_bounds(self):
         bad_regions = self.subset_plugin.import_region([my_reg], return_bad_regions=True)
         assert len(bad_regions) == 0
         self.verify_region_loaded('Subset 1')
-        assert len(self.imviz.get_interactive_regions()) == 1
+        assert len(self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions()) == 1
 
     def test_regions_mask(self):
         mask = np.zeros((10, 10), dtype=np.bool_)
@@ -86,13 +86,13 @@ def test_regions_mask(self):
         bad_regions = self.subset_plugin.import_region([mask], return_bad_regions=True)
         assert len(bad_regions) == 0
         self.verify_region_loaded('MaskedSubset 1')
-        assert self.imviz.get_interactive_regions() == {}
+        assert self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions() == {}
 
         mask[1, 1] = True
         bad_regions = self.subset_plugin.import_region([mask], return_bad_regions=True)
         assert len(bad_regions) == 0
         self.verify_region_loaded('MaskedSubset 2')
-        assert self.imviz.get_interactive_regions() == {}
+        assert self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions() == {}
 
         # Also test deletion by label here.
         self.imviz._delete_region('MaskedSubset 1')
@@ -103,7 +103,7 @@ def test_regions_mask(self):
         bad_regions = self.subset_plugin.import_region([mask], return_bad_regions=True)
         assert len(bad_regions) == 0
         self.verify_region_loaded('MaskedSubset 3')
-        assert self.imviz.get_interactive_regions() == {}
+        assert self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions() == {}
 
         # Deletion of non-existent label is silent no-op.
         self.imviz._delete_region('foo')
@@ -115,7 +115,7 @@ def test_regions_pixel(self):
                                                        combination_mode='new')
         assert len(bad_regions) == 0
         self.verify_region_loaded('Subset 1')
-        assert len(self.imviz.get_interactive_regions()) == 1
+        assert len(self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions()) == 1
 
     def test_regions_sky_has_wcs(self):
         # Mimic interactive region (before)
@@ -141,7 +141,7 @@ def test_regions_sky_has_wcs(self):
         # Check interactive regions. We do not check if the translation is correct,
         # that check hopefully is already done in glue-astronomy.
         # Apparently, static region ate up one number...
-        subsets = self.imviz.get_interactive_regions()
+        subsets = self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions()
         assert list(subsets.keys()) == ['Subset 1', 'Subset 2', 'Subset 3', 'Subset 5', 'Subset 6'], subsets  # noqa: E501
         assert isinstance(subsets['Subset 1'], CirclePixelRegion)
         assert isinstance(subsets['Subset 2'], CirclePixelRegion)
@@ -161,7 +161,7 @@ def test_regions_annulus_from_load_data(self):
         # Test data is set up such that 1 pixel is 1 arcsec.
         subset_radii = {"Subset 1": [0.5, 1], "Subset 2": [1, 3]}
 
-        subsets = self.imviz.get_interactive_regions()
+        subsets = self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions()
         subset_names = sorted(subsets.keys())
         assert subset_names == ['Subset 1', 'Subset 2']
         for n in subset_names:
@@ -174,7 +174,7 @@ def test_photutils_pixel(self):
         bad_regions = self.subset_plugin.import_region([my_aper], return_bad_regions=True)
         assert len(bad_regions) == 0
         self.verify_region_loaded('Subset 1')
-        assert len(self.imviz.get_interactive_regions()) == 1
+        assert len(self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions()) == 1
 
     def test_photutils_sky_has_wcs(self):
         sky = SkyCoord(ra=337.5202808, dec=-20.833333059999998, unit='deg')
@@ -182,7 +182,7 @@ def test_photutils_sky_has_wcs(self):
         bad_regions = self.subset_plugin.import_region([my_aper_sky], return_bad_regions=True)
         assert len(bad_regions) == 0
         self.verify_region_loaded('Subset 1')
-        assert len(self.imviz.get_interactive_regions()) == 1
+        assert len(self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions()) == 1
 
 
 class TestLoadRegionsFromFile(BaseRegionHandler):
@@ -204,7 +204,7 @@ def test_ds9_load_all(self, imviz_helper):
         assert len(bad_regions) == 1
 
         # Will load 8/9 and 7 of that become ROIs.
-        subsets = imviz_helper.get_interactive_regions()
+        subsets = imviz_helper.plugins['Subset Tools']._obj.get_subsets_as_regions()
         assert list(subsets.keys()) == ['Subset 1', 'Subset 2', 'Subset 3',
                                         'Subset 4', 'Subset 5', 'Subset 6', 'Subset 7'], subsets
 
@@ -217,7 +217,7 @@ def test_ds9_load_two_good(self, imviz_helper):
         bad_regions = imviz_helper.plugins['Subset Tools'].import_region(
             self.region_file, combination_mode='new', max_num_regions=2, return_bad_regions=True)
         assert len(bad_regions) == 0
-        subsets = imviz_helper.get_interactive_regions()
+        subsets = imviz_helper.plugins['Subset Tools']._obj.get_subsets_as_regions()
         assert list(subsets.keys()) == ['Subset 1', 'Subset 2'], subsets
         self.verify_region_loaded('MaskedSubset 1', count=0)
 
@@ -227,7 +227,7 @@ def test_ds9_load_one_bad(self, imviz_helper):
         bad_regions = imviz_helper.plugins['Subset Tools'].import_region(
             self.raw_regions[6], return_bad_regions=True)
         assert len(bad_regions) == 1
-        assert imviz_helper.get_interactive_regions() == {}
+        assert imviz_helper.plugins['Subset Tools']._obj.get_subsets_as_regions() == {}
         self.verify_region_loaded('MaskedSubset 1', count=0)
 
     def test_ds9_load_one_good_one_bad(self, imviz_helper):
@@ -237,12 +237,12 @@ def test_ds9_load_one_good_one_bad(self, imviz_helper):
             [self.raw_regions[3], self.raw_regions[6]], return_bad_regions=True)
         assert len(bad_regions) == 1
 
-        subsets = imviz_helper.get_interactive_regions()
+        subsets = imviz_helper.plugins['Subset Tools']._obj.get_subsets_as_regions()
         assert list(subsets.keys()) == ['Subset 1'], subsets
         self.verify_region_loaded('MaskedSubset 1', count=0)
 
 
-class TestGetInteractiveRegions(BaseImviz_WCS_NoWCS):
+class TestGetRegions(BaseImviz_WCS_NoWCS):
     def test_annulus(self):
         # Outer circle
         self.imviz._apply_interactive_region('bqplot:truecircle', (0, 0), (9, 9))
@@ -250,22 +250,30 @@ def test_annulus(self):
         self.imviz._apply_interactive_region('bqplot:truecircle', (2, 2), (7, 7))
 
         # At this point, there should be two normal circles.
-        subsets = self.imviz.get_interactive_regions()
+        subsets = self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions()
         assert list(subsets.keys()) == ['Subset 1', 'Subset 2'], subsets
         assert isinstance(subsets['Subset 1'], CirclePixelRegion)
         assert isinstance(subsets['Subset 2'], CirclePixelRegion)
+        assert subsets['Subset 1'].center == PixCoord(4.5, 4.5)
+        assert subsets['Subset 2'].center == PixCoord(4.5, 4.5)
+        # ensure agreement between app.get_subsets and subset_tools.get_subsets_as_regions
+        ss = self.imviz.app.get_subsets()
+        assert ss['Subset 1'][0]['region'] == subsets['Subset 1']
+        assert ss['Subset 2'][0]['region'] == subsets['Subset 2']
 
         # Create a third subset that is an annulus.
         subset_groups = self.imviz.app.data_collection.subset_groups
         new_subset = subset_groups[0].subset_state & ~subset_groups[1].subset_state
         self.viewer.apply_subset_state(new_subset)
 
-        subsets = self.imviz.get_interactive_regions()
+        subsets = self.imviz.plugins['Subset Tools']._obj.get_subsets_as_regions()
         assert len(self.imviz.app.data_collection.subset_groups) == 3
         assert list(subsets.keys()) == ['Subset 1', 'Subset 2', 'Subset 3'], subsets
         assert isinstance(subsets['Subset 1'], CirclePixelRegion)
         assert isinstance(subsets['Subset 2'], CirclePixelRegion)
         assert isinstance(subsets['Subset 3'], CircleAnnulusPixelRegion)
+        # and check the new retrieved region
+        assert subsets['Subset 3'].center == PixCoord(4.5, 4.5)
 
         # Clear the regions for next test.
         self.imviz._delete_all_regions()
diff --git a/jdaviz/configs/specviz/helper.py b/jdaviz/configs/specviz/helper.py
index c13ac67175..627379ef27 100644
--- a/jdaviz/configs/specviz/helper.py
+++ b/jdaviz/configs/specviz/helper.py
@@ -1,6 +1,7 @@
 import warnings
 
 from astropy import units as u
+from astropy.utils.decorators import deprecated
 from regions.core.core import Region
 from glue.core.subset_group import GroupedSubset
 from specutils import SpectralRegion, Spectrum1D
@@ -147,6 +148,7 @@ def get_spectra(self, data_label=None, spectral_subset=None, apply_slider_redshi
 
             return output_spectra
 
+    @deprecated(since="4.1", alternative="subset_tools.get_subsets_as_regions")
     def get_spectral_regions(self, use_display_units=False):
         """
         A simple wrapper around the app-level call to retrieve only spectral
diff --git a/jdaviz/configs/specviz/tests/test_helper.py b/jdaviz/configs/specviz/tests/test_helper.py
index a3d943e430..d70547630e 100644
--- a/jdaviz/configs/specviz/tests/test_helper.py
+++ b/jdaviz/configs/specviz/tests/test_helper.py
@@ -134,7 +134,7 @@ def test_get_spectral_regions_two(self):
                                  7800*self.spec.spectral_axis.unit))
         self.spec_app.plugins['Subset Tools'].import_region(subset, combination_mode='or')
 
-        spec_region = self.spec_app.get_spectral_regions()
+        spec_region = self.spec_app.plugins['Subset Tools']._obj.get_subsets_as_regions()
 
         assert len(spec_region['Subset 1'].subregions) == 2
 
@@ -147,7 +147,7 @@ def test_get_spectral_regions_three(self):
                                  7800*self.spec.spectral_axis.unit))
         self.spec_app.plugins['Subset Tools'].import_region(subset, combination_mode='or')
 
-        spec_region = self.spec_app.get_spectral_regions()
+        spec_region = self.spec_app.plugins['Subset Tools']._obj.get_subsets_as_regions()
 
         assert len(spec_region['Subset 1'].subregions) == 3
         # Assert correct values for test with 3 subregions
@@ -173,7 +173,7 @@ def test_get_spectral_regions_does_not_raise_value_error(self):
                                  6*self.spec.spectral_axis.unit))
         self.spec_app.plugins['Subset Tools'].import_region(subset, combination_mode='or')
 
-        spec_region = self.spec_app.get_spectral_regions()
+        spec_region = self.spec_app.plugins['Subset Tools']._obj.get_subsets_as_regions()
         assert_quantity_allclose(spec_region['Subset 1'].subregions[0][0].value,
                                  1, atol=1e-5)
         assert_quantity_allclose(spec_region['Subset 1'].subregions[0][1].value,
@@ -194,7 +194,7 @@ def test_get_spectral_regions_composite_region(self):
         self.spec_app.plugins['Subset Tools'].import_region(
             subset, combination_mode=['new', 'andnot', 'and'])
 
-        spec_region = self.spec_app.get_spectral_regions()
+        spec_region = self.spec_app.plugins['Subset Tools']._obj.get_subsets_as_regions()
 
         assert len(spec_region['Subset 1'].subregions) == 1
         # Assert correct values for test with 3 subregions
@@ -213,7 +213,7 @@ def test_get_spectral_regions_composite_region_multiple_and_nots(self):
         self.spec_app.plugins['Subset Tools'].import_region(
             subset, combination_mode=['new', 'andnot', 'andnot'])
 
-        spec_region = self.spec_app.get_spectral_regions()
+        spec_region = self.spec_app.plugins['Subset Tools']._obj.get_subsets_as_regions()
 
         assert len(spec_region['Subset 1'].subregions) == 3
         # Assert correct values for test with 3 subregions
@@ -274,7 +274,7 @@ def test_get_spectral_regions_unit(specviz_helper, spectrum1d):
                             7000 * spectrum1d.spectral_axis.unit)
     specviz_helper.plugins['Subset Tools'].import_region(subset)
 
-    subsets = specviz_helper.get_spectral_regions()
+    subsets = specviz_helper.plugins['Subset Tools']._obj.get_subsets_as_regions()
     reg = subsets.get('Subset 1')
 
     assert spectrum1d.spectral_axis.unit == reg.lower.unit
@@ -317,13 +317,13 @@ def test_get_spectral_regions_unit_conversion(specviz_helper, spectrum1d):
     specviz_helper.plugins['Subset Tools'].import_region(subset)
 
     # Retrieve the Subset
-    subsets = specviz_helper.get_spectral_regions(use_display_units=False)
-    reg = subsets.get('Subset 1')
+    ss = specviz_helper.plugins['Subset Tools']._obj.get_subsets_as_regions(use_display_units=False)
+    reg = ss.get('Subset 1')
     assert reg.lower.unit == u.Angstrom
     assert reg.upper.unit == u.Angstrom
 
-    subsets = specviz_helper.get_spectral_regions(use_display_units=True)
-    reg = subsets.get('Subset 1')
+    ss = specviz_helper.plugins['Subset Tools']._obj.get_subsets_as_regions(use_display_units=True)
+    reg = ss.get('Subset 1')
     assert reg.lower.unit == u.um
     assert reg.upper.unit == u.um
 
diff --git a/jdaviz/core/helpers.py b/jdaviz/core/helpers.py
index 98583b5617..df42464aae 100644
--- a/jdaviz/core/helpers.py
+++ b/jdaviz/core/helpers.py
@@ -19,8 +19,9 @@
 from glue.config import data_translator
 from ipywidgets.widgets import widget_serialization
 
-import astropy.units as u
 from astropy.nddata import NDDataArray, CCDData, StdDevUncertainty
+import astropy.units as u
+from astropy.utils.decorators import deprecated
 from regions.core.core import Region
 from specutils import Spectrum1D, SpectralRegion
 
@@ -878,8 +879,10 @@ def load_regions(self, regions, max_num_regions=None, refdata_label=None,
         if return_bad_regions:
             return bad_regions
 
+    @deprecated(since="4.1", alternative="subset_tools.get_subsets_as_regions")
     def get_interactive_regions(self):
-        """Return spatial regions that can be interacted with in the viewer.
+        """
+        Return spatial regions that can be interacted with in the viewer.
         This does not return masked regions added via :meth:`load_regions`.
 
         Unsupported region shapes will be skipped. When that happens,