diff --git a/DESCRIPTION b/DESCRIPTION
index d145534..a74de2a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -9,7 +9,7 @@ Authors@R: c(person("Lam Si Tung", "Ho", role=c("aut", "cre"), email="lamho86@gm
       person("Rachel", "Feldman",role="ctb"),
       person("Qing", "Yu",role="ctb"))
 Depends: R (>= 3.0), ape
-Imports: parallel
+Imports: future.apply
 Description: Provides functions for fitting phylogenetic linear models and phylogenetic generalized linear models. The computation uses an algorithm that is linear in the number of tips in the tree. The package also provides functions for simulating continuous or binary traits along the tree. Other tools include functions to test the adequacy of a population tree.
 License: GPL (>= 2) | file LICENSE
 URL: https://CRAN.R-project.org/package=phylolm
diff --git a/R/phyloglm.R b/R/phyloglm.R
index 6721774..a70e3ea 100644
--- a/R/phyloglm.R
+++ b/R/phyloglm.R
@@ -1,6 +1,6 @@
 phyloglm <- function(formula, data=list(), phy, method=c("logistic_MPLE","logistic_IG10","poisson_GEE"),
                      btol = 10, log.alpha.bound = 4, start.beta=NULL, start.alpha=NULL, 
-                     boot = 0, full.matrix = TRUE, parallel = NULL)
+                     boot = 0, full.matrix = TRUE)
 {
   ### initialize	
   if (!inherits(phy, "phylo")) stop("object \"phy\" is not of class \"phylo\".")	
@@ -503,21 +503,7 @@ You can increase this bound by increasing 'btol'.")
       return(bootvector)
     }
     
-    ## set up cluster for parallel programming if needed
-    if (!is.null(parallel)) {
-      # create a cluster of type `parallel`, avoid using too many nodes
-      cl <- parallel::makeCluster(min(c(parallel::detectCores() - 1, boot)), parallel)
-      # load all data and variables on the nodes
-      parallel::clusterExport(cl, ls(), environment())
-      # guarantee safe exit
-      on.exit(parallel::stopCluster(cl))
-    }
-    
-    if (is.null(parallel)) {
-      bootmatrix <- lapply(as.data.frame(bootobject), boot_model)
-    } else {
-      bootmatrix <- parallel::parLapply(cl, as.data.frame(bootobject), boot_model)
-    }
+    bootmatrix <- future.apply::future_lapply(as.data.frame(bootobject), boot_model)
     bootmatrix <- do.call(rbind, bootmatrix)
     
     # summarize bootstrap estimates
