diff --git a/virtualizarr/backend.py b/virtualizarr/backend.py
index 84af17d4..09aa192b 100644
--- a/virtualizarr/backend.py
+++ b/virtualizarr/backend.py
@@ -256,11 +256,16 @@ def open_virtual_dataset(
             )
 
         case "fits":
-            from kerchunk.fits import process_file
-
-            # handle inconsistency in kerchunk, see GH issue https://github.com/zarr-developers/VirtualiZarr/issues/160
-            refs = {"refs": process_file(filepath, **reader_options)}
+            from virtualizarr.readers.fits import open_virtual_dataset
 
+            return open_virtual_dataset(
+                filepath,
+                group=group,
+                drop_variables=drop_variables,
+                loadable_variables=loadable_variables,
+                indexes=indexes,
+                reader_options=reader_options,
+            )
         case _:
             raise NotImplementedError(f"Unsupported file type: {filetype.name}")
 
diff --git a/virtualizarr/readers/fits.py b/virtualizarr/readers/fits.py
new file mode 100644
index 00000000..4a5b2793
--- /dev/null
+++ b/virtualizarr/readers/fits.py
@@ -0,0 +1,53 @@
+from typing import Iterable, Mapping, Optional
+
+from kerchunk.fits import process_file
+from xarray import Dataset
+from xarray.core.indexes import Index
+
+from virtualizarr.backend import (
+    construct_virtual_dataset,
+    open_loadable_vars_and_indexes,
+)
+from virtualizarr.translators.kerchunk import (
+    extract_group,
+    virtual_vars_and_metadata_from_kerchunk_refs,
+)
+
+
+def open_virtual_dataset(
+    filepath: str,
+    group: str | None = None,
+    drop_variables: Iterable[str] | None = None,
+    loadable_variables: Iterable[str] | None = None,
+    decode_times: bool | None = None,
+    indexes: Mapping[str, Index] | None = None,
+    reader_options: Optional[dict] = None,
+) -> Dataset:
+    # handle inconsistency in kerchunk, see GH issue https://github.com/zarr-developers/VirtualiZarr/issues/160
+    refs = {"refs": process_file(filepath, **reader_options)}
+
+    refs = extract_group(refs, group)
+
+    virtual_vars, attrs, coord_names = virtual_vars_and_metadata_from_kerchunk_refs(
+        refs,
+        loadable_variables,
+        drop_variables,
+    )
+
+    loadable_vars, indexes = open_loadable_vars_and_indexes(
+        filepath,
+        loadable_variables=loadable_variables,
+        reader_options=reader_options,
+        drop_variables=drop_variables,
+        indexes=indexes,
+        group=group,
+        decode_times=decode_times,
+    )
+
+    return construct_virtual_dataset(
+        virtual_vars=virtual_vars,
+        loadable_vars=loadable_vars,
+        indexes=indexes,
+        coord_names=coord_names,
+        attrs=attrs,
+    )