diff --git a/R/DRAGON.R b/R/DRAGON.R
index 2cc86774..9538c285 100644
--- a/R/DRAGON.R
+++ b/R/DRAGON.R
@@ -389,6 +389,14 @@ dragon = function(layer1,layer2,pval = FALSE,gradient = "finite_difference", ver
   precmat = get_precision_matrix_dragon(layer1, layer2, lambdas)
   ggm = get_partial_correlation_dragon(layer1, layer2, lambdas)
   
+  # Propagate variable names from input layers to output matrices
+  all_names = c(colnames(layer1), colnames(layer2))
+  if(!is.null(all_names) && length(all_names) == ncol(shrunken_cov)) {
+    rownames(shrunken_cov) = colnames(shrunken_cov) = all_names
+    rownames(precmat) = colnames(precmat) = all_names
+    rownames(ggm) = colnames(ggm) = all_names
+  }
+  
   # if pval, return pval approx with finite difference
   if(pval)
   {
diff --git a/tests/testthat/test-dragon.R b/tests/testthat/test-dragon.R
index 30247508..9aa880a1 100644
--- a/tests/testthat/test-dragon.R
+++ b/tests/testthat/test-dragon.R
@@ -112,6 +112,37 @@ test_that("[DRAGON] dragon() exported function works on simulated data",
 
 })
 
+# test that output matrices carry input variable names (issue #319)
+test_that("[DRAGON] output matrices preserve input column names",
+{
+  toy_layer1 = matrix(c(1,2,3,1,5,12),nrow=3,byrow=T)
+  toy_layer2 = matrix(c(9,7,8),nrow=3,byrow=T)
+  colnames(toy_layer1) = c("gene1","gene2")
+  colnames(toy_layer2) = c("methyl1")
+
+  res = dragon(layer1 = toy_layer1, layer2 = toy_layer2, pval = F)
+  expected_names = c("gene1","gene2","methyl1")
+
+  expect_equal(rownames(res$cov), expected_names)
+  expect_equal(colnames(res$cov), expected_names)
+  expect_equal(rownames(res$prec), expected_names)
+  expect_equal(colnames(res$prec), expected_names)
+  expect_equal(rownames(res$ggm), expected_names)
+  expect_equal(colnames(res$ggm), expected_names)
+})
+
+# test that dragon still works when inputs have no column names
+test_that("[DRAGON] output matrices work without input column names",
+{
+  toy_layer1 = matrix(c(1,2,3,1,5,12),nrow=3,byrow=T)
+  toy_layer2 = matrix(c(9,7,8),nrow=3,byrow=T)
+
+  res = dragon(layer1 = toy_layer1, layer2 = toy_layer2, pval = F)
+
+  expect_null(rownames(res$cov))
+  expect_null(colnames(res$cov))
+})
+
 # test log likelihood function
 # test_that("[DRAGON] Log likelihood function for estimation of kappa is correct",{
 #   # log_lik_shrunken = function(kappa, p, lambda, rhos)