diff --git a/Project.toml b/Project.toml
index 10b16db8d..0b4afa33d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -39,10 +39,12 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [weakdeps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 
 [extensions]
+SymbolicsGroebnerExt = "Groebner"
 SymbolicsPreallocationToolsExt = ["ForwardDiff", "PreallocationTools"]
 SymbolicsSymPyExt = "SymPy"
 
@@ -56,7 +58,7 @@ Distributions = "0.25"
 DocStringExtensions = "0.9"
 DomainSets = "0.6"
 DynamicPolynomials = "0.5"
-Groebner = "0.5"
+Groebner = "0.5, 0.6"
 IfElse = "0.1"
 LaTeXStrings = "1.3"
 LambertW = "0.4.5"
@@ -81,6 +83,7 @@ julia = "1.10"
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
@@ -91,4 +94,4 @@ SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "SafeTestsets", "Pkg", "PkgBenchmark", "PreallocationTools", "ForwardDiff", "BenchmarkTools", "ReferenceTests", "SymPy", "Random"]
+test = ["Test", "SafeTestsets", "Pkg", "PkgBenchmark", "PreallocationTools", "ForwardDiff", "Groebner", "BenchmarkTools", "ReferenceTests", "SymPy", "Random"]
diff --git a/ext/SymbolicsGroebnerExt.jl b/ext/SymbolicsGroebnerExt.jl
new file mode 100644
index 000000000..2ebb89b4a
--- /dev/null
+++ b/ext/SymbolicsGroebnerExt.jl
@@ -0,0 +1,73 @@
+module SymbolicsGroebnerExt
+
+using Groebner
+
+if isdefined(Base, :get_extension)
+    using Symbolics
+    using Symbolics: Num, symtype
+else
+    using ..Symbolics
+    using ..Symbolics: Num, symtype
+end
+
+"""
+    groebner_basis(polynomials; kwargs...)
+
+Computes a Groebner basis of the ideal generated by the given `polynomials`
+using Groebner.jl as the backend.
+
+The basis is guaranteed to be unique.
+The algorithm is randomized, and the output is correct with high probability.
+
+If a coefficient in the resulting basis becomes too large to be represented
+exactly, `DomainError` is thrown.
+
+## Optional Arguments
+
+The Groebner.jl backend provides a number of useful keyword arguments, which are
+also available for this function. See `?Groebner.groebner`.
+
+## Example
+
+```jldoctest
+julia> using Symbolics, Groebner
+
+julia> @variables x y;
+
+julia> groebner_basis([x*y^2 + x, x^2*y + y])
+```
+"""
+function Symbolics.groebner_basis(polynomials::Vector{Num}; kwargs...)
+    polynoms, pvar2sym, sym2term = Symbolics.symbol_to_poly(polynomials)
+    basis = Groebner.groebner(polynoms; kwargs...)
+    PolyType = symtype(first(polynomials))
+    Symbolics.poly_to_symbol(basis, pvar2sym, sym2term, PolyType)
+end
+
+"""
+    is_groebner_basis(polynomials; kwargs...)
+
+Checks whether the given `polynomials` forms a Groebner basis using Groebner.jl
+as the backend.
+
+## Optional Arguments
+
+The Groebner.jl backend provides a number of useful keyword arguments, which are
+also available for this function. See `?Groebner.isgroebner`.
+
+## Example
+
+```jldoctest
+julia> using Symbolics, Groebner
+
+julia> @variables x y;
+
+julia> is_groebner_basis([x^2 - y^2, x*y^2 + x, y^3 + y])
+```
+"""
+function Symbolics.is_groebner_basis(polynomials::Vector{Num}; kwargs...)
+    polynoms, _, _ = Symbolics.symbol_to_poly(polynomials)
+    Groebner.isgroebner(polynoms; kwargs...)
+end
+
+end # module
diff --git a/src/Symbolics.jl b/src/Symbolics.jl
index 82062a688..d453afd73 100644
--- a/src/Symbolics.jl
+++ b/src/Symbolics.jl
@@ -78,6 +78,8 @@ include("equations.jl")
 export Inequality, ≲, ≳
 include("inequality.jl")
 
+using Bijections
+import DynamicPolynomials
 include("utils.jl")
 
 using ConstructionBase
@@ -116,9 +118,8 @@ include("logexpfunctions-lib.jl")
 
 include("linear_algebra.jl")
 
-using Groebner
 include("groebner_basis.jl")
