diff --git a/holoviews/operation/element.py b/holoviews/operation/element.py
index 6b35b4432e..5bb20b10ae 100644
--- a/holoviews/operation/element.py
+++ b/holoviews/operation/element.py
@@ -519,7 +519,7 @@ class histogram(Operation):
     dimension = param.String(default=None, doc="""
       Along which dimension of the Element to compute the histogram.""")
 
-    frequency_label = param.String(default='{dim} Frequency', doc="""
+    frequency_label = param.String(default=None, doc="""
       Format string defining the label of the frequency dimension of the Histogram.""")
 
     groupby = param.ClassSelector(default=None, class_=(basestring, Dimension), doc="""
@@ -552,31 +552,32 @@ class histogram(Operation):
     style_prefix = param.String(default=None, allow_None=None, doc="""
       Used for setting a common style for histograms in a HoloMap or AdjointLayout.""")
 
-    def _process(self, view, key=None):
+    def _process(self, element, key=None):
         if self.p.groupby:
-            if not isinstance(view, Dataset):
+            if not isinstance(element, Dataset):
                 raise ValueError('Cannot use histogram groupby on non-Dataset Element')
-            grouped = view.groupby(self.p.groupby, group_type=Dataset, container_type=NdOverlay)
+            grouped = element.groupby(self.p.groupby, group_type=Dataset, container_type=NdOverlay)
             self.p.groupby = None
             return grouped.map(self._process, Dataset)
 
         if self.p.dimension:
             selected_dim = self.p.dimension
         else:
-            selected_dim = [d.name for d in view.vdims + view.kdims][0]
-        data = np.array(view.dimension_values(selected_dim))
+            selected_dim = [d.name for d in element.vdims + element.kdims][0]
+        dim = element.get_dimension(selected_dim)
+        data = np.array(element.dimension_values(selected_dim))
         if self.p.nonzero:
             mask = data > 0
             data = data[mask]
         if self.p.weight_dimension:
-            weights = np.array(view.dimension_values(self.p.weight_dimension))
+            weights = np.array(element.dimension_values(self.p.weight_dimension))
             if self.p.nonzero:
                 weights = weights[mask]
         else:
             weights = None
 
         data = data[isfinite(data)]
-        hist_range = self.p.bin_range or view.range(selected_dim)
+        hist_range = self.p.bin_range or element.range(selected_dim)
         # Avoids range issues including zero bin range and empty bins
         if hist_range == (0, 0) or any(not isfinite(r) for r in hist_range):
             hist_range = (0, 1)
@@ -625,27 +626,30 @@ def _process(self, view, key=None):
 
         params = {}
         if self.p.weight_dimension:
-            params['vdims'] = [view.get_dimension(self.p.weight_dimension)]
+            params['vdims'] = [element.get_dimension(self.p.weight_dimension)]
+        elif self.p.frequency_label:
+            label = self.p.frequency_label.format(dim=dim.pprint_label)
+            params['vdims'] = [Dimension('Frequency', label=label)]
         else:
-            label = self.p.frequency_label.format(dim=selected_dim)
-            params['vdims'] = [Dimension('{}_frequency'.format(selected_dim),
+            label = 'Frequency' if normed else 'Count'
+            params['vdims'] = [Dimension('{0}_{1}'.format(dim.name, label.lower()),
                                          label=label)]
 
-        if view.group != view.__class__.__name__:
-            params['group'] = view.group
+        if element.group != element.__class__.__name__:
+            params['group'] = element.group
 
         if self.p.cumulative:
             hist = np.cumsum(hist)
             if self.p.normed in (True, 'integral'):
                 hist *= edges[1]-edges[0]
-        return Histogram((edges, hist), kdims=[view.get_dimension(selected_dim)],
-                         label=view.label, **params)
+        return Histogram((edges, hist), kdims=[element.get_dimension(selected_dim)],
+                         label=element.label, **params)
 
 
 class decimate(Operation):
     """
     Decimates any column based Element to a specified number of random
-    rows if the current view defined by the x_range and y_range
+    rows if the current element defined by the x_range and y_range
     contains more than max_samples. By default the operation returns a
     DynamicMap with a RangeXY stream allowing dynamic downsampling.
     """
diff --git a/holoviews/operation/stats.py b/holoviews/operation/stats.py
index 54d1719324..51373a69bb 100644
--- a/holoviews/operation/stats.py
+++ b/holoviews/operation/stats.py
@@ -78,8 +78,7 @@ def _process(self, element, key=None):
             params['label'] = element.label
             vdim = element.vdims[0]
             vdim_name = '{}_density'.format(selected_dim.name)
-            vdim_label = '{} Density'.format(selected_dim.label)
-            vdims = [vdim(vdim_name, label=vdim_label) if vdim.name == 'Density' else vdim]
+            vdims = [vdim(vdim_name, label='Density') if vdim.name == 'Density' else vdim]
         else:
             if self.p.dimension:
                 selected_dim = element.get_dimension(self.p.dimension)
@@ -91,8 +90,7 @@ def _process(self, element, key=None):
                                      type(element).__name__)
                 selected_dim = dimensions[0]
             vdim_name = '{}_density'.format(selected_dim.name)
-            vdim_label = '{} Density'.format(selected_dim.label)
-            vdims = [Dimension(vdim_name, label=vdim_label)]
+            vdims = [Dimension(vdim_name, label='Density')]
 
         data = element.dimension_values(selected_dim)
         bin_range = self.p.bin_range or element.range(selected_dim)
diff --git a/holoviews/plotting/bokeh/chart.py b/holoviews/plotting/bokeh/chart.py
index e99bacf10c..7fb2587edb 100644
--- a/holoviews/plotting/bokeh/chart.py
+++ b/holoviews/plotting/bokeh/chart.py
@@ -14,7 +14,7 @@
 from ...element import Bars
 from ...operation import interpolate_curve
 from ...util.transform import dim
-from ..util import compute_sizes, get_min_distance, dim_axis_label, get_axis_padding
+from ..util import compute_sizes, get_min_distance, get_axis_padding
 from .element import ElementPlot, ColorbarPlot, LegendPlot
 from .styles import (expand_batched_style, line_properties, fill_properties,
                      mpl_to_bokeh, rgb2hex)
@@ -824,22 +824,13 @@ def _get_factors(self, element):
         return coords
 
 
-    def _get_axis_labels(self, *args, **kwargs):
-        """
-        Override axis mapping by setting the first key and value
-        dimension as the x-axis and y-axis labels.
-        """
-        element = self.current_frame
-        if self.batched:
-            element = element.last
-        xlabel = dim_axis_label(element.kdims[0])
+    def _get_axis_dims(self, element):
         if element.ndims > 1 and not (self.stacked or self.stack_index):
-            gdim = element.get_dimension(1)
+            xdims = element.kdims
         else:
-            gdim = None
-        if gdim and gdim in element.kdims:
-            xlabel = ', '.join([xlabel, dim_axis_label(gdim)])
-        return (xlabel, dim_axis_label(element.vdims[0]), None)
+            xdims = element.kdims[0]
+        return (xdims, element.vdims[0])
+
 
     def get_stack(self, xvals, yvals, baselines, sign='positive'):
         """
diff --git a/holoviews/plotting/bokeh/element.py b/holoviews/plotting/bokeh/element.py
index 3e60d1d4d0..b0e3364771 100644
--- a/holoviews/plotting/bokeh/element.py
+++ b/holoviews/plotting/bokeh/element.py
@@ -40,7 +40,7 @@
 from .util import (
     bokeh_version, decode_bytes, get_tab_title, glyph_order,
     py2js_tickformatter, recursive_model_update, theme_attr_json,
-    cds_column_replace, hold_policy
+    cds_column_replace, hold_policy, match_dim_specs
 )
 
 
@@ -217,43 +217,69 @@ def _get_hover_data(self, data, element, dimensions=None):
                 data[dim] = [v for _ in range(len(list(data.values())[0]))]
 
 
-    def _merge_ranges(self, plots, xlabel, ylabel):
+    def _merge_ranges(self, plots, xspecs, yspecs):
         """
         Given a list of other plots return axes that are shared
-        with another plot by matching the axes labels
+        with another plot by matching the dimensions specs stored
+        as tags on the dimensions.
         """
         plot_ranges = {}
         for plot in plots:
             if plot is None:
                 continue
-            if hasattr(plot, 'xaxis'):
-                if plot.xaxis[0].axis_label == xlabel:
+            if hasattr(plot, 'x_range') and plot.x_range.tags and xspecs is not None:
+                if match_dim_specs(plot.x_range.tags[0], xspecs):
                     plot_ranges['x_range'] = plot.x_range
-                if plot.xaxis[0].axis_label == ylabel:
+                if match_dim_specs(plot.x_range.tags[0], yspecs):
                     plot_ranges['y_range'] = plot.x_range
-            if hasattr(plot, 'yaxis'):
-                if plot.yaxis[0].axis_label == ylabel:
+            if hasattr(plot, 'y_range') and plot.y_range.tags and yspecs is not None:
+                if match_dim_specs(plot.y_range.tags[0], yspecs):
                     plot_ranges['y_range'] = plot.y_range
-                if plot.yaxis[0].axis_label == xlabel:
+                if match_dim_specs(plot.y_range.tags[0], xspecs):
                     plot_ranges['x_range'] = plot.y_range
         return plot_ranges
 
 
+    def _get_axis_dims(self, element):
+        """Returns the dimensions corresponding to each axis.
+
+        Should return a list of dimensions or list of lists of
+        dimensions, which will be formatted to label the axis
+        and to link axes.
+        """
+        dims = element.dimensions()[:2]
+        return dims+[None] if len(dims) < 2 else dims
+
+
     def _axes_props(self, plots, subplots, element, ranges):
         # Get the bottom layer and range element
         el = element.traverse(lambda x: x, [Element])
         el = el[0] if el else element
 
-        dims = el.nodes.dimensions() if isinstance(el, Graph) else el.dimensions()
+        dims = self._get_axis_dims(el)
         xlabel, ylabel, zlabel = self._get_axis_labels(dims)
         if self.invert_axes:
             xlabel, ylabel = ylabel, xlabel
+            dims = dims[:2][::-1]
+        xdims, ydims = dims[:2]
+        if xdims:
+            if not isinstance(xdims, list):
+                xdims = [xdims]
+            xspecs = tuple((xd.name, xd.label, xd.unit) for xd in xdims)
+        else:
+            xspecs = None
+        if ydims:
+            if not isinstance(ydims, list):
+                ydims = [ydims]
+            yspecs = tuple((yd.name, yd.label, yd.unit) for yd in ydims)
+        else:
+            yspecs = None
 
         plot_ranges = {}
         # Try finding shared ranges in other plots in the same Layout
         norm_opts = self.lookup_options(el, 'norm').options
         if plots and self.shared_axes and not norm_opts.get('axiswise', False):
-            plot_ranges = self._merge_ranges(plots, xlabel, ylabel)
+            plot_ranges = self._merge_ranges(plots, xspecs, yspecs)
 
         # Get the Element that determines the range and get_extents
         range_el = el if self.batched and not isinstance(self, OverlayPlot) else element
@@ -261,19 +287,37 @@ def _axes_props(self, plots, subplots, element, ranges):
         if self.invert_axes:
             l, b, r, t = b, l, t, r
 
-        xtype = el.get_dimension_type(0)
-        if ((xtype is np.object_ and type(l) in util.datetime_types) or
-            xtype in util.datetime_types):
-            x_axis_type = 'datetime'
-        else:
-            x_axis_type = 'log' if self.logx else 'auto'
+        categorical = any(self.traverse(lambda x: x._categorical))
+        categorical_x = any(isinstance(x, util.basestring) for x in (l, r))
+        categorical_y = any(isinstance(y, util.basestring) for y in (b, t))
+
+        x_axis_type = 'log' if self.logx else 'auto'
+        if xdims:
+            if len(xdims) > 1:
+                x_axis_type = 'auto'
+                categorical_x = True
+            else:
+                if isinstance(el, Graph):
+                    xtype = el.nodes.get_dimension_type(xdims[0])
+                else:
+                    xtype = el.get_dimension_type(xdims[0])
+                if ((xtype is np.object_ and type(l) in util.datetime_types) or
+                    xtype in util.datetime_types):
+                    x_axis_type = 'datetime'
 
         y_axis_type = 'log' if self.logy else 'auto'
-        if len(dims) > 1:
-            ytype = el.get_dimension_type(1)
-            if ((ytype is np.object_ and type(b) in util.datetime_types)
-                or ytype in util.datetime_types):
-                y_axis_type = 'datetime'
+        if ydims:
+            if len(ydims) > 1:
+                y_axis_type = 'auto'
+                categorical_y = True
+            else:
+                if isinstance(el, Graph):
+                    ytype = el.nodes.get_dimension_type(ydims[0])
+                else:
+                    ytype = el.get_dimension_type(ydims[0])
+                if ((ytype is np.object_ and type(b) in util.datetime_types)
+                    or ytype in util.datetime_types):
+                    y_axis_type = 'datetime'
 
         # Declare shared axes
         if 'x_range' in plot_ranges:
@@ -281,10 +325,6 @@ def _axes_props(self, plots, subplots, element, ranges):
         if 'y_range' in plot_ranges:
             self._shared['y'] = True
 
-        categorical = any(self.traverse(lambda x: x._categorical))
-        categorical_x = any(isinstance(x, util.basestring) for x in (l, r))
-        categorical_y = any(isinstance(y, util.basestring) for y in (b, t))
-
         range_types = (self._x_range_type, self._y_range_type)
         if self.invert_axes: range_types = range_types[::-1]
         x_range_type, y_range_type = range_types
@@ -300,6 +340,12 @@ def _axes_props(self, plots, subplots, element, ranges):
         elif 'y_range' not in plot_ranges:
             plot_ranges['y_range'] = y_range_type()
 
+        x_range, y_range = plot_ranges['x_range'], plot_ranges['y_range']
+        if not x_range.tags and xspecs is not None:
+            x_range.tags.append(xspecs)
+        if not y_range.tags and yspecs is not None:
+            y_range.tags.append(yspecs)
+
         return (x_axis_type, y_axis_type), (xlabel, ylabel, zlabel), plot_ranges
 
 
@@ -452,7 +498,7 @@ def _axis_properties(self, axis, key, plot, dimension=None,
                 formatter = PrintfTickFormatter(format=formatter)
             if formatter is not None:
                 axis_props['formatter'] = formatter
-        elif FuncTickFormatter is not None and ax_mapping and dimension:
+        elif FuncTickFormatter is not None and ax_mapping and isinstance(dimension, Dimension):
             formatter = None
             if dimension.value_format:
                 formatter = dimension.value_format
@@ -495,15 +541,14 @@ def _update_plot(self, key, plot, element=None):
         """
         el = element.traverse(lambda x: x, [Element])
         el = el[0] if el else element
-        dimensions = el.nodes.dimensions() if isinstance(el, Graph) else el.dimensions()
-        if not len(dimensions) >= 2:
-            dimensions = dimensions+[None]
+        dimensions = self._get_axis_dims(el)
         plot.update(**self._plot_properties(key, plot, element))
 
         props = {axis: self._axis_properties(axis, key, plot, dim)
                  for axis, dim in zip(['x', 'y'], dimensions)}
         xlabel, ylabel, zlabel = self._get_axis_labels(dimensions)
-        if self.invert_axes: xlabel, ylabel = ylabel, xlabel
+        if self.invert_axes:
+            xlabel, ylabel = ylabel, xlabel
         props['x']['axis_label'] = xlabel if 'x' in self.labelled else ''
         props['y']['axis_label'] = ylabel if 'y' in self.labelled else ''
         recursive_model_update(plot.xaxis[0], props.get('x', {}))
@@ -1695,6 +1740,13 @@ def _get_factors(self, overlay):
         return util.unique_array(xfactors), util.unique_array(yfactors)
 
 
+    def _get_axis_dims(self, element):
+        subplots = list(self.subplots.values())
+        if subplots:
+            return subplots[0]._get_axis_dims(element)
+        return super(OverlayPlot, self)._get_axis_dims(element)
+
+
     def initialize_plot(self, ranges=None, plot=None, plots=None):
         key = util.wrap_tuple(self.hmap.last_key)
         nonempty = [(k, el) for k, el in self.hmap.data.items() if el]
diff --git a/holoviews/plotting/bokeh/graphs.py b/holoviews/plotting/bokeh/graphs.py
index 2d61a7b968..3094b06867 100644
--- a/holoviews/plotting/bokeh/graphs.py
+++ b/holoviews/plotting/bokeh/graphs.py
@@ -94,13 +94,8 @@ def get_extents(self, element, ranges, range_type='combined'):
         return super(GraphPlot, self).get_extents(element.nodes, ranges, range_type)
 
 
-    def _get_axis_labels(self, *args, **kwargs):
-        """
-        Override axis labels to group all key dimensions together.
-        """
-        element = self.current_frame
-        xlabel, ylabel = [kd.pprint_label for kd in element.nodes.kdims[:2]]
-        return xlabel, ylabel, None
+    def _get_axis_dims(self, element):
+        return element.nodes.dimensions()[:2]
 
 
     def _get_edge_colors(self, element, ranges, edge_data, edge_mapping, style):
diff --git a/holoviews/plotting/bokeh/stats.py b/holoviews/plotting/bokeh/stats.py
index 8e6327ba61..432a7cbe2b 100644
--- a/holoviews/plotting/bokeh/stats.py
+++ b/holoviews/plotting/bokeh/stats.py
@@ -90,14 +90,8 @@ def get_extents(self, element, ranges, range_type='combined'):
             element, ranges, range_type, 'categorical', element.vdims[0]
         )
 
-    def _get_axis_labels(self, *args, **kwargs):
-        """
-        Override axis labels to group all key dimensions together.
-        """
-        element = self.current_frame
-        xlabel = ', '.join([kd.pprint_label for kd in element.kdims])
-        ylabel = element.vdims[0].pprint_label
-        return xlabel, ylabel, None
+    def _get_axis_dims(self, element):
+        return element.kdims, element.vdims[0]
 
     def _glyph_properties(self, plot, element, source, ranges, style, group=None):
         if element.ndims > 0:
diff --git a/holoviews/plotting/bokeh/util.py b/holoviews/plotting/bokeh/util.py
index 8fe0ffd978..145bd9a75e 100644
--- a/holoviews/plotting/bokeh/util.py
+++ b/holoviews/plotting/bokeh/util.py
@@ -620,3 +620,22 @@ def multi_polygons_data(element):
         xsh.append(multi_xs)
         ysh.append(multi_ys)
     return xsh, ysh
+
+
+def match_dim_specs(specs1, specs2):
+    """Matches dimension specs used to link axes.
+
+    Axis dimension specs consists of a list of tuples corresponding
+    to each dimension, each tuple spec has the form (name, label, unit).
+    The name and label must match exactly while the unit only has to
+    match if both specs define one.
+    """
+    if (specs1 is None or specs2 is None) or (len(specs1) != len(specs2)):
+        return False
+    for spec1, spec2 in zip(specs1, specs2):
+        for s1, s2 in zip(spec1, spec2):
+            if s1 is None or s2 is None:
+                continue
+            if s1 != s2:
+                return False
+    return True
diff --git a/holoviews/plotting/mpl/raster.py b/holoviews/plotting/mpl/raster.py
index 701173be8d..a41dc0508a 100644
--- a/holoviews/plotting/mpl/raster.py
+++ b/holoviews/plotting/mpl/raster.py
@@ -330,8 +330,7 @@ def _get_axis_kwargs(self):
         xticks = (self._xticks, [xdim.pprint_value(l) for l in self._xkeys])
         yticks = (self._yticks, [ydim.pprint_value(l) if ydim else ''
                                  for l in self._ykeys])
-        return dict(xlabel=xdim.pprint_label, ylabel=ydim.pprint_label if ydim else '',
-                    xticks=xticks, yticks=yticks)
+        return dict(dimensions=[xdim, ydim], xticks=xticks, yticks=yticks)
 
 
     def _compute_borders(self):
diff --git a/holoviews/plotting/plot.py b/holoviews/plotting/plot.py
index 2afd214cb0..02c0cf150d 100644
--- a/holoviews/plotting/plot.py
+++ b/holoviews/plotting/plot.py
@@ -1042,11 +1042,15 @@ def _get_axis_labels(self, dimensions, xlabel=None, ylabel=None, zlabel=None):
         if self.xlabel is not None:
             xlabel = self.xlabel
         elif dimensions and xlabel is None:
-            xlabel = dim_axis_label(dimensions[0]) if dimensions[0] else ''
+            xdims = dimensions[0]
+            xlabel = dim_axis_label(xdims) if xdims else ''
+
         if self.ylabel is not None:
             ylabel = self.ylabel
         elif len(dimensions) >= 2 and ylabel is None:
-            ylabel = dim_axis_label(dimensions[1]) if dimensions[1] else ''
+            ydims = dimensions[1]
+            ylabel = dim_axis_label(ydims) if ydims else ''
+
         if getattr(self, 'zlabel', None) is not None:
             zlabel = self.zlabel
         elif self.projection == '3d' and len(dimensions) >= 3 and zlabel is None:
diff --git a/holoviews/tests/element/teststatselements.py b/holoviews/tests/element/teststatselements.py
index a2670a49fa..e0208e0017 100644
--- a/holoviews/tests/element/teststatselements.py
+++ b/holoviews/tests/element/teststatselements.py
@@ -122,7 +122,7 @@ def test_distribution_composite(self):
         dist = Distribution(np.array([0, 1, 2]))
         area = Compositor.collapse_element(dist, backend='matplotlib')
         self.assertIsInstance(area, Area)
-        self.assertEqual(area.vdims, [Dimension(('Value_density', 'Value Density'))])
+        self.assertEqual(area.vdims, [Dimension(('Value_density', 'Density'))])
 
     def test_distribution_composite_transfer_opts(self):
         dist = Distribution(np.array([0, 1, 2])).opts(style=dict(color='red'))
@@ -146,13 +146,13 @@ def test_distribution_composite_not_filled(self):
         dist = Distribution(np.array([0, 1, 2]), ).opts(plot=dict(filled=False))
         curve = Compositor.collapse_element(dist, backend='matplotlib')
         self.assertIsInstance(curve, Curve)
-        self.assertEqual(curve.vdims, [Dimension(('Value_density', 'Value Density'))])
+        self.assertEqual(curve.vdims, [Dimension(('Value_density', 'Density'))])
 
     def test_distribution_composite_empty_not_filled(self):
         dist = Distribution([]).opts(plot=dict(filled=False))
         curve = Compositor.collapse_element(dist, backend='matplotlib')
         self.assertIsInstance(curve, Curve)
-        self.assertEqual(curve.vdims, [Dimension(('Value_density', 'Value Density'))])
+        self.assertEqual(curve.vdims, [Dimension(('Value_density', 'Density'))])
 
     def test_bivariate_composite(self):
         dist = Bivariate(np.random.rand(10, 2))
diff --git a/holoviews/tests/operation/testoperation.py b/holoviews/tests/operation/testoperation.py
index fcdcf0801e..1f1ed744ca 100644
--- a/holoviews/tests/operation/testoperation.py
+++ b/holoviews/tests/operation/testoperation.py
@@ -128,42 +128,24 @@ def test_points_histogram(self):
         points = Points([float(i) for i in range(10)])
         op_hist = histogram(points, num_bins=3)
 
-        # Make sure that the name and label are as desired
-        op_freq_dim = op_hist.get_dimension('x_frequency')
-        self.assertEqual(op_freq_dim.label, 'x Frequency')
-
-        # Because the operation labels are now different from the
-        #  default Element label, change back before comparing.
-        op_hist = op_hist.redim(x_frequency='Frequency')
-        hist = Histogram(([0.1, 0.1, 0.133333], [0, 3, 6, 9]))
+        hist = Histogram(([0.1, 0.1, 0.133333], [0, 3, 6, 9]),
+                         vdims=('x_frequency', 'Frequency'))
         self.assertEqual(op_hist, hist)
 
     def test_points_histogram_bin_range(self):
         points = Points([float(i) for i in range(10)])
         op_hist = histogram(points, num_bins=3, bin_range=(0, 3))
 
-        # Make sure that the name and label are as desired
-        op_freq_dim = op_hist.get_dimension('x_frequency')
-        self.assertEqual(op_freq_dim.label, 'x Frequency')
-
-        # Because the operation labels are now different from the
-        #  default Element label, change back before comparing.
-        op_hist = op_hist.redim(x_frequency='Frequency')
-        hist = Histogram(([0.25, 0.25, 0.5], [0., 1., 2., 3.]))
+        hist = Histogram(([0.25, 0.25, 0.5], [0., 1., 2., 3.]),
+                         vdims=('x_frequency', 'Frequency'))
         self.assertEqual(op_hist, hist)
 
     def test_points_histogram_explicit_bins(self):
         points = Points([float(i) for i in range(10)])
         op_hist = histogram(points, bins=[0, 1, 3], normed=False)
 
-        # Make sure that the name and label are as desired
-        op_freq_dim = op_hist.get_dimension('x_frequency')
-        self.assertEqual(op_freq_dim.label, 'x Frequency')
-
-        # Because the operation labels are now different from the
-        #  default Element label, change back before comparing.
-        op_hist = op_hist.redim(x_frequency='Frequency')
-        hist = Histogram(([0, 1, 3], [1, 3]))
+        hist = Histogram(([0, 1, 3], [1, 3]),
+                         vdims=('x_count', 'Count'))
         self.assertEqual(op_hist, hist)
 
     def test_points_histogram_cumulative(self):
@@ -171,28 +153,16 @@ def test_points_histogram_cumulative(self):
         points = Points(arr)
         op_hist = histogram(points, cumulative=True, num_bins=3, normed=False)
 
-        # Make sure that the name and label are as desired
-        op_freq_dim = op_hist.get_dimension('x_frequency')
-        self.assertEqual(op_freq_dim.label, 'x Frequency')
-
-        # Because the operation labels are now different from the
-        #  default Element label, change back before comparing.
-        op_hist = op_hist.redim(x_frequency='Frequency')
-        hist = Histogram(([0, 1, 2, 3], [1, 2, 4]))
+        hist = Histogram(([0, 1, 2, 3], [1, 2, 4]),
+                         vdims=('x_count', 'Count'))
         self.assertEqual(op_hist, hist)
 
     def test_points_histogram_not_normed(self):
         points = Points([float(i) for i in range(10)])
         op_hist = histogram(points, num_bins=3, normed=False)
 
-        # Make sure that the name and label are as desired
-        op_freq_dim = op_hist.get_dimension('x_frequency')
-        self.assertEqual(op_freq_dim.label, 'x Frequency')
-
-        # Because the operation labels are now different from the
-        #  default Element label, change back before comparing.
-        op_hist = op_hist.redim(x_frequency='Frequency')
-        hist = Histogram(([3, 3, 4], [0, 3, 6, 9]))
+        hist = Histogram(([3, 3, 4], [0, 3, 6, 9]),
+                         vdims=('x_count', 'Count'))
         self.assertEqual(op_hist, hist)
 
     def test_histogram_operation_datetime(self):
@@ -203,7 +173,7 @@ def test_histogram_operation_datetime(self):
                                        '2017-01-04T00:00:00.000000'], dtype='datetime64[us]'),
                      'Date_frequency': np.array([  3.85802469e-18,   3.85802469e-18,   3.85802469e-18,
                                                    3.85802469e-18])}
