diff --git a/package/MDAnalysis/analysis/contacts.py b/package/MDAnalysis/analysis/contacts.py
index 50ea47f4e90..a39ad96138d 100644
--- a/package/MDAnalysis/analysis/contacts.py
+++ b/package/MDAnalysis/analysis/contacts.py
@@ -339,9 +339,17 @@ def __init__(self, u, selection, refgroup, method="cutoff", radius=4.5,
         self.grB = u.select_atoms(selection[1])
 
         # contacts formed in reference
-        refA, refB = refgroup
-        self.r0 = distance_array(refA.positions, refB.positions)
-        self.initial_contacts = contact_matrix(self.r0, radius)
+        self.r0 = []
+        self.initial_contacts = []
+
+        if isinstance(refgroup[0], mda.AtomGroup.AtomGroup):
+            refA, refB = refgroup
+            self.r0.append(distance_array(refA.positions, refB.positions))
+            self.initial_contacts.append(contact_matrix(self.r0, radius))
+        else:
+            for refA, refB in refgroup:
+                self.r0.append(distance_array(refA.positions, refB.positions))
+                self.initial_contacts.append(contact_matrix(self.r0, radius))
 
         self.fraction_kwargs = kwargs
         self.timeseries = []
@@ -350,11 +358,13 @@ def _single_frame(self):
         # compute distance array for a frame
         d = distance_array(self.grA.positions, self.grB.positions)
 
+        y = np.empty(len(self.r0))
+        for i, (initial_contacts, r0) in enumerate(zip(self.initial_contacts,
+                                                       self.r0)):
         # select only the contacts that were formed in the reference state
-        r = d[self.initial_contacts]
-        r0 = self.r0[self.initial_contacts]
-
-        y = self.fraction_contacts(r, r0, **self.fraction_kwargs)
+            r = d[initial_contacts]
+            r0 = r0[initial_contacts]
+            y[i] = self.fraction_contacts(r, r0, **self.fraction_kwargs)
 
         self.timeseries.append((self._ts.frame, y))