diff --git a/docs/release-notes/3362.doc.md b/docs/release-notes/3362.doc.md
new file mode 100644
index 0000000000..1dae77b3e2
--- /dev/null
+++ b/docs/release-notes/3362.doc.md
@@ -0,0 +1 @@
+Improve {func}`~scanpy.external.pp.harmony_integrate` docs {smaller}`D Kühl`
diff --git a/src/scanpy/external/pp/_harmony_integrate.py b/src/scanpy/external/pp/_harmony_integrate.py
index 27c4d2ac8f..1104690d53 100644
--- a/src/scanpy/external/pp/_harmony_integrate.py
+++ b/src/scanpy/external/pp/_harmony_integrate.py
@@ -4,6 +4,7 @@
 
 from __future__ import annotations
 
+from collections.abc import Sequence  # noqa: TCH003
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -19,7 +20,7 @@
 @doctest_needs("harmonypy")
 def harmony_integrate(
     adata: AnnData,
-    key: str,
+    key: str | Sequence[str],
     *,
     basis: str = "X_pca",
     adjusted_basis: str = "X_pca_harmony",
@@ -42,7 +43,9 @@ def harmony_integrate(
         The annotated data matrix.
     key
         The name of the column in ``adata.obs`` that differentiates
-        among experiments/batches.
+        among experiments/batches. To integrate over two or more covariates,
+        you can pass multiple column names as a list. See ``vars_use``
+        parameter of the ``harmonypy`` package for more details.
     basis
         The name of the field in ``adata.obsm`` where the PCA table is
         stored. Defaults to ``'X_pca'``, which is the default for