-        hist = Histogram(hist_data, kdims='Date', vdims=('Date_frequency', 'Date Frequency'))
+        hist = Histogram(hist_data, kdims='Date', vdims=('Date_frequency', 'Frequency'))
         self.assertEqual(op_hist, hist)
 
     def test_histogram_operation_datetime64(self):
@@ -214,7 +184,7 @@ def test_histogram_operation_datetime64(self):
                                        '2017-01-04T00:00:00.000000'], dtype='datetime64[us]'),
                      'Date_frequency': np.array([  3.85802469e-18,   3.85802469e-18,   3.85802469e-18,
                                                    3.85802469e-18])}
-        hist = Histogram(hist_data, kdims='Date', vdims=('Date_frequency', 'Date Frequency'))
+        hist = Histogram(hist_data, kdims='Date', vdims=('Date_frequency', 'Frequency'))
         self.assertEqual(op_hist, hist)
 
     @attr(optional=1) # Requires matplotlib
@@ -226,7 +196,7 @@ def test_histogram_operation_pd_period(self):
                                        '2017-01-04T00:00:00.000000'], dtype='datetime64[us]'),
                      'Date_frequency': np.array([  3.85802469e-18,   3.85802469e-18,   3.85802469e-18,
                                                    3.85802469e-18])}
