diff --git a/R/plotAbundanceDensity.R b/R/plotAbundanceDensity.R
index e5cd5185..8cdf050c 100644
--- a/R/plotAbundanceDensity.R
+++ b/R/plotAbundanceDensity.R
@@ -17,6 +17,10 @@
 #'
 #' @param colour_by Deprecated. Use \code{colour.by} instead.
 #'
+#' @param colour.features \code{Character vector}. A named vector of colors
+#'   for features. Names should match feature names. Useful for matching
+#'   colors across multiple plots. (Default: \code{NULL})
+#'
 #' @param shape.by \code{Character scalar}. Defines a column from
 #'   \code{colData}, that is used to group observations to different point shape
 #'   groups. Must be a value of \code{colData()} function. \code{shape.by} is
@@ -165,6 +169,7 @@ setMethod("plotAbundanceDensity", signature = c(x = "SummarizedExperiment"),
         assay.type = assay_name, assay_name = "counts",
         n = min(nrow(x), 25L), colour.by = colour_by,
         colour_by = NULL,
+        colour.features = NULL,
         shape.by = shape_by,
         shape_by = NULL,
         size.by = size_by,
@@ -219,6 +224,7 @@ setMethod("plotAbundanceDensity", signature = c(x = "SummarizedExperiment"),
             layout = layout,
             xlab = assay.type,
             colour_by = colour_by,
+            colour_features = colour.features,
             shape_by = shape_by,
             size_by = size_by,
             ...)
@@ -295,6 +301,7 @@ setMethod("plotAbundanceDensity", signature = c(x = "SummarizedExperiment"),
         xlab,
         ylab = NULL,
         colour_by = NULL,
+        colour_features = NULL,
         shape_by = NULL,
         size_by = NULL,
         point_shape = point.shape,
@@ -356,6 +363,9 @@ setMethod("plotAbundanceDensity", signature = c(x = "SummarizedExperiment"),
             size = point_size,
             colour = point_colour)
         point_args$args$mapping$y <- sym("Y")
+        if (!is.null(colour_features)) {
+            point_args$args$fill <- NULL
+        }
         if (layout == "point"){
             plot_out <- plot_out +
                 do.call(geom_point, point_args$args)
@@ -412,6 +422,12 @@ setMethod("plotAbundanceDensity", signature = c(x = "SummarizedExperiment"),
                     axis.line.y = element_blank()) # Removes y-axis
         }
     }
+    # Apply custom feature colours if provided
+    if (!is.null(colour_features)) {
+        plot_out <- plot_out +
+            aes(fill = .data[["Y"]]) +
+            scale_fill_manual(values = colour_features, guide = "none")
+    }
     # add additional guides
     plot_out <- .add_extra_guide(plot_out, shape_by, size_by)
     # add legend
diff --git a/man/plotAbundanceDensity.Rd b/man/plotAbundanceDensity.Rd
index eaa72f8f..fea78c75 100644
--- a/man/plotAbundanceDensity.Rd
+++ b/man/plotAbundanceDensity.Rd
@@ -15,6 +15,7 @@ plotAbundanceDensity(x, ...)
   n = min(nrow(x), 25L),
   colour.by = colour_by,
   colour_by = NULL,
+  colour.features = NULL,
   shape.by = shape_by,
   shape_by = NULL,
   size.by = size_by,
@@ -87,6 +88,10 @@ taxa to show. (Default: \code{min(nrow(x), 25L)})}
 
 \item{colour_by}{Deprecated. Use \code{colour.by} instead.}
 
+\item{colour.features}{\code{Character vector}. A named vector of colors
+for features. Names should match feature names. Useful for matching
+colors across multiple plots. (Default: \code{NULL})}
+
 \item{shape.by}{\code{Character scalar}. Defines a column from
 \code{colData}, that is used to group observations to different point shape
 groups. Must be a value of \code{colData()} function. \code{shape.by} is