diff --git a/Project.toml b/Project.toml
index 2930dc3b..6875c9bb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ForwardDiff"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.30"
+version = "0.10.31"
 
 [deps]
 CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950"
@@ -24,7 +24,7 @@ DiffTests = "0.0.1, 0.1"
 LogExpFunctions = "0.3"
 NaNMath = "0.2.2, 0.3, 1"
 Preferences = "1"
-SpecialFunctions = "0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.10, 1.0, 2"
+SpecialFunctions = "0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.8.3, 0.9, 0.10, 0.11, 0.12, 1.0"
 julia = "1"
 
diff --git a/src/dual.jl b/src/dual.jl
index 8e99c9fe..42682b42 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -755,9 +755,9 @@ function LinearAlgebra.eigen(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:Rea
     Eigen(λ,Dual{Tg}.(Q, tuple.(parts...)))
 end
 
-# SpecialFunctions.logabsgamma           #
-# Derivative is not defined in DiffRules #
-#----------------------------------------#
+# Functions in SpecialFunctions which return tuples #
+# Their derivatives are not defined in DiffRules    #
+#---------------------------------------------------#
 
 function SpecialFunctions.logabsgamma(d::Dual{T,<:Real}) where {T}
     x = value(d)
@@ -765,6 +765,14 @@ function SpecialFunctions.logabsgamma(d::Dual{T,<:Real}) where {T}
     return (Dual{T}(y, SpecialFunctions.digamma(x) * partials(d)), s)
 end
 
+# Derivatives wrt to first parameter and precision setting are not supported
+function SpecialFunctions.gamma_inc(a::Real, d::Dual{T,<:Real}, ind::Integer) where {T}
+    x = value(d)
+    p, q = SpecialFunctions.gamma_inc(a, x, ind)
+    ∂p = exp(-x) * x^(a - 1) / SpecialFunctions.gamma(a) * partials(d)
+    return (Dual{T}(p, ∂p), Dual{T}(q, -∂p))
+end
+
 ###################
 # Pretty Printing #
 ###################
diff --git a/test/DualTest.jl b/test/DualTest.jl
index 66a1866a..eca0da2f 100644
--- a/test/DualTest.jl
+++ b/test/DualTest.jl
@@ -540,8 +540,29 @@ for N in (0,3), M in (0,4), V in (Int, Float32)
         @test dual_isapprox(f(PRIMAL, PRIMAL2, FDNUM3), Dual{TestTag()}(f(PRIMAL, PRIMAL2, PRIMAL3), PARTIALS3))
     end
 
+    # Functions in Specialfunctions that return tuples and
+    # therefore are not supported by DiffRules
     @test dual_isapprox(logabsgamma(FDNUM)[1], loggamma(abs(FDNUM)))
     @test dual_isapprox(logabsgamma(FDNUM)[2], sign(gamma(FDNUM)))
+
+    a = rand(float(V))
+    fdnum = Dual{TestTag()}(1 + PRIMAL, PARTIALS) # 1 + PRIMAL avoids issues with finite differencing close to 0
+    for ind in ((), (0,), (1,), (2,))
+        # Only test if primal method exists
+        # (e.g., older versions of SpecialFunctions don't define `gamma_inc(a, x)` but only `gamma_inc(a, x, ind)`
+        hasmethod(gamma_inc, typeof((a, 1 + PRIMAL, ind...))) || continue
+
+        pq = gamma_inc(a, fdnum, ind...)
+        @test pq isa Tuple{Dual{TestTag()},Dual{TestTag()}}
+        # We have to adjust tolerances if lower accuracy is requested
+        # Therefore we don't use `dual_isapprox`
+        tol = V === Float32 ? 5f-4 : 1e-6
+        tol = tol^(one(tol) / 2^(isempty(ind) ? 0 : first(ind)))
+        for i in 1:2
+            @test value(pq[i]) ≈ gamma_inc(a, 1 + PRIMAL, ind...)[i] rtol=tol
+            @test partials(pq[i]) ≈ PARTIALS * Calculus.derivative(x -> gamma_inc(a, x, ind...)[i], 1 + PRIMAL) rtol=tol
+        end
+    end
 end
 
 @testset "Exponentiation of zero" begin