-        hist = Histogram(hist_data, kdims='Date', vdims=('Date_frequency', 'Date Frequency'))
+        hist = Histogram(hist_data, kdims='Date', vdims=('Date_frequency', 'Frequency'))
         self.assertEqual(op_hist, hist)
 
     def test_points_histogram_weighted(self):
diff --git a/holoviews/tests/operation/teststatsoperations.py b/holoviews/tests/operation/teststatsoperations.py
index e5aa580fe2..0237917938 100644
--- a/holoviews/tests/operation/teststatsoperations.py
+++ b/holoviews/tests/operation/teststatsoperations.py
@@ -32,20 +32,20 @@ def test_univariate_kde(self):
         kde = univariate_kde(self.dist, n_samples=5, bin_range=(0, 4))
         xs = np.arange(5)
         ys = [0.17594505, 0.23548218, 0.23548218, 0.17594505, 0.0740306]
-        area = Area((xs, ys), 'Value', ('Value_density', 'Value Density'))
+        area = Area((xs, ys), 'Value', ('Value_density', 'Density'))
         self.assertEqual(kde, area)
 
     def test_univariate_kde_flat_distribution(self):
         dist = Distribution([1, 1, 1])
         kde = univariate_kde(dist, n_samples=5, bin_range=(0, 4))
