diff --git a/.cruft.json b/.cruft.json
index 1bba414c..cfbef221 100644
--- a/.cruft.json
+++ b/.cruft.json
@@ -11,7 +11,7 @@
       "project_slug": "xscen",
       "project_short_description": "A climate change scenario-building analysis framework, built with xclim/xarray.",
       "pypi_username": "RondeauG",
-      "version": "0.10.1-dev.4",
+      "version": "0.10.1-dev.5",
       "use_pytest": "y",
       "use_black": "y",
       "use_conda": "y",
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 03a4f0f3..b642f67e 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -17,10 +17,16 @@ New features and enhancements
 Breaking changes
 ^^^^^^^^^^^^^^^^
 * ``xs.get_warming_level`` has been renamed to ``xs.get_period_from_warming_level``. Its argument `return_horizon` was reversed and renamed `return_central_year` (:pull:`474`).
+* Removed support for the deprecated `xclim` function `change_significance` in `ensemble_stats`. (:pull:`482`).
 
 Bug fixes
 ^^^^^^^^^
 * ``xs.io.save_to_table`` now correctly handles the case where the input is a `DataArray` or a `Dataset` with a single variable. (:pull:`473`).
+* Fixed a bug in ``xs.utils.change_units`` where the original dataset was also getting modified. (:pull:`482`).
+
+Internal changes
+^^^^^^^^^^^^^^^^
+* Bumped the version of `xclim` to 0.53.2. (:pull:`482`).
 
 v0.10.0 (2024-09-30)
 --------------------
diff --git a/docs/notebooks/4_ensembles.ipynb b/docs/notebooks/4_ensembles.ipynb
index cf6a2576..dad0407f 100644
--- a/docs/notebooks/4_ensembles.ipynb
+++ b/docs/notebooks/4_ensembles.ipynb
@@ -53,9 +53,7 @@
     "}\n",
     "\n",
     "for d in datasets:\n",
-    "    ds = open_dataset(datasets[d], branch=\"v2023.12.14\").isel(\n",
-    "        lon=slice(0, 4), lat=slice(0, 4)\n",
-    "    )\n",
+    "    ds = open_dataset(datasets[d]).isel(lon=slice(0, 4), lat=slice(0, 4))\n",
     "    ds = xs.climatological_op(\n",
     "        ds,\n",
     "        op=\"mean\",\n",
@@ -66,6 +64,15 @@
     "    datasets[d] = xs.compute_deltas(ds, reference_horizon=\"1981-2010\")"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "datasets"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -263,7 +270,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.2"
+   "version": "3.12.7"
   }
  },
  "nbformat": 4,
diff --git a/environment-dev.yml b/environment-dev.yml
index 9d4434da..5b06e841 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -30,7 +30,7 @@ dependencies:
   - sparse
   - toolz
   - xarray >=2023.11.0, !=2024.6.0
-  - xclim >=0.52.2, <0.53
+  - xclim >=0.53.2, <0.54
   - xesmf >=0.7
   - zarr >=2.13
   # Opt
@@ -40,13 +40,13 @@ dependencies:
   - babel
   - black ==24.8.0
   - blackdoc ==0.3.9
-  - bump-my-version >=0.25.1
+  - bump-my-version >=0.26.8
   - coverage>=7.5.0
   - coveralls>=4.0.1
-  - flake8 >=6.1.0
+  - flake8 >=7.1.0
   - flake8-rst-docstrings>=0.3.0
   - ipykernel
-  - ipython
+  - ipython >=8.5.0
   - isort ==5.13.2
   - jupyter_client
   - nbsphinx
diff --git a/environment.yml b/environment.yml
index e5336ab8..3e4232f3 100644
--- a/environment.yml
+++ b/environment.yml
@@ -30,7 +30,7 @@ dependencies:
   - sparse
   - toolz
   - xarray >=2023.11.0, !=2024.6.0
-  - xclim >=0.52.2, <0.53
+  - xclim >=0.53.2, <0.54
   - xesmf >=0.7
   - zarr >=2.13
   # To install from source
@@ -39,4 +39,4 @@ dependencies:
   # Opt
   - nc-time-axis >=1.3.1
   - pyarrow >=10.0.1
-  - pip
+  - pip >=24.2.0
diff --git a/pyproject.toml b/pyproject.toml
index d8771658..21e1d4b9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -65,7 +65,7 @@ dependencies = [
   "sparse",
   "toolz",
   "xarray >=2023.11.0, !=2024.6.0",
-  "xclim >=0.52.2, <0.53",
+  "xclim >=0.53.2, <0.54",
   "zarr >=2.13"
 ]
 
