diff --git a/NEWS.md b/NEWS.md
index 9e3b8f0cae..ff6a7e5578 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 * Some methods related to `get_vector` for `GeneralUnitaryMatrices` now have `AbstractVector` upper bound for coefficients.
 * Minimum Julia version is now 1.10 (the LTS which replaced 1.6)
+* The dependency ManifoldDiff.jl has been upgraded from v0.3 to v0.4, to bring compatibility with DifferentiationInterface.jl.
 
 ## [0.10.7] – 2024-11-16
 
diff --git a/Project.toml b/Project.toml
index 3ab0bb342e..48563914a7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -53,7 +53,7 @@ Graphs = "1.4"
 HybridArrays = "0.4"
 Kronecker = "0.4, 0.5"
 LinearAlgebra = "1.6"
-ManifoldDiff = "0.3.13"
+ManifoldDiff = "0.4.0"
 ManifoldsBase = "0.15.18"
 Markdown = "1.6"
 MatrixEquations = "2.2"
@@ -74,6 +74,7 @@ StatsBase = "0.32, 0.33, 0.34"
 julia = "1.10"
 
 [extras]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
@@ -98,4 +99,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 VisualRegressionTests = "34922c18-7c2a-561c-bac1-01e79b2c4c92"
 
 [targets]
-test = ["Test", "BoundaryValueDiffEq", "Colors", "DiffEqCallbacks", "Distributions", "DoubleFloats", "FiniteDifferences", "Gtk", "HybridArrays", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PythonPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "RecursiveArrayTools"]
+test = ["Test", "ADTypes", "BoundaryValueDiffEq", "Colors", "DiffEqCallbacks", "Distributions", "DoubleFloats", "FiniteDifferences", "Gtk", "HybridArrays", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PythonPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "RecursiveArrayTools"]
diff --git a/src/Manifolds.jl b/src/Manifolds.jl
index 73ed5d7d93..e43e215606 100644
--- a/src/Manifolds.jl
+++ b/src/Manifolds.jl
@@ -345,7 +345,6 @@ using ManifoldDiff:
     _jacobian!,
     set_default_differential_backend!
 using ManifoldDiff:
-    AbstractDiffBackend,
     AbstractRiemannianDiffBackend,
     CoprojectorOntoVector,
     ExplicitEmbeddedBackend,
@@ -365,6 +364,8 @@ using Statistics
 using StatsBase
 using StatsBase: AbstractWeights
 
+const AbstractDiffBackend = Any  # TODO: remove
+
 include("utils.jl")
 
 include("manifold_fallbacks.jl")