diff --git a/R/phylolm.R b/R/phylolm.R
index 38c0bfb..8aee883 100644
--- a/R/phylolm.R
+++ b/R/phylolm.R
@@ -1,7 +1,7 @@
 phylolm <- function(formula, data=list(), phy, 
 	model=c("BM","OUrandomRoot","OUfixedRoot","lambda","kappa","delta","EB","trend"),
 	lower.bound=NULL, upper.bound=NULL, starting.value=NULL, measurement_error = FALSE,
-	boot=0,full.matrix = TRUE, parallel = NULL, ...)
+	boot=0,full.matrix = TRUE, ...)
 {
 
   ## initialize	
@@ -385,21 +385,7 @@ phylolm <- function(formula, data=list(), phy,
       return(bootvector)
     }
     
-    ## set up cluster for parallel programming if needed
-    if (!is.null(parallel)) {
-      # create a cluster of type `parallel`, avoid using too many nodes
-      cl <- parallel::makeCluster(min(c(parallel::detectCores() - 1, boot)), parallel)
-      # load all data and variables on the nodes
-      parallel::clusterExport(cl, ls(), environment())
-      # guarantee safe exit
-      on.exit(parallel::stopCluster(cl))
-    }
-    
-    if (is.null(parallel)) {
-      bootmatrix <- lapply(as.data.frame(booty), boot_model)
-    } else {
-      bootmatrix <- parallel::parLapply(cl, as.data.frame(booty), boot_model)
-    }
+    bootmatrix <- future.apply::future_lapply(as.data.frame(booty), boot_model)
     bootmatrix <- do.call(rbind, bootmatrix)
     
     # summarize bootstrap estimates
diff --git a/man/phyloglm.Rd b/man/phyloglm.Rd
index 5c6d5aa..4f163f6 100644
--- a/man/phyloglm.Rd
+++ b/man/phyloglm.Rd
@@ -10,7 +10,7 @@
 phyloglm(formula, data, phy, method = c("logistic_MPLE","logistic_IG10",
          "poisson_GEE"), btol = 10, log.alpha.bound = 4,
          start.beta=NULL, start.alpha=NULL,
-         boot = 0, full.matrix = TRUE, parallel = NULL)
+         boot = 0, full.matrix = TRUE)
 }
 
 \arguments{
@@ -25,10 +25,14 @@ phyloglm(formula, data, phy, method = c("logistic_MPLE","logistic_IG10",
   \item{start.alpha}{(logistic regression only) starting values for alpha (phylogenetic correlation).}
   \item{boot}{number of independent bootstrap replicates, \code{0} means no bootstrap.}
   \item{full.matrix}{if \code{TRUE}, the full matrix of bootstrap estimates (coefficients and alpha) will be returned.}
-  \item{parallel}{bootstrapping can be performed on a cluster by specifying the type of cluster (e.g. \code{"SOCK"})}
 }
 \details{
-This function uses an algorithm that is linear in the number of tips in the tree.
+  This function uses an algorithm that is linear in the number of tips in the tree.
+
+  Bootstrapping can be parallelized using the \code{future} package on any future 
+  compatible back-end. For example, run \code{library(future); plan(multiprocess))}, 
+  after which bootstrapping will automatically occur in parallel. See 
+  \code{\link[future]{plan}} for options.
 }
 \value{
   \item{coefficients}{the named vector of coefficients.}
diff --git a/man/phylolm.Rd b/man/phylolm.Rd
index 27e1bc5..4c64c95 100644
--- a/man/phylolm.Rd
+++ b/man/phylolm.Rd
@@ -9,7 +9,7 @@ phylolm(formula, data = list(), phy, model = c("BM", "OUrandomRoot",
        "OUfixedRoot", "lambda", "kappa", "delta", "EB", "trend"),
        lower.bound = NULL, upper.bound = NULL,
        starting.value = NULL, measurement_error = FALSE,
-       boot=0,full.matrix = TRUE,parallel=NULL, ...)
+       boot=0,full.matrix = TRUE, ...)
 }
 
 \arguments{
@@ -24,7 +24,6 @@ phylolm(formula, data = list(), phy, model = c("BM", "OUrandomRoot",
   \item{measurement_error}{a logical value indicating whether there is measurement error \code{sigma2_error} (see Details).}
   \item{boot}{number of independent bootstrap replicates, 0 means no bootstrap.}
   \item{full.matrix}{if \code{TRUE}, the full matrix of bootstrap estimates (coefficients and covariance parameters) will be returned.}
-  \item{parallel}{bootstrapping can be performed on a cluster by specifying the type of cluster (e.g. \code{"SOCK"})}
   \item{\dots}{further arguments to be passed to the function \code{optim}.}
 }
 \details{This function uses an algorithm that is linear in the number of
@@ -53,6 +52,11 @@ phylolm(formula, data = list(), phy, model = c("BM", "OUrandomRoot",
   \eqn{[10^{-5},3]} for \eqn{\delta}{delta} and
   \eqn{[-3/T,0]} for \code{rate}, where \eqn{T} is the mean root-to-tip distance.
   \eqn{[10^{-16}, 10^{16}]} for the ratio \code{sigma2_error}/\code{sigma2} (if measurement errors is used).
+  
+  Bootstrapping can be parallelized using the \code{future} package on any future 
+  compatible back-end. For example, run \code{library(future); plan(multiprocess))}, 
+  after which bootstrapping will automatically occur in parallel. See 
+  \code{\link[future]{plan}} for options.
 }
 \value{
    \item{coefficients}{the named vector of coefficients.}