@@ -84,7 +84,7 @@ dev = [
   "isort ==5.13.2",
   "mypy",
   "numpydoc >=1.8.0",
-  "pooch",
+  "pooch >=1.8.0",
   "pre-commit >=3.3.2",
   "pytest-cov >=5.0.0",
   "pytest >=8.3.2",
@@ -97,7 +97,7 @@ dev = [
 docs = [
   # Documentation and examples
   "ipykernel",
-  "ipython",
+  "ipython >=8.5.0",
   "jupyter_client",
   "nbsphinx",
   "nbval",
@@ -133,7 +133,7 @@ target-version = [
 ]
 
 [tool.bumpversion]
-current_version = "0.10.1-dev.4"
+current_version = "0.10.1-dev.5"
 commit = true
 commit_args = "--no-verify"
 tag = false
diff --git a/src/xscen/__init__.py b/src/xscen/__init__.py
index 3b2f7704..d9cea8d2 100644
--- a/src/xscen/__init__.py
+++ b/src/xscen/__init__.py
@@ -71,7 +71,7 @@
 
 __author__ = """Gabriel Rondeau-Genesse"""
 __email__ = "rondeau-genesse.gabriel@ouranos.ca"
-__version__ = "0.10.1-dev.4"
+__version__ = "0.10.1-dev.5"
 
 
 def warning_on_one_line(
diff --git a/src/xscen/ensembles.py b/src/xscen/ensembles.py
index 5d9181a4..cc55fb99 100644
--- a/src/xscen/ensembles.py
+++ b/src/xscen/ensembles.py
@@ -108,7 +108,7 @@ def ensemble_stats(  # noqa: C901
     # if input files are .zarr, change the engine automatically
     if isinstance(datasets, list) and isinstance(datasets[0], str | os.PathLike):
         path = Path(datasets[0])
-        if path.suffix == ".zarr":
+        if path.suffix in [".zarr", ".zip"]:
             create_kwargs.setdefault("engine", "zarr")
 
     if not isinstance(datasets, xr.Dataset):
@@ -137,6 +137,7 @@ def ensemble_stats(  # noqa: C901
 
         # Workaround for robustness_categories
         real_stat = None
+        categories_kwargs = {}
         if stat == "robustness_categories":
             real_stat = "robustness_categories"
             stat = "robustness_fractions"
@@ -155,20 +156,10 @@ def ensemble_stats(  # noqa: C901
                     f"Weighting is not supported for '{stat}'. The results may be incorrect."
                 )
 
-        # FIXME: change_significance is deprecated and will be removed in xclim 0.49.
         if stat in [
-            "change_significance",
             "robustness_fractions",
             "robustness_categories",
         ]:
-            # FIXME: This can be removed once change_significance is removed.
-            #  It's here because the 'ref' default was removed for change_significance in xclim 0.47.
-            stats_kwargs.setdefault("ref", None)
-            if (stats_kwargs.get("ref") is not None) and len(statistics_to_compute) > 1:
-                raise ValueError(
-                    f"The input requirements for '{stat}' when 'ref' is specified are not compatible with other statistics."
-                )
-
             # These statistics only work on DataArrays
             for v in ens.data_vars:
                 with xr.set_options(keep_attrs=True):
@@ -188,22 +179,6 @@ def ensemble_stats(  # noqa: C901
                     # Call the function
                     tmp = getattr(ensembles, stat)(ens_v, **stats_kwargs)
 
-                    # Manage the multiple outputs of change_significance
-                    # FIXME: change_significance is deprecated and will be removed in xclim 0.49.
-                    if (
-                        stat == "change_significance"
-                        and stats_kwargs.get("p_vals", False) is False
-                    ):
-                        ens_stats[f"{v}_change_frac"], ens_stats[f"{v}_pos_frac"] = tmp
-                    elif stat == "change_significance" and stats_kwargs.get(
-                        "p_vals", False
-                    ):
-                        (
-                            ens_stats[f"{v}_change_frac"],
-                            ens_stats[f"{v}_pos_frac"],
-                            ens_stats[f"{v}_p_vals"],
-                        ) = tmp
-
                     # Robustness categories
                     if real_stat == "robustness_categories":
                         categories = ensembles.robustness_categories(
diff --git a/src/xscen/utils.py b/src/xscen/utils.py
index 0dcc8871..d68bebec 100644
--- a/src/xscen/utils.py
+++ b/src/xscen/utils.py
@@ -812,20 +812,26 @@ def change_units(ds: xr.Dataset, variables_and_units: dict) -> xr.Dataset:
                 ).dimensionality.get("[time]")
 
                 if time_in_ds == time_in_out:
-                    ds[v] = units.convert_units_to(ds[v], variables_and_units[v])
+                    ds = ds.assign(
+                        {v: units.convert_units_to(ds[v], variables_and_units[v])}
+                    )
                 elif time_in_ds - time_in_out == 1:
                     # ds is an amount
-                    ds[v] = units.amount2rate(ds[v], out_units=variables_and_units[v])
+                    ds = ds.assign(
+                        {v: units.amount2rate(ds[v], out_units=variables_and_units[v])}
+                    )
                 elif time_in_ds - time_in_out == -1:
                     # ds is a rate
-                    ds[v] = units.rate2amount(ds[v], out_units=variables_and_units[v])
+                    ds = ds.assign(
+                        {v: units.rate2amount(ds[v], out_units=variables_and_units[v])}
+                    )
                 else:
                     raise ValueError(
                         f"No known transformation between {ds[v].units} and {variables_and_units[v]} (temporal dimensionality mismatch)."
                     )
             elif (v in ds) and (ds[v].units != variables_and_units[v]):
                 # update unit name if physical units are equal but not their name (ex. degC vs °C)
-                ds[v] = ds[v].assign_attrs(units=variables_and_units[v])
+                ds = ds.assign({v: ds[v].assign_attrs(units=variables_and_units[v])})
 
     return ds
 
diff --git a/tests/test_ensembles.py b/tests/test_ensembles.py
index 6aa596e2..e21a409a 100644
--- a/tests/test_ensembles.py
+++ b/tests/test_ensembles.py
@@ -258,22 +258,6 @@ def test_errors(self):
                 },
             )
 
-        # Error if you try to use a robustness_fractions with a reference dataset, but also specify other statistics
-        with pytest.raises(
-            ValueError, match="The input requirements for 'robustness_fractions'"
-        ):
-            xs.ensemble_stats(
-                ens,
-                statistics={
-                    "robustness_fractions": {
-                        "test": "threshold",
-                        "abs_thresh": 2.5,
-                        "ref": ref,
-                    },
-                    "ensemble_mean_std_max_min": None,
-                },
-            )
-
 
 class TestGenerateWeights:
     @staticmethod