-        area = Area([], 'Value', ('Value_density', 'Value Density'))
+        area = Area([], 'Value', ('Value_density', 'Density'))
         self.assertEqual(kde, area)
 
     def test_univariate_kde_nans(self):
         kde = univariate_kde(self.dist_nans, n_samples=5, bin_range=(0, 4))
         xs = np.arange(5)
         ys = [0, 0, 0, 0, 0]
-        area = Area((xs, ys), 'Value', ('Value_density', 'Value Density'))
+        area = Area((xs, ys), 'Value', ('Value_density', 'Density'))
         self.assertEqual(kde, area)
 
     def test_bivariate_kde(self):
diff --git a/holoviews/tests/plotting/bokeh/testlayoutplot.py b/holoviews/tests/plotting/bokeh/testlayoutplot.py
index 75754cf6d8..70a263e99d 100644
--- a/holoviews/tests/plotting/bokeh/testlayoutplot.py
+++ b/holoviews/tests/plotting/bokeh/testlayoutplot.py
@@ -1,7 +1,7 @@
 import numpy as np
 
 from holoviews.core import (HoloMap, GridSpace, Layout, Empty, Dataset,
-                            NdOverlay, DynamicMap)
+                            NdOverlay, DynamicMap, Dimension)
 from holoviews.element import Curve, Image, Points
 from holoviews.streams import Stream
 
