diff --git a/DESCRIPTION b/DESCRIPTION
index 698acbd..f0e6ba2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -25,5 +25,9 @@ RoxygenNote: 7.2.2
 Imports:
     colorspace (>= 1.4), parallel, grDevices, graphics,
     stats, utils
-Suggests: knitr, rmarkdown
+Suggests: 
+    knitr,
+    rmarkdown,
+    testthat (>= 3.0.0)
 Maintainer: Giorgio Tumino <giorgio.tumino@wur.nl>
+Config/testthat/edition: 3
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..977a00a
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,12 @@
+# This file is part of the standard setup for testthat.
+# It is recommended that you do not modify it.
+#
+# Where should you do additional test configuration?
+# Learn more about the roles of various files in:
+# * https://r-pkgs.org/testing-design.html#sec-tests-files-overview
+# * https://testthat.r-lib.org/articles/special-files.html
+
+library(testthat)
+library(mpQTL)
+
+test_check("mpQTL")
diff --git a/tests/testthat/test-mpQTL.R b/tests/testthat/test-mpQTL.R
new file mode 100644
index 0000000..286810c
--- /dev/null
+++ b/tests/testthat/test-mpQTL.R
@@ -0,0 +1,13 @@
+
+
+test_that(desc = 'linear model works',
+          code = {expect_no_error(map.QTL(
+                      phenotypes = mppheno,
+                      genotypes = mpsnpdose,
+                      ploidy = 4,
+                      map = mpmap))
+          }
+)
+
+
+
diff --git a/tests/testthat/test-test_z.R b/tests/testthat/test-test_z.R
new file mode 100644
index 0000000..58351ad
--- /dev/null
+++ b/tests/testthat/test-test_z.R
@@ -0,0 +1,64 @@
+
+
+## create test data frame with relication in 2 environments
+testdf <- data.frame(genotype=rep(rownames(mppheno),2),
+                     site = rep(c('S1','S2'),each=nrow(mppheno)),
+                     phenotype1=c(mppheno,mppheno * 1.1))
+## add  row names for tabulation
+rownames(testdf) <- paste(testdf$genotype,testdf$site,sep='_')
+## create pheotype matrix
+pheno_matrix <- matrix(testdf[,3])
+rownames(pheno_matrix) <- rownames(testdf)
+## cofactor matrix for site factor
+cof_matrix <- matrix(rep(c('S1','S2'),each=nrow(mppheno)))
+rownames(cof_matrix) <- rownames(testdf)
+
+## z matrix with dimensions aligned to y matrix and  dosage mmatrix
+z_matrix <- as.matrix(table(rownames(testdf),testdf$genotype))[rownames(testdf),colnames(mpsnpdose)]
+
+test_that(desc='z matrix colnames match snp matrix',
+          code = {expect_identical(colnames(mpsnpdose),colnames(z_matrix))})
+
+test_that(desc='z matrix rownames match phenotype matrix',
+          code = {expect_identical(rownames(pheno_matrix),rownames(z_matrix))})
+
+test_that(desc='cofactor matrix rownames match phenotype matrix',
+          code = {expect_identical(rownames(cof_matrix),rownames(pheno_matrix))})
+
+
+test_that(desc = 'linear model with z matrix  and cofactors works',
+          code = {expect_no_error(
+            map.QTL(
+            phenotypes = pheno_matrix,
+            genotypes = mpsnpdose,
+            Z = z_matrix,
+            cofactor = cof_matrix,
+            cofactor.type = 'cat',
+            ploidy = 4,
+            map = mpmap)
+            )
+          }
+)
+
+## Test Simple Case of a Single Replicate
+##
+mppheno2 <- rbind(54,mppheno)
+rownames(mppheno2) <- c(c("A1_P1_1","A1_P1_2"),rownames(mppheno)[-1])
+# form diag matrix
+Z <- diag(nrow(mppheno))
+z_matrix2 <- rbind(Z[1,],Z[1,],Z[-1,])
+
+rownames(z_matrix2) <- rownames(mppheno2)
+colnames(z_matrix2) <- rownames(mppheno)
+
+test_that(desc = 'linear model with z matrix  no cofactors works',
+          code = {expect_no_error(
+            map.QTL(
+              phenotypes = mppheno2,
+              genotypes = mpsnpdose,
+              Z = z_matrix2,
+              ploidy = 4,
+              map = mpmap)
+          )
+          }
+)