diff --git a/docs/references.bib b/docs/references.bib
index 5a155dd05..db0f6e105 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -40,6 +40,17 @@ @article{omnipath
      doi = {10.1038/nmeth.4077}
 }
 
+@article{rustworkx,
+  author = {Treinish, Matthew and Carvalho, Ivan and Tsilimigkounakis, Georgios and Sá, Nahum},
+  title = {rustworkx: A High-Performance Graph Library for Python},
+  journal = {Journal of Open Source Software},
+  volume = {7},
+  number = {79},
+  pages = {3968},
+  year = {2022},
+  doi = {10.21105/joss.03968},
+}
+
 @InProceedings{networkx,
   author = {Aric A. Hagberg and Daniel A. Schult and Pieter J. Swart},
   title = {Exploring Network Structure, Dynamics, and Function using NetworkX},
diff --git a/pyproject.toml b/pyproject.toml
index 2a9967bb3..0baf4dc8a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -64,6 +64,7 @@ dependencies = [
   "pillow>=8",
   "pooch>=1.6",
   "pyyaml>=6",
+  "rustworkx>=0.18",
   "scanpy>=1.9.3",
   "scikit-image>=0.25",
   # due to https://github.com/scikit-image/scikit-image/issues/6850 breaks rescale ufunc
diff --git a/src/squidpy/_docs.py b/src/squidpy/_docs.py
index 5a98af9b9..5ad515b2a 100644
--- a/src/squidpy/_docs.py
+++ b/src/squidpy/_docs.py
@@ -213,7 +213,7 @@ def decorator2(obj: Any) -> Any:
 scalebar_kwargs
     Keyword arguments for :class:`matplotlib_scalebar.ScaleBar`.
 edges_kwargs
-    Keyword arguments for :func:`networkx.drawing.nx_pylab.draw_networkx_edges`.
+    Keyword arguments for :func:`networkx.drawing.nx_pylab.draw_networkx_edges` :cite:`networkx`.
 kwargs
     Keyword arguments for :func:`matplotlib.pyplot.scatter` or :func:`matplotlib.pyplot.imshow`.
 """
diff --git a/src/squidpy/gr/_nhood.py b/src/squidpy/gr/_nhood.py
index a92f8a5ea..64e1a1741 100644
--- a/src/squidpy/gr/_nhood.py
+++ b/src/squidpy/gr/_nhood.py
@@ -6,15 +6,16 @@
 from functools import partial
 from typing import Any, NamedTuple
 
-import networkx as nx
 import numba.types as nt
 import numpy as np
 import pandas as pd
+import rustworkx as rx
 from anndata import AnnData
-from numba import njit
+from numba import njit, prange
 from numpy.typing import NDArray
 from pandas import CategoricalDtype
 from scanpy import logging as logg
+from scipy.sparse import csr_matrix
 from spatialdata import SpatialData
 
 from squidpy._constants._constants import Centrality
@@ -266,7 +267,7 @@ def centrality_scores(
     %(table_key)s
     %(cluster_key)s
     score
-        Centrality measures as described in :mod:`networkx.algorithms.centrality` :cite:`networkx`.
+        Group centrality measures as implemented in ``rustworkx`` :cite:`rustworkx`.
         If `None`, use all the options below. Valid options are:
 
             - `{c.CLOSENESS.s!r}` - measure of how close the group is to other nodes.
@@ -296,7 +297,9 @@ def centrality_scores(
 
     centralities = [Centrality(c) for c in centrality]
 
-    graph = nx.Graph(adata.obsp[connectivity_key])
+    # a rustworkx graph mirrors the undirected connectivity graph for the group closeness/degree
+    # measures; a symmetric, self-loop-free CSR feeds the clustering-coefficient kernel.
+    graph, adj = _build_graph(adata.obsp[connectivity_key])
 
     cat = adata.obs[cluster_key].cat.categories.values
     clusters = adata.obs[cluster_key].values
@@ -304,11 +307,13 @@ def centrality_scores(
     fun_dict = {}
     for c in centralities:
         if c == Centrality.CLOSENESS:
-            fun_dict[c.s] = partial(nx.algorithms.centrality.group_closeness_centrality, graph)
+            fun_dict[c.s] = partial(rx.group_closeness_centrality, graph)
         elif c == Centrality.DEGREE:
-            fun_dict[c.s] = partial(nx.algorithms.centrality.group_degree_centrality, graph)
+            fun_dict[c.s] = partial(rx.group_degree_centrality, graph)
         elif c == Centrality.CLUSTERING:
-            fun_dict[c.s] = partial(nx.algorithms.cluster.average_clustering, graph)
+            # average the per-node clustering coefficients over the group (0 if the group is empty).
+            node_clustering = _local_clustering(adj.indptr, adj.indices, adj.shape[0])
+            fun_dict[c.s] = lambda idx, cc=node_clustering: float(cc[idx].mean()) if len(idx) else 0.0
         else:
             raise NotImplementedError(f"Centrality `{c}` is not yet implemented.")
 
@@ -424,6 +429,68 @@ def _interaction_matrix(
     return output
 
 
+def _build_graph(conn: Any) -> tuple[rx.PyGraph, csr_matrix]:
+    """Build the graph representations used by :func:`centrality_scores`.
+
+    Returns a :class:`rustworkx.PyGraph` mirroring the undirected connectivity graph
+    (used by the group closeness/degree measures) and a symmetric, self-loop-free,
+    index-sorted CSR matrix feeding the clustering-coefficient kernel.
+    """
+    from scipy.sparse import triu
+
+    adj = csr_matrix(conn)
+    # undirected, unweighted, no self-loops: matches ``networkx.Graph(conn)`` topology.
+    adj = (adj + adj.T).tocsr()
+    adj.setdiag(0)
+    adj.eliminate_zeros()
+    adj.sort_indices()  # the clustering kernel merges neighbor lists, which must be sorted.
+
+    n = adj.shape[0]
+    graph = rx.PyGraph(multigraph=False)
+    graph.add_nodes_from(range(n))
+    # the strict upper triangle lists each undirected edge exactly once.
+    rows, cols = triu(adj, k=1).nonzero()
+    graph.add_edges_from_no_data([(int(i), int(j)) for i, j in zip(rows, cols, strict=True)])
+    return graph, adj
+
+
+@njit(parallel=True, cache=True)
+def _local_clustering(indptr: NDArrayA, indices: NDArrayA, n: int) -> NDArrayA:
+    """Local clustering coefficient per node over a symmetric, index-sorted CSR graph.
+
+    ``c_v = 2 * T(v) / (k_v * (k_v - 1))`` where ``T(v)`` is the number of edges among the
+    neighbors of ``v`` and ``k_v`` its degree; ``c_v = 0`` when ``k_v < 2``. Triangles are
+    counted by intersecting sorted neighbor lists, so only existing edges are ever visited
+    (no length-2-path matrix is materialized). Matches :func:`networkx.clustering`.
+    """
+    out = np.zeros(n, dtype=np.float64)
+    for v in prange(n):
+        start = indptr[v]
+        end = indptr[v + 1]
+        k = end - start
+        if k < 2:
+            continue
+        # summing |N(u) ∩ N(v)| over u in N(v) counts each neighbor-neighbor edge twice,
+        # so it already equals 2 * T(v) -> c_v = that sum / (k * (k - 1)).
+        two_triangles = 0
+        for a in range(start, end):
+            u = indices[a]
+            i = start
+            j = indptr[u]
+            u_end = indptr[u + 1]
+            while i < end and j < u_end:
+                if indices[i] == indices[j]:
+                    two_triangles += 1
+                    i += 1
+                    j += 1
+                elif indices[i] < indices[j]:
+                    i += 1
+                else:
+                    j += 1
+        out[v] = two_triangles / (k * (k - 1))
+    return out
+
+
 def _centrality_scores_helper(
     cat: Iterable[Any],
     clusters: Sequence[str],
diff --git a/tests/graph/test_nhood.py b/tests/graph/test_nhood.py
index f5c721d4d..fb19fc220 100644
--- a/tests/graph/test_nhood.py
+++ b/tests/graph/test_nhood.py
@@ -88,6 +88,23 @@ def test_centrality_scores(nhood_data: AnnData):
     assert adata.uns[key]["closeness_centrality"].dtype == np.dtype("float64")
 
 
+def test_centrality_scores_networkx_parity(nhood_data: AnnData):
+    # centrality_scores swapped networkx for rustworkx (+ a numba clustering kernel); pin the
+    # numeric parity of all three group measures against networkx (still a dependency).
+    import networkx as nx
+
+    adata = nhood_data
+    df = centrality_scores(adata, cluster_key=_CK, connectivity_key="spatial", copy=True)
+
+    graph = nx.Graph(adata.obsp["spatial_connectivities"])
+    clusters = adata.obs[_CK].values
+    for cat in df.index:
+        idx = list(np.where(clusters == cat)[0])
+        np.testing.assert_allclose(df.loc[cat, "closeness_centrality"], nx.group_closeness_centrality(graph, idx))
+        np.testing.assert_allclose(df.loc[cat, "degree_centrality"], nx.group_degree_centrality(graph, idx))
+        np.testing.assert_allclose(df.loc[cat, "average_clustering"], nx.average_clustering(graph, idx))
+
+
 @pytest.mark.parametrize("copy", [True, False])
 def test_interaction_matrix_copy(nhood_data: AnnData, copy: bool):
     adata = nhood_data