@@ -259,3 +259,29 @@ def test_layout_dimensioned_stream_title_update(self):
         self.assertIn('test: 1', plot.handles['title'].text)
         plot.cleanup()
         self.assertEqual(stream._subscribers, [])
+
+    def test_layout_axis_link_matching_name_label(self):
+        layout = Curve([1, 2, 3], vdims=('a', 'A')) + Curve([1, 2, 3], vdims=('a', 'A'))
+        plot = bokeh_renderer.get_plot(layout)
+        p1, p2 = (sp.subplots['main'] for sp in plot.subplots.values())
+        self.assertIs(p1.handles['y_range'], p2.handles['y_range'])
+
+    def test_layout_axis_not_linked_mismatching_name(self):
+        layout = Curve([1, 2, 3], vdims=('b', 'A')) + Curve([1, 2, 3], vdims=('a', 'A'))
+        plot = bokeh_renderer.get_plot(layout)
+        p1, p2 = (sp.subplots['main'] for sp in plot.subplots.values())
+        self.assertIsNot(p1.handles['y_range'], p2.handles['y_range'])
+
+    def test_layout_axis_linked_unit_and_no_unit(self):
+        layout = (Curve([1, 2, 3], vdims=Dimension('length', unit='m')) +
+                  Curve([1, 2, 3], vdims='length'))
+        plot = bokeh_renderer.get_plot(layout)
+        p1, p2 = (sp.subplots['main'] for sp in plot.subplots.values())
+        self.assertIs(p1.handles['y_range'], p2.handles['y_range'])
+
+    def test_layout_axis_not_linked_mismatching_unit(self):
+        layout = (Curve([1, 2, 3], vdims=Dimension('length', unit='m')) +
+                  Curve([1, 2, 3], vdims=Dimension('length', unit='cm')))
+        plot = bokeh_renderer.get_plot(layout)
+        p1, p2 = (sp.subplots['main'] for sp in plot.subplots.values())
+        self.assertIsNot(p1.handles['y_range'], p2.handles['y_range'])