diff --git a/test/differentiation.jl b/test/differentiation.jl
index 3db2a3745c..ffa5b108b3 100644
--- a/test/differentiation.jl
+++ b/test/differentiation.jl
@@ -32,18 +32,19 @@ function ManifoldDiff.gradient(::AbstractManifold, f, p, ::TestRiemannianBackend
     return collect(1.0:length(p))
 end
 
+using ADTypes
 using FiniteDifferences
 using LinearAlgebra: Diagonal, dot
 
 @testset "Differentiation backend" begin
-    fd51 = ManifoldDiff.FiniteDifferencesBackend()
+    fd51 = AutoFiniteDifferences(central_fdm(5, 1))
     @testset "default_differential_backend" begin
-        @test default_differential_backend() isa ManifoldDiff.FiniteDifferencesBackend
+        @test default_differential_backend() isa AutoFiniteDifferences
 
-        @test length(fd51.method.grid) == 5
+        @test length(fd51.fdm.grid) == 5
         # check method order
-        @test typeof(fd51.method).parameters[2] == 1
-        fd71 = ManifoldDiff.FiniteDifferencesBackend(central_fdm(7, 1))
+        @test typeof(fd51.fdm).parameters[2] == 1
+        fd71 = AutoFiniteDifferences(central_fdm(7, 1))
         @test set_default_differential_backend!(fd71) == fd71
         @test default_differential_backend() == fd71
     end
@@ -105,7 +106,7 @@ rb_onb_default = TangentDiffBackend(
     DefaultOrthonormalBasis(),
 )
 
-rb_onb_fd51 = TangentDiffBackend(ManifoldDiff.FiniteDifferencesBackend())
+rb_onb_fd51 = TangentDiffBackend(AutoFiniteDifferences(central_fdm(5, 1)))
 
 rb_onb_default2 = TangentDiffBackend(
     default_differential_backend();
diff --git a/test/manifolds/symplectic.jl b/test/manifolds/symplectic.jl
index deaa1c8767..7fe4accc65 100644
--- a/test/manifolds/symplectic.jl
+++ b/test/manifolds/symplectic.jl
@@ -192,7 +192,7 @@ using ManifoldDiff
         analytical_grad_f(p) = (1 / 2) * (p * J * p * J + p * p')
 
         p_grad = points[1]
-        fd_diff = RiemannianProjectionBackend(ManifoldDiff.FiniteDifferencesBackend())
+        fd_diff = RiemannianProjectionBackend(AutoFiniteDifferences(central_fdm(5, 1)))
 
         @test isapprox(
             Manifolds.gradient(Sp_6, test_f, p_grad, fd_diff),
diff --git a/test/manifolds/symplecticstiefel.jl b/test/manifolds/symplecticstiefel.jl
index 23250d153e..c23cc6c9f0 100644
--- a/test/manifolds/symplecticstiefel.jl
+++ b/test/manifolds/symplecticstiefel.jl
@@ -291,7 +291,7 @@ end
             return Q_grad * p * (euc_grad_f') * Q_grad * p + euc_grad_f * p' * p
         end
         p_grad = convert(Array{Float64}, points[1])
-        fd_diff = RiemannianProjectionBackend(ManifoldDiff.FiniteDifferencesBackend())
+        fd_diff = RiemannianProjectionBackend(AutoFiniteDifferences(central_fdm(5, 1)))
 
         @test isapprox(
             Manifolds.gradient(M, test_f, p_grad, fd_diff),
diff --git a/test/metric.jl b/test/metric.jl
index 8422dbda9b..f1b9f59c18 100644
--- a/test/metric.jl
+++ b/test/metric.jl
@@ -5,7 +5,7 @@ using ManifoldsBase: TraitList
 import ManifoldsBase: default_retraction_method
 import Manifolds: solve_exp_ode
 using Manifolds: InducedBasis, connection, get_chart_index, induced_basis, mean!, median!
-using ManifoldDiff: FiniteDifferencesBackend
+using ADTypes
 include("header.jl")
 
 struct TestEuclidean{N} <: AbstractManifold{ℝ} end
@@ -372,7 +372,7 @@ Manifolds.inner(::MetricManifold{ℝ,<:AbstractManifold{ℝ},Issue539Metric}, p,
             @test gaussian_curvature(M, p, B_chart_p) ≈ 0 atol = 1e-6
             @test einstein_tensor(M, p, B_chart_p) ≈ zeros(n, n) atol = 1e-6
 
-            fdm = FiniteDifferencesBackend(forward_fdm(2, 1))
+            fdm = AutoFiniteDifferences(forward_fdm(2, 1))
             @test christoffel_symbols_first(M, p, B_chart_p; backend=fdm) ≈ zeros(n, n, n) atol =
                 1e-6
             @test christoffel_symbols_second(M, p, B_chart_p; backend=fdm) ≈ zeros(n, n, n) atol =
@@ -384,7 +384,7 @@ Manifolds.inner(::MetricManifold{ℝ,<:AbstractManifold{ℝ},Issue539Metric}, p,
             @test gaussian_curvature(M, p, B_chart_p; backend=fdm) ≈ 0 atol = 1e-6
             @test einstein_tensor(M, p, B_chart_p; backend=fdm) ≈ zeros(n, n) atol = 1e-6
 
-            fd_diff = FiniteDifferencesBackend()
+            fd_diff = AutoFiniteDifferences(central_fdm(5, 1))
             @test christoffel_symbols_first(M, p, B_chart_p; backend=fd_diff) ≈
                   zeros(n, n, n) atol = 1e-6
             @test christoffel_symbols_second(M, p, B_chart_p; backend=fd_diff) ≈