diff --git a/RELEASES.md b/RELEASES.md
index 30c248044..387527987 100644
--- a/RELEASES.md
+++ b/RELEASES.md
@@ -9,6 +9,7 @@
 + The `linspace` method of the backends now has the `type_as` argument to convert to the same dtype and device. (PR #533)
 + The `convolutional_barycenter2d` and `convolutional_barycenter2d_debiased` functions now work with different devices.. (PR #533)
 + New API for Gromov-Wasserstein solvers with `ot.solve_gromov` function (PR #536)
++ New LP solvers from scipy used by default for LP barycenter (PR #537)
 
 #### Closed issues
 - Fix line search evaluating cost outside of the interpolation range (Issue #502, PR #504)
diff --git a/examples/barycenters/plot_barycenter_lp_vs_entropic.py b/examples/barycenters/plot_barycenter_lp_vs_entropic.py
index 6502f16d8..9c4ca4551 100644
--- a/examples/barycenters/plot_barycenter_lp_vs_entropic.py
+++ b/examples/barycenters/plot_barycenter_lp_vs_entropic.py
@@ -83,7 +83,7 @@
 
 
 ot.tic()
-bary_wass2 = ot.lp.barycenter(A, M, weights, solver='interior-point', verbose=True)
+bary_wass2 = ot.lp.barycenter(A, M, weights)
 ot.toc()
 
 pl.figure(2)
@@ -149,7 +149,7 @@
 
 
 ot.tic()
-bary_wass2 = ot.lp.barycenter(A, M, weights, solver='interior-point', verbose=True)
+bary_wass2 = ot.lp.barycenter(A, M, weights)
 ot.toc()
 
 
@@ -223,7 +223,7 @@
 
 
 ot.tic()
-bary_wass2 = ot.lp.barycenter(A, M, weights, solver='interior-point', verbose=True)
+bary_wass2 = ot.lp.barycenter(A, M, weights)
 ot.toc()
 
 
diff --git a/ot/lp/cvx.py b/ot/lp/cvx.py
index 3f7eb36c6..f9572962a 100644
--- a/ot/lp/cvx.py
+++ b/ot/lp/cvx.py
@@ -25,7 +25,7 @@ def scipy_sparse_to_spmatrix(A):
     return SP
 
 
-def barycenter(A, M, weights=None, verbose=False, log=False, solver='interior-point'):
+def barycenter(A, M, weights=None, verbose=False, log=False, solver='highs-ipm'):
     r"""Compute the Wasserstein barycenter of distributions A
 
      The function solves the following optimization problem [16]:
@@ -115,13 +115,13 @@ def barycenter(A, M, weights=None, verbose=False, log=False, solver='interior-po
     A_eq = sps.vstack((A_eq1, A_eq2))
     b_eq = np.concatenate((b_eq1, b_eq2))
 
-    if not cvxopt or solver in ['interior-point']:
+    if not cvxopt or solver in ['interior-point', 'highs', 'highs-ipm', 'highs-ds']:
         # cvxopt not installed or interior point
 
         if solver is None:
             solver = 'interior-point'
 
-        options = {'sparse': True, 'disp': verbose}
+        options = {'disp': verbose}
         sol = sp.optimize.linprog(c, A_eq=A_eq, b_eq=b_eq, method=solver,
                                   options=options)
         x = sol.x
diff --git a/requirements.txt b/requirements.txt
index fd39cbab4..6ac25eb3c 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
 numpy>=1.20
-scipy>=1.3
+scipy>=1.6
 matplotlib
 autograd
 pymanopt
diff --git a/setup.py b/setup.py
index dc9066d62..bf59c7d6d 100644
--- a/setup.py
+++ b/setup.py
@@ -69,7 +69,7 @@
     scripts=[],
     data_files=[],
     setup_requires=["oldest-supported-numpy", "cython>=0.23"],
-    install_requires=["numpy>=1.16", "scipy>=1.0"],
+    install_requires=["numpy>=1.16", "scipy>=1.6"],
     python_requires=">=3.6",
     classifiers=[
         'Development Status :: 5 - Production/Stable',