-export groebner_basis
+export groebner_basis, is_groebner_basis
 
 import Libdl
 include("build_function.jl")
@@ -182,6 +183,9 @@ end
 
 @static if !isdefined(Base,:get_extension)
     function __init__()
+        @require Groebner="0b43b601-686d-58a3-8a1c-6623616c7cd4" begin
+            include("../ext/SymbolicsGroebnerExt.jl")
+        end
         @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" begin
             include("../ext/SymbolicsSymPyExt.jl")
         end
diff --git a/src/groebner_basis.jl b/src/groebner_basis.jl
index 6446ab4f7..5bb9c5e5f 100644
--- a/src/groebner_basis.jl
+++ b/src/groebner_basis.jl
@@ -1,104 +1,25 @@
-import DynamicPolynomials
-using Bijections
 
-const DP = DynamicPolynomials
-# extracting underlying polynomial and coefficient type from Polyforms
-underlyingpoly(x::Number) = x
-underlyingpoly(pf::PolyForm) = pf.p
-coefftype(x::Number) = typeof(x)
-coefftype(pf::PolyForm) = DP.coefficienttype(underlyingpoly(pf))
-
-#=
-Converts an array of symbolic polynomials
-into an array of DynamicPolynomials.Polynomials
-=#
-function symbol_to_poly(sympolys::AbstractArray)
-    @assert !isempty(sympolys) "Empty input."
-
-    # standardize input
-    stdsympolys = map(unwrap, sympolys)
-    sort!(stdsympolys, lt=(<ₑ))
-
-    pvar2sym = Bijections.Bijection{Any,Any}()
-    sym2term = Dict{BasicSymbolic,Any}()
-    polyforms = map(f -> PolyForm(f, pvar2sym, sym2term), stdsympolys)
-
-    # Discover common coefficient type
-    commontype = mapreduce(coefftype, promote_type, polyforms, init=Int)
-    @assert commontype <: Union{Integer,Rational} "Only integer and rational coefficients are supported as input."
-
-    # Convert all to DP.Polynomial, so that coefficients are of same type,
-    # and constants are treated as polynomials
-    # We also need this because Groebner does not support abstract types as input
-    polynoms = Vector{DP.Polynomial{DP.Commutative{DP.CreationOrder}, DP.Graded{DP.LexOrder}, commontype}}(undef, length(sympolys))
-    for (i, pf) in enumerate(polyforms)
-        polynoms[i] = underlyingpoly(pf)
-    end
-
-    polynoms, pvar2sym, sym2term
-end
-
-#=
-Converts an array of AbstractPolynomialLike`s into an array of
-symbolic expressions mapping variables w.r.t pvar2sym
-=#
-function poly_to_symbol(polys, pvar2sym, sym2term, ::Type{T}) where {T}
-    map(f -> PolyForm{T}(f, pvar2sym, sym2term), polys)
-end
+@noinline __throw_absent_groebner_engine() = throw(
+    """Groebner bases engine is required. Execute `using Groebner` to enable this functionality.""")
 
 """
-
-groebner_basis(polynomials)
+    groebner_basis(polynomials)
 
 Computes a Groebner basis of the ideal generated by the given `polynomials`.
-The basis is reduced, thus guaranteed to be unique.
-
-# Example
-
-```jldoctest
-julia> using Symbolics
-
-julia> @variables x y;
-
-julia> groebner_basis([x*y^2 + x, x^2*y + y])
-```
-
-The coefficients in the resulting basis are in the same domain as for input polynomials.
-Hence, if the coefficient becomes too large to be represented exactly, `DomainError` is thrown.
-
-The algorithm is randomized, so the basis will be correct with high probability.
 
+This function requires a Groebner bases backend (such as Groebner.jl) to be loaded.
 """
-function groebner_basis(polynomials)
-    polynoms, pvar2sym, sym2term = symbol_to_poly(polynomials)
-
-    basis = groebner(polynoms)
-
-    # polynomials is nonemtpy
-    T = symtype(first(polynomials))
-    poly_to_symbol(basis, pvar2sym, sym2term, T)
+function groebner_basis(args; kwargs...)
+    __throw_absent_groebner_engine()
 end
 
 """
-
-Do not document this for now.
-
-is_groebner_basis(polynomials)
+    groebner_basis(polynomials)
 
 Checks whether the given `polynomials` forms a Groebner basis.
 
-# Example
-
-```jldoctest
-julia> using Symbolics
-
-julia> @variables x y;
-
-julia> is_groebner_basis([x^2 - y^2, x*y^2 + x, y^3 + y])
-```
-
+This function requires a Groebner bases backend (such as Groebner.jl) to be loaded.
 """
-function is_groebner_basis(polynomials)
-    polynoms, _, _ = symbol_to_poly(polynomials)
-    isgroebner(polynoms)
+function is_groebner_basis(args; kwargs...)
+    __throw_absent_groebner_engine()
 end
diff --git a/src/utils.jl b/src/utils.jl
index 94b9433a7..f86ed12c3 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -331,3 +331,50 @@ function coeff(p, sym=nothing)
         throw(DomainError(p, "Datatype $(typeof(p)) not accepted."))
     end
 end
+
+### Nums <--> Polys
+
+const DP = DynamicPolynomials
+# extracting underlying polynomial and coefficient type from Polyforms
+underlyingpoly(x::Number) = x
+underlyingpoly(pf::PolyForm) = pf.p
+coefftype(x::Number) = typeof(x)
+coefftype(pf::PolyForm) = DP.coefficient_type(underlyingpoly(pf))
+
+#=
+Converts an array of symbolic polynomials
+into an array of DynamicPolynomials.Polynomials
+=#
+function symbol_to_poly(sympolys::AbstractArray)
+    @assert !isempty(sympolys) "Empty input."
+
+    # standardize input
+    stdsympolys = map(unwrap, sympolys)
+    sort!(stdsympolys, lt=(<ₑ))
+
+    pvar2sym = Bijections.Bijection{Any,Any}()
+    sym2term = Dict{BasicSymbolic,Any}()
+    polyforms = map(f -> PolyForm(f, pvar2sym, sym2term), stdsympolys)
+
+    # Discover common coefficient type
+    commontype = mapreduce(coefftype, promote_type, polyforms, init=Int)
+    @assert commontype <: Union{Integer,Rational} "Only integer and rational coefficients are supported as input."
+
+    # Convert all to DP.Polynomial, so that coefficients are of same type,
+    # and constants are treated as polynomials
+    # We also need this because Groebner does not support abstract types as input
+    polynoms = Vector{DP.Polynomial{DP.Commutative{DP.CreationOrder},DP.Graded{DP.LexOrder},commontype}}(undef, length(sympolys))
+    for (i, pf) in enumerate(polyforms)
+        polynoms[i] = underlyingpoly(pf)
+    end
+
+    polynoms, pvar2sym, sym2term
+end
+
+#=
+Converts an array of AbstractPolynomialLike`s into an array of
+symbolic expressions mapping variables w.r.t pvar2sym
+=#
+function poly_to_symbol(polys, pvar2sym, sym2term, ::Type{T}) where {T}
+    map(f -> PolyForm{T}(f, pvar2sym, sym2term), polys)
+end
diff --git a/test/extensions/groebner.jl b/test/extensions/groebner.jl
new file mode 100644
index 000000000..93b9d02f9
--- /dev/null
+++ b/test/extensions/groebner.jl
@@ -0,0 +1,76 @@
+using Symbolics
+using Symbolics: symbol_to_poly, poly_to_symbol
+using Groebner
+using Test
+
+@variables x y z
+
+syms = [
+    [1], [1, BigInt(2)], [x], [x^2 + 2], [x + (5 // 8)y],
+    [1 - x * y * z], [(x + y + z)^5], [0, BigInt(0)y^2, Rational(1)z^3],
+    [x, sin((44 // 31)y) * z]
+]
+for sym in syms
+    polynoms, pvar2sym, sym2term = Symbolics.symbol_to_poly(sym)
+    sym2 = Symbolics.poly_to_symbol(polynoms, pvar2sym, sym2term, Real)
+    @test isequal(expand.(sym2), expand.(sym))
+end
+
+@test isequal(expand.(groebner_basis([x, y])), [y, x])
+@test isequal(expand.(groebner_basis([y, x])), [y, x])
+@test isequal(expand.(groebner_basis([x])), [x])
+@test isequal(expand.(groebner_basis([x, x^2])), [x])
+@test isequal(expand.(groebner_basis([BigInt(1)x + 2 // 3])), [x + 2 // 3])
+
+@test Symbolics.is_groebner_basis([x, y, z])
+@test Symbolics.is_groebner_basis([x^2 - x, y^2 - y])
+@test !Symbolics.is_groebner_basis([x^2 + y, x * y^3 - 1, y^4 - 1])
+
+@variables x1 x2 x3 x4
+@test isequal(expand.(groebner_basis([x1, x, y])), [y, x1, x])
+
+# input unchanged
+f1 = [x2, x1, x4, x3]
+f2 = [x2, x1, x4, x3]
+groebner_basis(f1)
+@test isequal(expand.(f1), expand.(f2))
+
+@variables x1 x2 x3 x4 x5
+system = [
+    x1 + x2 + x3 + x4 + x5,
+    x1 * x2 + x1 * x3 + x1 * x4 + x1 * x5 + x2 * x3 + x2 * x4 + x2 * x5 + x3 * x4 + x3 * x5 + x4 * x5,
+    x1 * x2 * x3 + x1 * x2 * x4 + x1 * x2 * x5 + x1 * x3 * x4 + x1 * x3 * x5 + x1 * x4 * x5 + x2 * x3 * x4 + x2 * x3 * x5 + x2 * x4 * x5 + x3 * x4 * x5,
+    x1 * x2 * x3 * x4 + x1 * x2 * x3 * x5 + x1 * x2 * x4 * x5 + x1 * x3 * x4 * x5 + x2 * x3 * x4 * x5,
+    x1 * x2 * x3 * x4 * x5 - 1
+]
+truebasis = [
+    x1 + x2 + x3 + x4 + x5,
+    x2^2 + x2 * x3 + x2 * x4 + x2 * x5 + x3^2 + x3 * x4 + x3 * x5 + x4^2 + x4 * x5 + x5^2,
+    x3^3 + x3 * (x4^2) + x3 * (x5^2) + x4^3 + x4 * (x3^2) + x5 * (x4^2) + x4 * (x5^2) + x5^3 + x5 * (x3^2) + x3 * x4 * x5,
+    x4^4 + x4 * (x5^3) + x5^4 + x5 * (x4^3) + (x4^2) * (x5^2),
+    x5^5 - 1
+]
+basis = expand.(groebner_basis(system))
+@test isequal(basis, truebasis)
+
+basis = expand.(groebner_basis(system, linalg=:deterministic))
+@test isequal(basis, truebasis)
+
+N = 45671930739135174346839766056203605080877915151
+system = [
+    x1 + x2 + x3 + x4,
+    x1 * x2 + x1 * x3 + x1 * x4 + x2 * x3 + x2 * x4 + x3 * x4,
+    x1 * x2 * x3 + x1 * x2 * x4 + x1 * x3 * x4 + x2 * x3 * x4,
+    x1 * x2 * x3 * x4 + N
+]
+truebasis = [
+    x1 + x2 + x3 + x4,
+    x2^2 + x2 * x3 + x3^2 + x2 * x4 + x3 * x4 + x4^2,
+    x3^3 + x3^2 * x4 + x3 * x4^2 + x4^3,
+    x4^4 - N
+]
+basis = groebner_basis(system)
+@test isequal(expand.(basis), truebasis)
+
+# Groebner does not yet work with constant ideals
+@test_broken groebner_basis([1])
diff --git a/test/groebner_basis.jl b/test/groebner_basis.jl
deleted file mode 100644
index 638f5c3d1..000000000
--- a/test/groebner_basis.jl
+++ /dev/null
@@ -1,75 +0,0 @@
-using Symbolics
-using Symbolics: symbol_to_poly, poly_to_symbol
-
-using Test
-
-@variables x y z
-
-syms = [
-    [1], [1, BigInt(2)], [x], [x^2 + 2], [x + (5//8)y],
-    [1 - x*y*z], [(x + y + z)^5], [0, BigInt(0)y^2, Rational(1)z^3],
-    [x, sin((44//31)y)*z]
-]
-for sym in syms
-    polynoms, pvar2sym, sym2term = Symbolics.symbol_to_poly(sym)
-    sym2 = Symbolics.poly_to_symbol(polynoms, pvar2sym, sym2term, Real)
-    @test isequal(expand.(sym2), expand.(sym))
-end
-
-@test isequal(expand.(groebner_basis([x, y])), [y, x])
-@test isequal(expand.(groebner_basis([y, x])), [y, x])
-@test isequal(expand.(groebner_basis([x])), [x])
-@test isequal(expand.(groebner_basis([x, x^2])), [x])
-@test isequal(expand.(groebner_basis([BigInt(1)x + 2//3])), [x  + 2//3])
-
-@test Symbolics.is_groebner_basis([x, y, z])
-@test Symbolics.is_groebner_basis([x^2 - x, y^2 - y])
-@test !Symbolics.is_groebner_basis([x^2 + y, x*y^3 - 1, y^4 - 1])
-
-@variables x1 x2 x3 x4
-@test isequal(expand.(groebner_basis([x1, x, y])), [y, x1, x])
-
-# input unchanged
-f1 = [x2, x1, x4, x3]
-f2 = [x2, x1, x4, x3]
-groebner_basis(f1)
-@test isequal(expand.(f1), expand.(f2))
-
-@variables x1 x2 x3 x4 x5
-system = [
-   x1 + x2 + x3 + x4 + x5,
-   x1*x2 + x1*x3 + x1*x4 + x1*x5 + x2*x3 + x2*x4 + x2*x5 + x3*x4 + x3*x5 + x4*x5,
-   x1*x2*x3 + x1*x2*x4 + x1*x2*x5 + x1*x3*x4 + x1*x3*x5 + x1*x4*x5 + x2*x3*x4 + x2*x3*x5 + x2*x4*x5 + x3*x4*x5,
-   x1*x2*x3*x4 + x1*x2*x3*x5 + x1*x2*x4*x5 + x1*x3*x4*x5 + x2*x3*x4*x5,
-   x1*x2*x3*x4*x5 - 1
-]
-truebasis = [
-        x1 + x2 + x3 + x4 + x5,
-        x2^2 + x2*x3 + x2*x4 + x2*x5 + x3^2 + x3*x4 + x3*x5 + x4^2 + x4*x5 + x5^2,
-        x3^3 + x3*(x4^2) + x3*(x5^2) + x4^3 + x4*(x3^2) + x5*(x4^2) + x4*(x5^2) + x5^3 + x5*(x3^2) + x3*x4*x5,
-        x4^4 + x4*(x5^3) + x5^4 + x5*(x4^3) + (x4^2)*(x5^2),
-        x5^5 - 1
-]
-basis = expand.(groebner_basis(system))
-@test isequal(basis, truebasis)
-
-
-N = 45671930739135174346839766056203605080877915151
-system = [
-    x1 + x2 + x3 + x4,
-    x1*x2 + x1*x3 + x1*x4 + x2*x3 + x2*x4 + x3*x4,
-    x1*x2*x3 + x1*x2*x4 + x1*x3*x4 + x2*x3*x4,
-    x1*x2*x3*x4 + N
-]
-truebasis = [
-    x1 + x2 + x3 + x4,
-    x2^2 + x2*x3 + x3^2 + x2*x4 + x3*x4 + x4^2,
-    x3^3 + x3^2*x4 + x3*x4^2 + x4^3,
-    x4^4 - N
-]
-basis = groebner_basis(system)
-@test isequal(expand.(basis), truebasis)
-
-
-# Groebner does not yet work with constant ideals
-@test_broken groebner_basis([1])
diff --git a/test/runtests.jl b/test/runtests.jl
index 7ec01e6c4..ecfecc37a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -32,7 +32,6 @@ if GROUP == "All" || GROUP == "Core"
     @safetestset "Is Linear or Affine Test" begin include("islinear_affine.jl") end
     @safetestset "Linear Solver Test" begin include("linear_solver.jl") end
     @safetestset "Algebraic Solver Test" begin include("solver.jl") end
-    @safetestset "Groebner Bases Test" begin include("groebner_basis.jl") end
     @safetestset "Overloading Test" begin include("overloads.jl") end
     @safetestset "Nested ForwardDiff Sparsity Test" begin include("nested_forwarddiff_sparsity.jl") end
     @safetestset "Build Function Test" begin include("build_function.jl") end
@@ -47,6 +46,12 @@ if GROUP == "All" || GROUP == "Core"
     @safetestset "CartesianIndex Test" begin include("cartesianindex.jl") end
     @safetestset "LogExpFunctions Test" begin include("logexpfunctions.jl") end
     @safetestset "Registration without using Test" begin include("registration_without_using.jl") end
+
+    @safetestset "Groebner extension Test" begin include("extensions/groebner.jl") end
+end
+
+if GROUP == "GroebnerExt"
+    @safetestset "Groebner extension Test" begin include("extensions/groebner.jl") end
 end
 
 if GROUP == "All" || GROUP == "Core" || GROUP == "SymbolicIndexingInterface"